diff --git a/+misc/+fundamentalMatrix/rudolphWoittennekNumeric.m b/+misc/+fundamentalMatrix/rudolphWoittennekNumeric.m
index 6c11f8fd73623cb5a3379022a56a8e0470fd9064..c8693e71665082f41a08618b68fc901d5c0eb55e 100644
--- a/+misc/+fundamentalMatrix/rudolphWoittennekNumeric.m
+++ b/+misc/+fundamentalMatrix/rudolphWoittennekNumeric.m
@@ -43,16 +43,12 @@ pbar = misc.ProgressBar(...
F0 = A.fundamentalMatrix0(z);
-Phi(:,:,1) = quantity.Discrete( F0.subs('zeta', 0) );
+Phi(:,:,1) = F0.subs('zeta', 0);
B = B.subs('z', 'zeta');
A = A.subs('z', 'zeta');
Psi(:,:,1) = F0.volterra(B);
-% TODO Implement the intPhiM function based on quantity class
-
-
-
for kN = 2 : N
pbar.raise();
@@ -63,8 +59,8 @@ for kN = 2 : N
% compute the temporal matrices:
for i = 2 : dA
%M = A_i(z) * F_{kN-i+1}(z, xi); xi = 0;
- M = M + A(:, :, i) * Phi(:, :, kN-i+1);
- N = N + A(:, :, i) * Psi(:, :, kN-i+1);
+ M = M + A(:, :, i) * Phi(:, :, kN-i+1).subs('z', 'zeta');
+ N = N + A(:, :, i) * Psi(:, :, kN-i+1).subs('z', 'zeta');
end
if size(B, 3) >= kN
N = N + B(:, :, kN);
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 4aae3ae791f72884671ce9d1d76045ee468d9bf5..f2fbb7734b2f3e2abe111d78b2d2da2eb1c99b6e 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -1352,7 +1352,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if ip.isDefault('variable') && x.nargin > 1
error(['For x quantities wrt. two independent' ...
' variables, the name of the variable has to be' ...
- 'parametrized explicitly']);
+ ' parametrized explicitly']);
end
diff --git a/+quantity/Operator.m b/+quantity/Operator.m
index b0bededdaba6767e93b734cc68681a6bae3c2648..ceae683f5ee27f65fa780a74e2ba282f6278678e 100644
--- a/+quantity/Operator.m
+++ b/+quantity/Operator.m
@@ -72,7 +72,7 @@ classdef Operator < quantity.Symbolic
[Phi, F, Psi, P] = ...
misc.fundamentalMatrix.rudolphWoittennek(obj.sym, ...
'N', N, 'zeta', obj(1).grid{1}(1), 's', obj(1).operatorVar, ...
- 'B', B);
+ 'B', B.sym);
if isempty(B)
Psi = zeros(size(B));
@@ -84,9 +84,8 @@ classdef Operator < quantity.Symbolic
function [Phi, Psi] = fundamentalMatrixNumeric(obj, varargin)
ip = misc.Parser();
ip.addOptional('N', 5);
- ip.addOptional('B', quantity.Discrete.zeros([size(obj,1), 1], ...
- obj(1).grid) );
- ip.addOptional('spatialDomain', []);
+ ip.addOptional('B', quantity.Discrete.empty([size(obj,1), 0]));
+ ip.addOptional('spatialDomain', obj(1).grid{:});
ip.parse(varargin{:});
misc.struct2ws(ip.Results)
@@ -122,7 +121,7 @@ classdef Operator < quantity.Symbolic
f0 = expm(obj(:,:,1).on()*(z - zeta));
F0 = quantity.Symbolic(f0, 'grid', {spatialDomain, spatialDomain});
else
- f0 = misc.fundamentalMatrix.odeSolver_par( Adisc, spatialDomain);
+ f0 = misc.fundamentalMatrix.odeSolver_par( Adisc, spatialDomain );
F0 = quantity.Discrete(f0, ...
'size', [size(obj, 1), size(obj, 2)], ...
'gridName', {'z', 'zeta'}, ...
diff --git a/+quantity/Symbolic.m b/+quantity/Symbolic.m
index 069e044c1832a621e34a9f0d50b64b67ac9932bd..8d0b6a21564757a868060de381e800fc473e82ba 100644
--- a/+quantity/Symbolic.m
+++ b/+quantity/Symbolic.m
@@ -22,7 +22,8 @@ classdef Symbolic < quantity.Function
if nargin > 0
variableParser = misc.Parser();
variableParser.addParameter('variable', quantity.Symbolic.getVariable(valueContinuous));
- variableParser.parse2ws(varargin{:});
+ variableParser.parse(varargin{:});
+ variable = variableParser.Results.variable;
numVar = numel(variable);
defaultGrid = cell(1, numVar);
for it = 1:numVar
@@ -35,7 +36,6 @@ classdef Symbolic < quantity.Function
gridParser = misc.Parser();
gridParser.addParameter('grid', defaultGrid);
gridParser.parse(varargin{:});
- misc.struct2ws(gridParser.Results);
s = size(valueContinuous);
S = num2cell(s);
@@ -78,8 +78,8 @@ classdef Symbolic < quantity.Function
error('If the gridName(s) are set explicitly, they have to be the same as the variable names!')
end
- parentVarargin = {fun, varargin{:}, 'grid', gridParser.Results.grid, ...
- 'gridName', gridName};
+ parentVarargin = [{fun}, varargin(:)', {'grid'}, {gridParser.Results.grid}, ...
+ {'gridName'}, {gridName}];
end
% call parent class
obj@quantity.Function(parentVarargin{:});
diff --git a/+unittests/+quantity/testOperator.m b/+unittests/+quantity/testOperator.m
index 27bf9807cd144673b4809df5645f6f3c97687168..46b9acbdb183c763dc7c30292dd345385884a23b 100644
--- a/+unittests/+quantity/testOperator.m
+++ b/+unittests/+quantity/testOperator.m
@@ -42,13 +42,13 @@ 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);
+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);
+PsiSym = quantity.Operator(P, 'grid', Z, 'variable', z);
% plot(PhiSym, 'fig', 4)
@@ -58,19 +58,15 @@ if numeric.near(a, 1)
else
f = (exp(z) - exp(a * z) - (1 - a) * z * exp( a * z ) ) / (1 - a)^2;
end
-% plot(PhiNum, 'fig', 1);
-
-PhiA = quantity.Symbolic([exp(z), f; 0, exp(a*z)], 'grid', Z, 'name', 'Phi_A');
-% plot(PhiA, 'fig', 2);
-[D, sup] = relativeErrorSupremum(PhiNum, PhiA);
+PhiExact = quantity.Symbolic([exp(z), f; 0, exp(a*z)], 'grid', Z, 'name', 'Phi_A');
-% plot(D, 'fig', 3);
-
-%%
-F0 = A.fundamentalMatrix0;
-F0 = squeeze(F0(:,1,:,:));
-F0 = quantity.Discrete(quantity.Discrete.value2cell(F0, length(Z), [2, 2]), 'grid', Z);
+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);
end
@@ -95,34 +91,25 @@ A = quantity.Operator(cat(3, ...
0 0 0 0]...
), 'name', 'A', 'grid', Z);
-B = quantity.Operator([1 1; zeros(3,2)], 'name', 'B', 'grid', Z);
+B = quantity.Operator([1 1; zeros(3,2)], 'name', 'B', 'variable', z, 'grid', Z);
-[PhiNum, PsiNum] = A.fundamentalMatrixNumeric(N, B);
+[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);
-%
-% plot(D);
-%
-% figure(1); clf;
-% plot(1:N, max(D));
-% title(num2str(length(Z)));
+delta = D.on();
+testCase.verifyLessThan(max(delta(:)), 1e-2);
+D = relativeErrorSupremum(PsiNum, PsiSym);
+delta = D.on();
+testCase.verifyLessThan(max(delta(:)), 1e-2);
-%%
-verifyLessThan(testCase, max(max(D)), 1e-2);
-% #todo:
-% 2) implement the computation of the multiplication for polynomial matrices
-% 3) implement the computation of the determinant for polynomial matrices
end
@@ -150,13 +137,10 @@ A = quantity.Operator(cat(3, ...
Nz = A.gridSize;
n = size(A, 2);
d = size(A, 3); % As order of a polynom is zero based
-zz = A(1).grid{1};
-
-F = misc.fundamentalMatrix.odeSolver( A(:,:,1).on(), A(1).grid{1});
-phi = quantity.Discrete(quantity.Discrete.value2cell(squeeze(F(:,1,:,:)), [Nz], [n, n]), 'grid', {zz}, 'name', 'Phi');
-% plot(int(phi))
-
+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) ) ;