From 14bab8de4caea55683c76ad6f0722f17638f92b4 Mon Sep 17 00:00:00 2001
From: Ferdinand Fischer <ferdinand.fischer@fau.de>
Date: Thu, 23 May 2019 17:59:15 +0200
Subject: [PATCH] Some changes: *) Allow to set empty quantity.Discrete objects
with certain size *) added the isConstant property for the quantity.Discrete
*) Implemented a cast to the quantity.Operator object *) Implemented a cast
of empty quantity.Discretes to quantity.Symbolics *) Fixed bug in "on()" for
constant quantities *) Fixed some functions to work with empty quantities *)
Added a settings object to use global settings *) plot function returns now
an array of the figure handles *) Re-implementation of the quantity.Operator
object
---
+quantity/BasicVariable.m | 394 +++++++++++++-------------
+quantity/Discrete.m | 130 +++++++--
+quantity/Operator.m | 412 ++++++++++++++++++++--------
+quantity/Settings.m | 27 ++
+quantity/Symbolic.m | 43 +--
+unittests/+quantity/testDiscrete.m | 11 +
+unittests/+quantity/testOperator.m | 275 ++++++++++---------
+unittests/+quantity/testSymbolic.m | 6 +-
8 files changed, 818 insertions(+), 480 deletions(-)
create mode 100644 +quantity/Settings.m
diff --git a/+quantity/BasicVariable.m b/+quantity/BasicVariable.m
index 30cb439..8069312 100644
--- a/+quantity/BasicVariable.m
+++ b/+quantity/BasicVariable.m
@@ -1,199 +1,199 @@
classdef BasicVariable < quantity.Function
- %GEVREYFUNCTION class for desription of basic variables
- % The basic variable can be of the form
- % g(t) = [phi_1(t); phi(2); ... ]
- %
- properties
- % The constant gain of the Gevrey-function
- K double = 1;
- % Shifts the evaluation of the derivative of the Gevrey-function.
- % Can be used if the series, computed using the gevrey function,
- % should not start at the 0-th derivative.
- diffShift double {mustBeNonnegative, mustBeInteger};
- % Offset to raise the whole Gevrey-function g_ = g + offset.
- offset double;
- % Number of discretization points for the time variable at which
- % the Gevrey-function is evaluated.
- N_t double;
- % Number of derivatives that need to be considered for this gevrey
- % function
- N_diff double;
- % Length of the time interval the Gevrey-function is defined on.
- T double = -1;
- % Gevrey-order
- order double = -1;
- end
-
- properties (Dependent)
- sigma;
- dt;
- end
-
- methods
- %% CONSTRUCTOR
- function obj = BasicVariable(valueContinuous, varargin)
-
- parentVarargin = {};
-
- if nargin > 0
-
- % make default grid:
- preParser = misc.Parser();
- preParser.addParameter('T', 1);
- preParser.addParameter('N_t', 101);
- preParser.addParameter('dt', []);
- preParser.addParameter('N_diff', 1);
- preParser.parse(varargin{:});
-
- if ~isempty(preParser.Results.dt)
- N_t = quantity.BasicVariable.setDt(preParser.Results.T, preParser.Results.dt);
- preResults.T = preParser.Results.T;
- preResults.dt = preParser.Results.dt;
- preResults.N_diff = preParser.Results.N_diff;
- else
- N_t = preParser.Results.N_t;
-
- preResults.T = preParser.Results.T;
- preResults.N_t = preParser.Results.N_t;
- preResults.N_diff = preParser.Results.N_diff;
- end
-
- grid = {linspace(0, preParser.Results.T, N_t)'};
- varargin = [varargin{:}, 'grid', {grid}];
- parentVarargin = {valueContinuous, varargin{:}, 'gridName', {'t'}};
- end
-
- obj@quantity.Function(parentVarargin{:});
-
- if nargin > 0
- % first parser
- p1 = misc.Parser();
- p1.addParameter('K', 1);
- p1.addParameter('diffShift', 0);
- p1.addParameter('offset', 0);
- p1.addParameter('order', 1.5);
- p1.parse(varargin{:});
- for parameter = fieldnames(p1.Results)'
- obj.(parameter{1}) = p1.Results.(parameter{1});
- end
- for parameter = fieldnames(preResults)'
- obj.(parameter{1}) = preParser.Results.(parameter{1});
- end
- end
- end
-
- function n = nargin(obj)
- n = 1;
- end
-
- function D = diff(obj, K)
- if nargin == 1
- i = 1;
- end
-
- for i = 1:numel(K)
- % create a new object for the derivative of obj:
- D(:,i) = obj.copy();
- [D.name] = deal(['d/dt (' num2str(K(i)) ') ']);
- [D.valueDiscrete] = deal([]);
- for l = 1:numel(obj)
- D(l,i).diffShift = D(l, 1).diffShift + K(i);
- end
- end
- end
-
- %% PROPERTIES
- function dt = get.dt(obj)
- dt = obj.T / obj.N_t;
- end
- function set.dt(obj, dt)
- obj.N_t = quantity.BasicVariable.setDt(obj.T, dt);
- end
-
- function obj = set.diffShift(obj, n)
- obj.diffShift = n;
- obj.valueDiscrete = [];
- end
-
- function s = get.sigma(obj)
- s = 1 / ( obj.order - 1);
- end
- function obj = set.sigma(obj, s)
- obj.order = 1 + 1/s;
- end
- function obj = set.order(obj, order)
- obj.order = obj.set_order(order);
- end
- function order = get.order(obj)
- order = obj.get_order(obj.order);
- end
- function obj = set.N_t(obj, value)
- obj.N_t = obj.set_N_t(value);
- end
- function N = get.N_t(obj)
- N = obj.get_N_t(obj.N_t);
- end
- function obj = set.T(obj, value)
- obj.T = obj.set_T(value);
- end
- function T = get.T(obj)
- T = obj.get_T(obj.T);
- end
-
- function b = copy_K(obj, K)
- b = obj.copy();
- b.valueDiscrete = K / obj.K * obj.valueDiscrete();
- b.K = K;
- end
-
- end
-
- methods (Access = protected)
-
- function f = getValueContinuous(obj, f)
- f = @(z) obj.K * f(z, obj.diffShift, obj.T, obj.sigma);
- end
-
- function order = set_order(obj, order)
- end
- function order = get_order(obj, order)
- end
- function T = get_T(obj, T)
- end
- function T = set_T(obj, T)
- end
- function n = get_n(obj, n)
- end
- function n = set_n(obj, n)
- end
- function N = get_N_t(obj, N)
- end
- function N = set_N_t(obj, N)
- end
- end
-
- methods (Access = protected, Static)
- function Nt = setDt(T, dt)
- t = 0:dt:T;
- Nt = length(t);
- if ~numeric.near(t(end), T)
- warning(['Sampling with dt=' num2str(dt) ' leads to t(end)= ' num2str(t(end))]);
- end
- end
- end
-
- methods (Static)
-
-
- function b = makeBasicVariable(b0, K)
-
- for k = 1:numel(K)
- for l = 1:size(b0, 2)
- b(k, l) = b0(l).copy_K(K(k));
- end
- end
-
- end
-
- end
+ %BasicVariable class for desription of basic variables
+ % The basic variable can be of the form
+ % g(t) = [phi_1(t); phi(2); ... ]
+ %
+ properties
+ % The constant gain of the Gevrey-function
+ K double = 1;
+ % Shifts the evaluation of the derivative of the Gevrey-function.
+ % Can be used if the series, computed using the gevrey function,
+ % should not start at the 0-th derivative.
+ diffShift double {mustBeNonnegative, mustBeInteger};
+ % Offset to raise the whole Gevrey-function g_ = g + offset.
+ offset double;
+ % Number of discretization points for the time variable at which
+ % the Gevrey-function is evaluated.
+ N_t double;
+ % Number of derivatives that need to be considered for this gevrey
+ % function
+ N_diff double;
+ % Length of the time interval the Gevrey-function is defined on.
+ T double = -1;
+ % Gevrey-order
+ order double = -1;
+ end
+
+ properties (Dependent)
+ sigma;
+ dt;
+ end
+
+ methods
+ %% CONSTRUCTOR
+ function obj = BasicVariable(valueContinuous, varargin)
+
+ parentVarargin = {};
+
+ if nargin > 0
+
+ % make default grid:
+ preParser = misc.Parser();
+ preParser.addParameter('T', 1);
+ preParser.addParameter('N_t', 101);
+ preParser.addParameter('dt', []);
+ preParser.addParameter('N_diff', 1);
+ preParser.parse(varargin{:});
+
+ if ~isempty(preParser.Results.dt)
+ N_t = quantity.BasicVariable.setDt(preParser.Results.T, preParser.Results.dt);
+ preResults.T = preParser.Results.T;
+ preResults.dt = preParser.Results.dt;
+ preResults.N_diff = preParser.Results.N_diff;
+ else
+ N_t = preParser.Results.N_t;
+
+ preResults.T = preParser.Results.T;
+ preResults.N_t = preParser.Results.N_t;
+ preResults.N_diff = preParser.Results.N_diff;
+ end
+
+ grid = {linspace(0, preParser.Results.T, N_t)'};
+ varargin = [varargin{:}, 'grid', {grid}];
+ parentVarargin = {valueContinuous, varargin{:}, 'gridName', {'t'}};
+ end
+
+ obj@quantity.Function(parentVarargin{:});
+
+ if nargin > 0
+ % first parser
+ p1 = misc.Parser();
+ p1.addParameter('K', 1);
+ p1.addParameter('diffShift', 0);
+ p1.addParameter('offset', 0);
+ p1.addParameter('order', 1.5);
+ p1.parse(varargin{:});
+ for parameter = fieldnames(p1.Results)'
+ obj.(parameter{1}) = p1.Results.(parameter{1});
+ end
+ for parameter = fieldnames(preResults)'
+ obj.(parameter{1}) = preParser.Results.(parameter{1});
+ end
+ end
+ end
+
+ function n = nargin(obj)
+ n = 1;
+ end
+
+ function D = diff(obj, K)
+ if nargin == 1
+ i = 1;
+ end
+
+ for i = 1:numel(K)
+ % create a new object for the derivative of obj:
+ D(:,i) = obj.copy();
+ [D.name] = deal(['d/dt (' num2str(K(i)) ') ']);
+ [D.valueDiscrete] = deal([]);
+ for l = 1:numel(obj)
+ D(l,i).diffShift = D(l, 1).diffShift + K(i);
+ end
+ end
+ end
+
+ %% PROPERTIES
+ function dt = get.dt(obj)
+ dt = obj.T / obj.N_t;
+ end
+ function set.dt(obj, dt)
+ obj.N_t = quantity.BasicVariable.setDt(obj.T, dt);
+ end
+
+ function obj = set.diffShift(obj, n)
+ obj.diffShift = n;
+ obj.valueDiscrete = [];
+ end
+
+ function s = get.sigma(obj)
+ s = 1 / ( obj.order - 1);
+ end
+ function obj = set.sigma(obj, s)
+ obj.order = 1 + 1/s;
+ end
+ function obj = set.order(obj, order)
+ obj.order = obj.set_order(order);
+ end
+ function order = get.order(obj)
+ order = obj.get_order(obj.order);
+ end
+ function obj = set.N_t(obj, value)
+ obj.N_t = obj.set_N_t(value);
+ end
+ function N = get.N_t(obj)
+ N = obj.get_N_t(obj.N_t);
+ end
+ function obj = set.T(obj, value)
+ obj.T = obj.set_T(value);
+ end
+ function T = get.T(obj)
+ T = obj.get_T(obj.T);
+ end
+
+ function b = copy_K(obj, K)
+ b = obj.copy();
+ b.valueDiscrete = K / obj.K * obj.valueDiscrete();
+ b.K = K;
+ end
+
+ end
+
+ methods (Access = protected)
+
+ function f = getValueContinuous(obj, f)
+ f = @(z) obj.K * f(z, obj.diffShift, obj.T, obj.sigma);
+ end
+
+ function order = set_order(obj, order)
+ end
+ function order = get_order(obj, order)
+ end
+ function T = get_T(obj, T)
+ end
+ function T = set_T(obj, T)
+ end
+ function n = get_n(obj, n)
+ end
+ function n = set_n(obj, n)
+ end
+ function N = get_N_t(obj, N)
+ end
+ function N = set_N_t(obj, N)
+ end
+ end
+
+ methods (Access = protected, Static)
+ function Nt = setDt(T, dt)
+ t = 0:dt:T;
+ Nt = length(t);
+ if ~numeric.near(t(end), T)
+ warning(['Sampling with dt=' num2str(dt) ' leads to t(end)= ' num2str(t(end))]);
+ end
+ end
+ end
+
+ methods (Static)
+
+
+ function b = makeBasicVariable(b0, K)
+
+ for k = 1:numel(K)
+ for l = 1:size(b0, 2)
+ b(k, l) = b0(l).copy_K(K(k));
+ end
+ end
+
+ end
+
+ end
end
\ No newline at end of file
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 9973f3d..9e12a59 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -59,7 +59,13 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
valueOriginalSize = size(valueOriginal);
S = num2cell(valueOriginalSize);
if any(valueOriginalSize == 0)
- obj = quantity.Discrete.empty(S{:});
+ % If the size is specified in the arguements, it should
+ % be chosen instead of the default size from the
+ % valueOriginal.
+ myParser = misc.Parser();
+ myParser.addParameter('size', [S{:}]);
+ myParser.parse(varargin{:});
+ obj = quantity.Discrete.empty(myParser.Results.size);
return;
end
@@ -151,7 +157,12 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
%---------------------------
% --- getter and setters ---
- %---------------------------
+ %---------------------------
+ function i = isConstant(obj)
+ % the quantity is interpreted as constant if it has no grid or
+ % it has a grid that is only one point.
+ i = isempty(obj.gridSize) || prod(obj.gridSize) == 1;
+ end
function doNotCopy = get.doNotCopy(obj)
doNotCopy = obj.doNotCopyPropertiesName();
end
@@ -178,7 +189,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
obj.gridName = name;
end
-
+
function set.grid(obj, grid)
if ~iscell(grid)
grid = {grid};
@@ -203,6 +214,21 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
function d = double(obj)
d = obj.on();
end
+
+ function o = quantity.Operator(obj)
+ for k = 1:size(obj, 3)
+ A{k} = obj(:,:,k);
+ end
+ o = quantity.Operator(A, 'grid', obj(1).grid);
+ end
+
+ function o = quantity.Symbolic(obj)
+ if isempty(obj)
+ o = quantity.Symbolic.empty(size(obj));
+ else
+ error('Not yet implemented')
+ end
+ end
end
methods (Access = public)
@@ -241,6 +267,9 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
tempInterpolant = numeric.interpolant(...
[obj(1).grid, indexGrid{:}], value);
value = tempInterpolant.evaluate(myGrid{:}, indexGrid{:});
+ elseif obj.isConstant
+ value = repmat(value, [cellfun(@(v) numel(v), myGrid), ones(1, length(size(obj)))]);
+ gridPermuteIdx = 1:numel(myGrid);
end
value = permute(reshape(value, [cellfun(@(v) numel(v), myGrid), size(obj)]), ...
[gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
@@ -252,8 +281,13 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% further quantites are inputs via varargin, it is verified if
% that quantity has same grid and gridName as quantity a
% as well.
- referenceGridName = a(1).gridName;
- referenceGrid= a(1).grid;
+ if isempty(a)
+ referenceGridName = '';
+ referenceGrid = {};
+ else
+ referenceGridName = a(1).gridName;
+ referenceGrid= a(1).grid;
+ end
for it = 1 : numel(a)
assert(isequal(referenceGridName, a(it).gridName), ...
'All elements of a quantity must have same gridNames');
@@ -279,10 +313,21 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
function [referenceGrid, referenceGridName] = getFinestGrid(a, varargin)
% find the finest grid of all input quantities by comparing
% gridSize for each iteratively.
- referenceGridName = a(1).gridName;
- referenceGrid = a(1).grid;
- referenceGridSize = a(1).gridSize(referenceGridName);
+
+ if isempty(a)
+ referenceGridName = '';
+ referenceGrid = {};
+ referenceGridSize = [];
+ else
+ referenceGridName = a(1).gridName;
+ referenceGrid = a(1).grid;
+ referenceGridSize = a(1).gridSize(referenceGridName);
+ end
+
for it = 1 : numel(varargin)
+ if isempty(varargin{it})
+ continue;
+ end
assert(numel(referenceGridName) == numel(varargin{it}(1).gridName), ...
['For getFinestGrid, the gridName of all objects must be equal', ...
'. Maybe gridJoin() does what you want?']);
@@ -416,7 +461,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% Copyright 1984-2005 The MathWorks, Inc.
% Built-in function.
if nargin == 1
- objCell = {a};
+ objCell = {a};
else
objCell = [{a}, varargin(:)'];
@@ -599,6 +644,9 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if nargin == 1 || isempty(gridName2Replace)
% if gridName2Replace is empty, then nothing must be done.
solution = obj;
+ elseif isempty(obj);
+ % if the object is empty, nothing must be done.
+ solution = obj;
else
% input checks
assert(nargin == 3, ['Wrong number of input arguments. ', ...
@@ -776,21 +824,23 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
end
- function h = plot(obj, varargin)
-
+ function H = plot(obj, varargin)
+ H = [];
p = misc.Parser();
p.addParameter('fig', []);
- p.addParameter('dock', false);
+ p.addParameter('dock', quantity.Settings.instance().dockThePlot);
p.parse(varargin{:});
- %misc.struct2ws(p.Results);
+ additionalPlotOptions = misc.struct2namevaluepair(p.Unmatched);
+
fig = p.Results.fig;
dock = p.Results.dock;
for figureIdx = 1:size(obj, 3)
if isempty(p.Results.fig)
h = figure();
else
- h = figure(fig + figureIdx - 1);
+ h = figure(fig + figureIdx - 1);
end
+ H = [H, h];
if dock
set(h, 'WindowStyle', 'Docked');
@@ -805,6 +855,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
i = 1: numel(obj(:,:,figureIdx));
i = reshape(i, size(obj, 2), size(obj, 1))';
+ if obj.gridSize == 1
+ additionalPlotOptions = [additionalPlotOptions(:)', ...
+ {'x'}];
+ end
+
for rowIdx = subplotRowIdx
for columnIdx = subpotColumnIdx
subplot(size(obj, 1), size(obj, 2), i(rowIdx, columnIdx));
@@ -813,11 +868,13 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
plot(...
obj(rowIdx, columnIdx, figureIdx).grid{1}(:), ...
- obj(rowIdx, columnIdx, figureIdx).valueDiscrete );
+ obj(rowIdx, columnIdx, figureIdx).valueDiscrete, ...
+ additionalPlotOptions{:});
elseif obj.nargin() == 2
misc.isurf(obj(rowIdx, columnIdx, figureIdx).grid{1}(:), ...
obj(rowIdx, columnIdx, figureIdx).grid{2}(:), ...
- obj(rowIdx, columnIdx, figureIdx).valueDiscrete);
+ obj(rowIdx, columnIdx, figureIdx).valueDiscrete, ...
+ additionalPlotOptions{:});
ylabel(labelHelper(2), 'Interpreter','latex');
else
error('number inputs not supported');
@@ -854,7 +911,6 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
myText = strrep(myText, 'gamma', '\gamma');
myText = strrep(myText, 'Delta', '\Delta');
myText = strrep(myText, 'delta', '\delta');
- myText = strrep(myText, 'int', '\int');
if ~contains(myText, '\zeta') && ~contains(myText, '\Zeta')
myText = strrep(myText, 'eta', '\eta');
end
@@ -900,6 +956,9 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% change the grid of the obj quantity. The order of grid and
% gridName in the obj properties remains unchanged, only the
% data points are exchanged.
+ if isempty(obj)
+ return;
+ end
gridIndexNew = obj(1).gridIndex(gridNameNew);
myGrid = cell(1, numel(obj(1).grid));
myGridName = cell(1, numel(obj(1).grid));
@@ -966,6 +1025,18 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
ftr = false; % flag for the special case a == const.
+ if isempty(b) || isempty(a)
+ % TODO: actually this only holds for multiplication of
+ % matrices. If higher dimensional arrays are multiplied
+ % this can lead to wrong results.
+ sa = size(a);
+ sb = size(b);
+
+ P = quantity.Discrete.empty([sa(1:end-1) sb(2:end)]);
+ return
+ end
+
+
if isa(a, 'double')
if numel(a) == 1
% simple multiplication in scalar case
@@ -979,7 +1050,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
if a.nargin() == 0
- % TODO: Whats the use of this?
+ % If the first argument a is constant value, then bad
+ % things will happen. To avoid this, we calculate
+ % a * b = (b' * a')'
+ % instead. Thus we have to exchange both values and set the
+ % flag, that we have to transpose the result in the end.
a_ = a;
a = b';
b = a_';
@@ -1362,6 +1437,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
%% PLUS is the sum of two quantities.
assert(isequal(size(A), size(B)), 'plus() not supports mismatching sizes')
+ if isempty(A) || isempty(B)
+ C = quantity.Discrete.empty(size(A));
+ return
+ end
+
% combine both domains with finest grid
[gridJoined, gridNameJoined] = gridJoin(A, B);
gridJoinedLength = cellfun(@(v) numel(v), gridJoined);
@@ -1413,7 +1493,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
P = A.copy();
- if ~isa(B, 'quantity.quantiy')
+ if ~isa(B, 'quantity.Discrete')
B = quantity.Discrete(B);
end
supremum = nan(size(A));
@@ -1431,6 +1511,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
maxValue = reshape(max(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
end
+ function maxValue = MAX(obj)
+ maxValue = obj.max();
+ maxValue = max(maxValue(:));
+ end
+
function minValue = min(obj)
% max returns the minimum value of all elements of A over all
% variables as a double array.
@@ -1465,6 +1550,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% variable, this one is chosen. The 'fixedLowerBound' option
% defines if a = 0 or a = s.
+ if isempty(x)
+ f = quantity.Discrete.empty(size(K, 1), 0);
+ return
+ end
+
ip = misc.Parser();
ip.addParameter('variable', x(1).gridName{1});
ip.addParameter('fixedLowerBound', true);
@@ -1673,7 +1763,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
function [idx, permuteGrid] = computePermutationVectors(a, b)
% Computes the required permutation vectors to use
% misc.multArray for multiplication of matrices
-
+
uA = unique(a(1).gridName);
assert(numel(uA) == numel(a(1).gridName), 'Gridnames have to be unique!');
uB = unique(b(1).gridName);
diff --git a/+quantity/Operator.m b/+quantity/Operator.m
index ceae683..976dafa 100644
--- a/+quantity/Operator.m
+++ b/+quantity/Operator.m
@@ -1,127 +1,305 @@
-classdef Operator < quantity.Symbolic
- %OPERATOR system operator of the form A[x(t)](z)
- % This class describes system opertors of a system of first order
- % PDE wrt. space of the form
- % dz x(z,t) = A[x(t)](z)
- % where A can be written as
- % A[x(t)](z) = A0(z) x(z,t) + A1(z) dz x(z,t) + A2(z) dz^2 x(z,t)
- % ... + Ak dz^k x(z,t)
-
- properties
- operatorVar;
- end
- properties (SetAccess = protected)
+classdef Operator < handle
+%OPERATOR define operator matrices
+% This class describes operators of the form
+% A[x(t)](z) = A0(z) x(z,t) + A1(z) dz x(z,t) + A2(z) dz^2 x(z,t)
+% ... + Ak dz^k x(z,t)
+
+ properties
+ % The variable that is used for the algebraic representation of the
+ % operator
+ s;
+ % continuous variable
+ variable;
+ % The values of the operator matrices in form of a
+ % quantity.Discrete object
+ coefficient quantity.Discrete;
+ % The grid, on which the operator is defined
+ grid;
+ end
+ properties (Hidden, SetAccess = protected)
% Fundamental matrix of this operator, that is the solution of
% dz phi(z,z0) = A(z,s) phi(z,z0)
- % phi(z0, z0) = I
- % A(z,s) = sum_k Ak(z) * s^k
+ % phi(z0, z0) = I A(z,s) = sum_k Ak(z) * s^k
Phi0;
end
-
- methods
- function obj = Operator(A, varargin)
-
- parentVarargin = {};
- if nargin > 0
- parentVarargin = {A, varargin{:}};
- end
-
- %OPERATOR
- obj@quantity.Symbolic(parentVarargin{:});
-
- if nargin > 0
- p = misc.Parser();
- p.addParameter('operatorVar', sym('s'));
- p.parse(varargin{:});
- misc.struct2ws(p.Results);
+
+ methods
+ function obj = Operator(A, varargin)
+
+ if nargin > 0
+
+ gridParser = misc.Parser();
+ gridParser.addOptional('grid', linspace(0,1)');
+ gridParser.addOptional('variable', sym('z', 'real'));
+ gridParser.addOptional('s', sym('s'));
+ gridParser.parse(varargin{:});
+
+ grid = gridParser.Results.grid;
+ variable = gridParser.Results.variable;
+ s = gridParser.Results.s;
+
+ %TODO assert that A has entries with the same size
+ for k = 1 : numel(A)
+ obj(k).coefficient = A{k}; %#ok<AGROW>
+ end
+
+ [obj.grid] = deal(grid);
+ [obj.variable] = deal(variable);
+ [obj.s] = deal(s);
+ end
+ end
+
+ function set.grid(obj, grid)
+ if ~iscell(grid)
+ grid = {grid};
+ end
+
+ isV = cellfun(@(v) (sum(size(v)>1) == 1) || (numel(v) == 1), grid); % also allow 1x1x1x41 grids
+ assert(all(isV(:)), 'Please use vectors for the grid entries!');
+
+ [obj.grid] = deal(grid);
+ end
+
+ end
+
+ methods (Access = public)
+
+ function c = copy(obj, W)
+ c = quantity.Operator(W, 'grid', obj(1).grid, 'variable', obj(1).variable, 's', obj(1).s);
+ end
+
+ function m = M(obj) % TODO: This function should be protected actually.
+ m = builtin('cat', 3, obj.coefficient);
+ end
+
+ function h = plot(obj, varargin)
+ h = obj.M.plot(varargin{:});
+ end
+
+ function S = powerSeries(obj)
+ S = misc.powerSeries2Sum(obj.M.sym, obj(1).s);
+ end
+
+ function A = ctranspose(obj)
+ for k = 1:length(obj)
+ M{k} = (-1)^(k-1) * obj(k).coefficient'; %#ok<AGROW>
+ end
+ A = obj.copy(M);
+ end
+
+ function [Phi, Psi] = stateTransitionByOde(A, varargin)
+ % STATETRANSTITIONBYODE
+ % [Phi, Psi] = stateTransitionByOde(A, varargin) computes the
+ % state-transition matrix of the boundary-value-problem
+ % dz x(z,s) = A(z,s) x(xi,s) + B(z,s) * u(s)
+ % so that Phi is the solution of
+ % dz Phi(z,xi,s) = A(z,s) Phi(z,xi,s)
+ % and Psi is given by
+ % Psi(z,xi,s) = int_xi_z Phi(z,zeta,s) B(zeta,s) d zeta
+ % if B(z,s) is defined as quantity.Operator object by
+ % name-value-pair.
+ % To solve this problem the recursion formula from [1]
+ % dz w_k(z) = sum_j=0 ^d A_j(z) * w_{k-j}(z) + B_k nu
+ % is used. This recursion of odes can be reformulated to a
+ % large ode:
+ %
+ % dz w_0 = A_0 w_0 + B_0 nu
+ % dz w_1 = A_1 w_0 + A_0 w_0 + B_1 nu
+ % dz w_2 = A_2 w_0 + A_1 w_1 + A_0 w_2 + B_2 nu
+ % ...
+ %
+ % [dz w_0] [ A_0 0 0 0 ... ][w_0] [B_0]
+ % [dz w_1] [ A_1 A_0 0 0 ... ][w_1] [B_1]
+ % [dz w_2] = [ A_2 A_1 A_0 0 ... ][w_2] + [B_2] nu
+ % [ ... ] [ ... [[...] [...]
+ %
+ % By some sorting of the terms, as shown in [2], the solution
+ % of the recursion can be described by
+ % w_k(z) = Phi_k(z, xi) gamma + Psi_k(z,xi) nu
+ % with w_0(xi) = gamma and w_k(xi) = 0, k > 0. From the
+ % solution for the ode of the recursion:
+ % W(z) = PHI(z, xi) W(0) + PSI(z,xi) nu
+ % and W(z) = [w_0'(z), w_1'(z), w_2'(z), ... ]' it is obvious,
+ % that the required Phi_k(z,xi) and Psi_k(z,xi) can be
+ % determined by appropriate partitioning of PHI(z,xi),
+ % PSI(z,xi):
+ %
+ % [w_0] [ Phi_0 * * ... ][gamma] [Psi_0]
+ % [w_1] [ Phi_1 * * ... ][ 0 ] [Psi_1]
+ % [w_2] = [ Phi_2 * * ... ][ 0 ] + [Psi_2] nu
+ % [...] [ ... ][ ... ] [ ... ]
+ %
+
+ n = size(A(1).coefficient, 2);
+
+ ip = misc.Parser();
+ ip.addOptional('N', 1);
+ ip.addOptional('B', quantity.Operator.empty(1, 0));
+ ip.parse(varargin{:});
+ B = ip.Results.B;
+ N = ip.Results.N;
+
+ if isempty(B)
+ m = 0;
+ else
+ m = size(B(1).coefficient, 2);
+ end
+
+ n = size(A(1).coefficient, 1);
+ d = length(A);
+ spatialDomain = A(1).coefficient(1).grid{:};
+
+ % build the large system of odes from the recursion
+ A_full = zeros(length(spatialDomain), d*n, d*n);
+ B_full = zeros(length(spatialDomain), d*n, m);
+
+ idx = 1:n;
+ for k = 1:N
+ idxk = idx + (k-1)*n;
+ if length(B) >= k
+ B_full(:, idxk, :) = B(k).coefficient.on();
+ end
+ for j = 1:d
+ idxj = idx + (j-1)*n;
+ if k - (j-1) > 0
+ A_full(:, idxk, idxj) = A(k - (j-1)).coefficient.on();
+ end
+ end
+ end
+ f = quantity.Discrete(...
+ misc.fundamentalMatrix.odeSolver_par( A_full, spatialDomain ), ...
+ 'grid', {spatialDomain, spatialDomain}, 'gridName', {'z', 'zeta'});
+
+ p = f.volterra(quantity.Discrete(B_full, 'grid', spatialDomain, ...
+ 'gridName', 'zeta'));
+
+ for k = 1:N
+ idxk = idx + (k-1)*n;
+ Phi(:,:,k) = f(idxk, idx).subs('zeta', 0);
+ Psi(:,:,k) = p(idxk,:);
+ end
+
+ end
+
+ function [Phi, Psi] = stateTransitionMatrix(A, varargin)
+ % STATETRANSITIONMATRIX computation of the
+ % state-transition-matrix
+ % [Phi, Psi] = state-transition-matrix(obj, varargin) computes the
+ % state-transition-matrix of the boundary value problem
+ % dz x(z,s) = A(z,s) x(xi,s) + B(z,s) * u(s)
+ % so that Phi is the solution of
+ % dz Phi(z,xi,s) = A(z,s) Phi(z,xi,s)
+ % and Psi is given by
+ % Psi(z,xi,s) = int_xi_z Phi(z,zeta,s) B(zeta,s) d zeta
+ % if B(z,s) is defined as quantity.Operator object by
+ % name-value-pair.
+
+ n = size(A(1).coefficient, 2);
+
+ ip = misc.Parser();
+ ip.addOptional('N', 1);
+ ip.addOptional('B', quantity.Operator.empty(1, 0));
+ ip.parse(varargin{:});
+ B = ip.Results.B;
+ N = ip.Results.N;
+
+ if isempty(B)
+ m = 0;
+ else
+ m = size(B(1).coefficient, 2);
+ end
+ d = length(A); % As order of a polynom is zero based
+ z = A(1).grid{1};
- for par = fieldnames(p.Results)'
- [obj.(par{1})] = deal(p.Results.(par{1}));
- end
- end
- end
- end
-
- methods (Access = public)
- function S = powerSeries(obj)
- S = misc.powerSeries2Sum(obj.valueSymbolic, obj.operatorVar);
- end
-
-
- function A = ctranspose(obj)
-
- argin = struct(obj(1));
- argin.name = [argin.name, ''''];
- argin = misc.struct2namevaluepair(argin);
-
- a = obj.sym;
- for k = 1:size(obj,3)
- a(:,:,k) = (-1)^(k-1) * a(:,:,k)';
- end
- A = quantity.Operator(a, argin{:});
- end
-
- function [Phi, F, Psi, P] = fundamentalMatrix(obj, varargin)
-
- ip = misc.Parser();
- ip.addOptional('N', 5);
- ip.addOptional('B', []);
- ip.parse(varargin{:});
- misc.struct2ws(ip.Results)
-
- [Phi, F, Psi, P] = ...
- misc.fundamentalMatrix.rudolphWoittennek(obj.sym, ...
- 'N', N, 'zeta', obj(1).grid{1}(1), 's', obj(1).operatorVar, ...
- 'B', B.sym);
-
- if isempty(B)
- Psi = zeros(size(B));
- P = zeros([size(B), N]);
- end
-
- end
-
- function [Phi, Psi] = fundamentalMatrixNumeric(obj, varargin)
- ip = misc.Parser();
- ip.addOptional('N', 5);
- ip.addOptional('B', quantity.Discrete.empty([size(obj,1), 0]));
- ip.addOptional('spatialDomain', obj(1).grid{:});
- ip.parse(varargin{:});
- misc.struct2ws(ip.Results)
-
- [Phi, Psi] = ...
- misc.fundamentalMatrix.rudolphWoittennekNumeric(...
- obj, N, B, spatialDomain);
+ % initialization of the ProgressBar
+ counter = N - 1;
+ pbar = misc.ProgressBar(...
+ 'name', 'Trajectory Planning (Woittennek)', ...
+ 'terminalValue', counter, ...
+ 'steps', counter, ...
+ 'initialStart', true);
+
+ %% computation of transition matrix as power series
+ % The fundamental matrix is considered in the laplace space as the state
+ % transition matrix of the ODE with complex variable s:
+ % dz w(z,s) = A(z, s) w(z,s) * B(z) v(s),
+ % Hence, fundamental matrix is a solution of the ODE:
+ % dz Phi(z,zeta) = A(Z) Phi(z,zeta)
+ % Phi(zeta,zeta) = I
+ %
+ % Psi(z,xi,s) = int_zeta^z Phi(z, xi, s) B(xi) d_xi
+ %
+ % At this, A(z, s) can be described by
+ % A(z,s) = A_0(z) + A_1(z) s + A_2(z) s^2 + ... + A_d(z) s^d
+ %
+ % With this, the recursion formula
+ % Phi_k(z, zeta) = int_zeta^z Phi_0(z, xi) * sum_i^d A_i(xi) *
+ % Phi_{k-i}(z, xi) d_xi
+ % Psi_k(z, zeta) = int_zeta^z Phi_0(z, xi) * sum_i^d A_i(xi) *
+ % Psi_{k-i}(z, xi) + B_k(xi) d_xi
+
+
+ F0 = A.computePhi0();
+ Phi(:,:,1) = F0.subs('zeta', 0);
+ Psi(:,:,1) = F0.volterra(B(1).coefficient.subs('z', 'zeta'));
+
+ On = quantity.Discrete.zeros([n n], z, 'gridName', 'z');
+ Om = quantity.Discrete.zeros([n m], z, 'gridName', 'z');
+ for k = 2 : N
+ pbar.raise();
+ dA = min(d, k);
+ % compute the temporal matrices:
+ M = On;
+ N = Om;
+ for l = 2 : dA
+ M = M + A(l).coefficient * Phi(:, :, k-l+1);
+ N = N + A(l).coefficient * Psi(:, :, k-l+1);
+ end
+ if size(B, 3) >= k
+ N = N + B(k).coefficient(:, :);
+ end
+
+ % compute the integration of the fundamental matrix with the temporal
+ % values:
+ % Phi_k(z) = int_0_z Phi_0(z, xi) * M(xi) d_xi
+ Phi(:, :, k) = F0.volterra( M.subs('z', 'zeta') );
+ Psi(:, :, k) = F0.volterra( N.subs('z', 'zeta') );
+ end
+ %
+ pbar.stop();
+
+
end
- function [F0] = fundamentalMatrix0(obj, spatialDomain)
- % Fundamental matrix of this operator, that is the solution of
- % dz phi(z,z0) = A(z,s) phi(z,z0)
- % phi(z0, z0) = I
- % A(z,s) = sum_k Ak(z) * s^k
+ function [F0] = computePhi0(obj, spatialDomain)
+ % COMPUTEPHI0 compute the static state-transition-matrix
+ % [F0] = computePhi0(obj, spatialDomain) computation of the
+ % state-transition matrix for this operator for the
+ % steady-state case, i.e., s=0.
+ % dz phi(z, z0) = A(z, s=0) phi(z,z0)
+ % phi(z0, z0) = I
% todo@ferdinand: Es sollte die möglichkeit geben den Operator
% mit neuen Eigenschaften (wie z.B. einem neuen grid) zu
% aktualisieren. Momentan wird er nur einmal berechnet.
- if isempty(obj(1).Phi0)
- if nargin == 1
- assert(obj(1).nargin == 1);
- spatialDomain = obj(1).grid{1};
- Adisc = obj(:,:,1).on();
- else
- Adisc = obj(:,:,1).on(spatialDomain);
- end
+ if isempty(obj(1).Phi0)
+ if nargin == 1
+ assert(numel(obj(1).variable) == 1);
+ spatialDomain = obj(1).grid{1};
+ end
- if all([obj.isConstant])
+ if all(cellfun(@(e) e.isConstant, {obj.coefficient}))
% for a constant system operator we can use the matrix
% exponential function:
z = sym('z', 'real');
zeta = sym('zeta', 'real');
- f0 = expm(obj(:,:,1).on()*(z - zeta));
+ f0 = expm(squeeze(obj(1).coefficient.on())*(z - zeta));
F0 = quantity.Symbolic(f0, 'grid', {spatialDomain, spatialDomain});
else
- f0 = misc.fundamentalMatrix.odeSolver_par( Adisc, spatialDomain );
+ f0 = misc.fundamentalMatrix.odeSolver_par( ...
+ obj(1).coefficient.on(spatialDomain), ...
+ spatialDomain );
F0 = quantity.Discrete(f0, ...
'size', [size(obj, 1), size(obj, 2)], ...
'gridName', {'z', 'zeta'}, ...
@@ -132,17 +310,17 @@ classdef Operator < quantity.Symbolic
F0 = obj(1).Phi0;
- end
-
- end
-
- methods (Static, Access = protected)
- function var = getVariable(symbolicFunction, obj)
- var = getVariable@quantity.Symbolic(symbolicFunction);
- idxS = arrayfun(@(s) isequal(s, obj.operatorVar), var);
- var = var(~idxS);
- end
- end
-
+ end
+
+ end
+
+ methods (Static, Access = protected)
+ function var = getVariable(symbolicFunction, obj)
+ var = getVariable@quantity.Symbolic(symbolicFunction);
+ idxS = arrayfun(@(s) isequal(s, obj.operatorVar), var);
+ var = var(~idxS);
+ end
+ end
+
end
diff --git a/+quantity/Settings.m b/+quantity/Settings.m
new file mode 100644
index 0000000..1093580
--- /dev/null
+++ b/+quantity/Settings.m
@@ -0,0 +1,27 @@
+classdef Settings < handle
+ %SETTINGS object to save some settings for the quantity class
+ %
+
+ properties (Access = public)
+ dockThePlot = false
+ end
+
+ methods (Access = private)
+ function obj = Settings()
+ end
+ end
+
+ methods (Static)
+ function obj = instance()
+ persistent uniqueObj
+ if isempty(uniqueObj)
+ obj = quantity.Settings();
+ uniqueObj = obj;
+ else
+ obj = uniqueObj;
+ end
+ end
+ end
+
+end
+
diff --git a/+quantity/Symbolic.m b/+quantity/Symbolic.m
index ace9702..33d381a 100644
--- a/+quantity/Symbolic.m
+++ b/+quantity/Symbolic.m
@@ -4,11 +4,6 @@ classdef Symbolic < quantity.Function
valueSymbolic sym;
variable sym;
end
-
- properties (Dependent)
- isConstant logical;
- end
-
properties (Constant)
defaultSymVar = sym('z', 'real');
end
@@ -98,8 +93,14 @@ classdef Symbolic < quantity.Function
end
end
- function i = get.isConstant(obj)
- i = isempty(symvar(obj.valueSymbolic));
+ function i = isConstant(obj)
+ i = true;
+ for k = 1:numel(obj)
+ i = i && isempty(symvar(obj(k).sym));
+ if ~i
+ break
+ end
+ end
end
function i = eq(A, B)
@@ -128,11 +129,6 @@ classdef Symbolic < quantity.Function
end
- function operator = quantity.Operator(obj)
- argin = getCopyParameters(obj(1));
- operator = quantity.Operator(obj.sym, argin{:});
- end
-
function Q = quantity.Discrete(obj, varargin)
% Cast of a quantity.Symbolic object into a quantity.Discrete
% object.
@@ -246,14 +242,22 @@ classdef Symbolic < quantity.Function
objIdx = find(isAquantityDiscrete, 1);
obj = objCell{objIdx};
+ if ~isempty(obj)
+ myVariable = obj(1).variable;
+ myGrid = obj(1).grid;
+ else
+ myVariable = '';
+ myGrid = {};
+ end
+
for k = 1:numel(objCell)
if isa(objCell{k}, 'quantity.Symbolic')
o = objCell{k};
elseif isa(objCell{k}, 'double') || isa(objCell{k}, 'sym')
o = quantity.Symbolic( objCell{k}, ...
- 'variable', obj(1).variable, ...
- 'grid', obj(1).grid);
+ 'variable', myVariable, ...
+ 'grid', myGrid);
end
objCell{k} = o;
end
@@ -493,9 +497,7 @@ classdef Symbolic < quantity.Function
if nargin == 1
k = 1;
end
- % TODO: if result is constant, delete gridName
% TODO: specify for which variable it should be differentiated
- % create a new object for the derivative of obj:
D = obj.copy();
[D.name] = deal(['(d_{' char(obj(1).variable) '} ' D(1).name, ')']);
[D.valueDiscrete] = deal([]);
@@ -509,12 +511,19 @@ classdef Symbolic < quantity.Function
mObj = uminus@quantity.Function(obj);
for k = 1:numel(obj)
mObj(k).valueSymbolic = - obj(k).valueSymbolic;
+ mObj(k).name = ['-', obj(k).name];
end
- [mObj.name] = deal(['-', obj(1).name]);
end
function C = mtimes(A, B)
% if one input is ordinary matrix, this is very simple.
+ if isempty(B) || isempty(A)
+ % TODO: actually this only holds for multiplication of
+ % matrices. If higher dimensional arrays are multiplied
+ % this can lead to wrong results.
+ C = quantity.Symbolic( mtimes@quantity.Discrete(A, B) );
+ return
+ end
if isnumeric(B)
C = quantity.Symbolic(A.sym() * B, ...
'grid', A(1).grid, 'variable', A(1).variable, ...
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index cd0a07a..0a75b91 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -684,6 +684,17 @@ testCase.verifyTrue(numeric.near(bc, BC.on));
end
+function testIsConstant(testCase)
+
+C = quantity.Discrete(rand(3,7), 'gridName', {});
+testCase.verifyTrue(C.isConstant());
+
+z = linspace(0, pi)';
+A = quantity.Discrete(sin(z), 'grid', z, 'gridName', 'z');
+testCase.verifyFalse(A.isConstant());
+
+end
+
function testMTimesConstant(testCase)
myGrid = linspace(0, 2*pi)';
diff --git a/+unittests/+quantity/testOperator.m b/+unittests/+quantity/testOperator.m
index 46b9acb..7ba6994 100644
--- a/+unittests/+quantity/testOperator.m
+++ b/+unittests/+quantity/testOperator.m
@@ -2,149 +2,168 @@ function [tests] = testOperator()
%TESTGEVREY Summary of this function goes here
% Detailed explanation goes here
tests = functiontests(localfunctions);
-
-% FIXME write useful unittests
-
end
-function testCastToQuantity(testCase)
+function setupOnce(testCase)
%%
-N = 8;
z = sym('z', 'real');
-Z = linspace(0, 1, 500)';
-A = quantity.Operator(cat(3, ...
- [ 0, 0, 0, 0; ...
- 1 0 0 0; ...
- 0 1 0 0; ...
- 0 0 1 0], ...
- [0 0 0, -3 + z; ...
- 0 0 0 0; ...
- 0 0 0 0; ...
- 0 0 0 0], ...
- [0 0 0, -4; ...
- 0 0 0 0; ...
- 0 0 0 0; ...
- 0 0 0 0]...
- ), 'name', 'A', 'grid', Z);
-
-%%
-X1 = quantity.Discrete(A);
-X2 = quantity.Discrete(A);
-
-verifyEqual(testCase, X1.on, A.on);
-verifyEqual(testCase, X2, X1);
+s = sym('s');
+Z = linspace(0, 1, 11)';
+
+% init with symbolic values
+A = cat(3, [ 0, 0, 0, 0; ...
+ 1, 0, 0, 0; ...
+ 0, 1, 0, 0; ...
+ 0, 0, 1, 0], ...
+ [0, 0, 0, -3 + z; ...
+ 0, 0, 0, 0; ...
+ 0, 0, 0, 0; ...
+ 0, 0, 0, 0], ...
+ [0, 0, 0, -4; ...
+ 0, 0, 0, 0; ...
+ 0, 0, 0, 0; ...
+ 0, 0, 0, 0]);
+
+B = cat(3, [1, 0, 0, 0]', zeros(4, 1), zeros(4, 1));
+
+for k = 1:3
+ a{k} = quantity.Symbolic(A(:,:,k), 'grid', Z, 'variable', z);
end
+
+A = quantity.Operator(a, 'grid', Z, 'variable', z, 's', s);
-function testPhiPsiSpaceDependent(testCase)
-%%
-a = 20;
-N = 1;
-z = sym('z', 'real');
-Z = linspace(0, 1, 12)';
-A = quantity.Operator([1, z; 0, a], 'name', 'A', 'grid', Z);
-B = quantity.Operator([0; 0], 'name', 'B', 'grid', Z, 'variable', z);
-
-[PhiNum, PsiNum] = A.fundamentalMatrixNumeric(N, B);
-[~, F, ~, P] = A.fundamentalMatrix(N, B);
-
-PhiSym = quantity.Operator(F, 'grid', Z, 'figureID', N + 1);
-PsiSym = quantity.Operator(P, 'grid', Z, 'variable', z);
-
-% plot(PhiSym, 'fig', 4)
-
-% use the exact solution from the peano-baker series paper.
-if numeric.near(a, 1)
- f = 0.5 * z^2 * exp(z);
-else
- f = (exp(z) - exp(a * z) - (1 - a) * z * exp( a * z ) ) / (1 - a)^2;
-end
-
-PhiExact = quantity.Symbolic([exp(z), f; 0, exp(a*z)], 'grid', Z, 'name', 'Phi_A');
-
-Dnum = relativeErrorSupremum(PhiNum, PhiExact);
-Dsym = relativeErrorSupremum(PhiSym, PhiExact);
-deltaNum = Dnum.on();
-deltaSym = Dsym.on();
-testCase.verifyLessThan(max(deltaNum(:)), 1e-6);
-testCase.verifyLessThan(max(deltaSym(:)), 1e-6);
+ testCase.TestData.A = A;
end
-function testPhiPsi(testCase)
-
-%%
-N = 3;
-z = sym('z', 'real');
-Z = linspace(0, 1, 15)';
-A = quantity.Operator(cat(3, ...
- [ 0, 0, 0, 0; ...
- 1 0 0 0; ...
- 0 1 0 0; ...
- 0 0 1 0], ...
- [0 0 0, -3 + z; ...
- 0 0 0 0; ...
- 0 0 0 0; ...
- 0 0 0 0], ...
- [0 0 0, -4; ...
- 0 0 0 0; ...
- 0 0 0 0; ...
- 0 0 0 0]...
- ), 'name', 'A', 'grid', Z);
-
-B = quantity.Operator([1 1; zeros(3,2)], 'name', 'B', 'variable', z, 'grid', Z);
-
-[PhiNum, PsiNum] = A.fundamentalMatrixNumeric(N, B, 'spatialDomain', Z);
-
-[PhiNum.figureID] = deal(6);
-[~, F, ~, P] = A.fundamentalMatrix(N, B);
-
-PhiSym = quantity.Operator(F, 'grid', Z);
-PsiSym = quantity.Operator(P, 'grid', Z);
-
-[D, sup] = relativeErrorSupremum(PhiNum, PhiSym);
-delta = D.on();
-testCase.verifyLessThan(max(delta(:)), 1e-2);
-
-D = relativeErrorSupremum(PsiNum, PsiSym);
-delta = D.on();
-testCase.verifyLessThan(max(delta(:)), 1e-2);
+function testCopy(testCase)
+ B = {quantity.Discrete(0, 'gridName', ''), ...
+ quantity.Discrete(1, 'gridName', ''), ...
+ quantity.Discrete(2, 'gridName', '')};
+ A = testCase.TestData.A;
+ C = A.copy(B);
+
+ testCase.verifyEqual(C(1).grid, A(1).grid);
+ testCase.verifyEqual(C(1).variable, A(1).variable);
+end
+function testPlot(testCase)
+ f = testCase.TestData.A.plot();
+ close(f);
+end
+function testPowerSeries(testCase)
+ syms z s
+ P = testCase.TestData.A.powerSeries();
+ testCase.verifyEqual(P(1, 4), - 4*s^2 + (z - 3)*s);
+end
+function testCTranspose(testCase)
+ At = testCase.TestData.A';
+ testCase.verifyTrue(At(3).coefficient == testCase.TestData.A(3).coefficient');
end
+%
+function testFundamentalMatrixSpaceDependent(testCase)
+%
+ a = 20;
+ z = sym('z', 'real');
+ Z = linspace(0, 1, 32)';
+ A = quantity.Operator(...
+ {quantity.Symbolic([1, z; 0, a], 'variable', 'z', 'grid', Z)}, ...
+ 'grid', Z);
+
+ F = A.computePhi0();
+ F0 = quantity.Discrete(squeeze(F.on({Z, 0})), 'gridName', 'z');
+
+ % compute the exact solution from the peano-baker series paper.
+ if numeric.near(a, 1)
+ f = 0.5 * z^2 * exp(z);
+ else
+ f = (exp(z) - exp(a * z) - (1 - a) * z * exp( a * z ) ) / (1 - a)^2;
+ end
+
+ PhiExact = quantity.Symbolic([exp(z), f; 0, exp(a*z)], 'grid', Z, 'name', 'Phi_A');
+
+ testCase.verifyTrue(numeric.near(F0.on() ./ F0.MAX(), PhiExact.on() ./ F0.MAX(), 1e-6));
+
+end
+%
+% function testPhiPsi(testCase)
+% N = 3;
+% z = sym('z', 'real');
+% Z = linspace(0, 1, 51)';
+% a = quantity.Symbolic(cat(3, ...
+% [ 0, 0, 0, 0; ...
+% 1 0 0 0; ...
+% 0 1 0 0; ...
+% 0 0 1 0], ...
+% [0 0 0, -3 + z; ...
+% 0 0 0 0; ...
+% 0 0 0 0; ...
+% 0 0 0 0], ...
+% [0 0 0, -4; ...
+% 0 0 0 0; ...
+% 0 0 0 0; ...
+% 0 0 0 0]...
+% ), 'name', 'A', 'grid', Z);
+% A = quantity.Operator(a);
+% B = quantity.Symbolic([1 1; zeros(3,2)], 'name', 'B', 'variable', z, 'grid', Z);
+% B = quantity.Operator(B);
+%
+% [Phi1, Psi1] = A.stateTransitionMatrix(N, B);
+% [Phi2, Psi2] = A.stateTransitionByOde(N, B);
+%
+% deltaPhi = Phi1.relativeErrorSupremum(Phi2);
+% deltaPsi = Psi1.relativeErrorSupremum(Psi2);
+%
+% %% TODO: write a test for the fundamental matrices.
+% gamma = ones(4,1);
+% nu = ones(2,1);
+%
+% for k = 1:N
+% wk_num(k,:) = Phi(:,:,k) * gamma + Psi(:,:,k) * nu;
+% end
+%
+% wk_ode = quantity.Discrete.empty(0,4);
+%
+% for k = 1:N
+% % odeStep = @(z, w, A, B, v) (A.on(z) * w + B.on(z) * v);
+% [~, tmp] = ode15s(@(z,w) odeStep(z, w, A, B, nu, wk_ode), Z, gamma * (k==1));
+% wk_ode(k, :) = quantity.Discrete(tmp, 'grid', Z, 'gridName', 'z');
+% end
+%
+% P = wk_num.relativeErrorSupremum(wk_ode);
+% P.plot
+%
+%
+% end
+
+%
function testFundamental(testCase)
-%%
-Z = sym('z', 'real');
-zeta = sym('zeta', 'real');
-A = quantity.Operator(cat(3, ...
- [ 0, 0, 0, 0; ...
- 1 0 0 0; ...
- 0 1 0 0; ...
- 0 0 1 0], ...
- [0 0 0, -3; ...
- 0 0 0 0; ...
- 0 0 0 0; ...
- 0 0 0 0], ...
- [0 0 0, -4; ...
- 0 0 0 0; ...
- 0 0 0 0; ...
- 0 0 0 0]...
- ), 'name', 'A');
-
-%%
-Nz = A.gridSize;
-n = size(A, 2);
-d = size(A, 3); % As order of a polynom is zero based
-zz = linspace(0, 1, 23)';
-
-F = misc.fundamentalMatrix.odeSolver( A(:,:,1).on(zz), zz);
-phi = quantity.Discrete( squeeze(F(:,1,:,:)), 'grid', {zz}, 'gridName', 'z', 'name', 'Phi');
-
-%%
-P = expm( A(:,:,1).at(1) * (Z - zeta) ) ;
-PHI0 = subs(P, zeta, 0);
-PHI = quantity.Symbolic(PHI0, 'grid', {zz}, 'figureID', 3);
+ Z = sym('z', 'real');
+ zeta = sym('zeta', 'real');
+ a0 = [ 0, 0, 0, 0; ...
+ 1, 0, 0, 0; ...
+ 0, 1, 0, 0; ...
+ 0, 0, 1, 0];
+ A0 = quantity.Discrete( reshape(a0, 1, 4, 4), 'gridName', '', 'grid', 0);
+
+ zz = linspace(0, 1, 51)';
+ A = quantity.Operator({A0}, 'grid', zz);
+ F0 = quantity.Discrete(squeeze(A.computePhi0().on({zz, 0})), 'grid', zz, 'gridName', 'z');
+
+ F = misc.fundamentalMatrix.odeSolver_par( A0.on(zz) , zz);
+ FF = quantity.Discrete(F, 'grid', {zz, zz}, 'gridName', {'z', 'zeta'});
+
+ F1 = quantity.Discrete(squeeze( F(:,1,:,:) ), 'grid', {zz}, 'gridName', {'z'}, 'name', 'Phi');
+
+ P = expm( squeeze(A0.on()) * (Z - zeta) ) ;
+ PHI0 = subs(P, zeta, 0);
+ F2 = quantity.Symbolic(PHI0, 'grid', {zz}, 'figureID', 3);
+
+ testCase.verifyTrue(numeric.near(F0.on(), F1.on(), 1e-7));
+ testCase.verifyTrue(numeric.near(F0.on(), F2.on(), 1e-7));
end
\ No newline at end of file
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index 4dc30e1..3320fc0 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -442,7 +442,11 @@ for k = 1:q.gridSize()
end
%%
-verifyEqual(testCase, magic700, q.on());
+testCase.verifyEqual(magic700, q.on());
+testCase.verifyTrue(q.isConstant());
+
+p = quantity.Symbolic([0 0 0 0 z], 'grid', linspace(0,1)', 'variable', z);
+testCase.verifyFalse(p.isConstant());
end
--
GitLab