Skip to content
GitLab
Commit b9fd6def authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

revised subs: works and all unittests are running.

parent 2cabaaab
Hide whitespace changes
Inline Side-by-side
+misc/Backstepping.m
View file @ b9fd6def
......@@ -24,39 +24,39 @@ classdef Backstepping < handle & matlab.mixin.Copyable
% char-arrays.
% 1. With the backstepping kernel:
% misc.Backstepping(...
% 'kernel', [], 'signOfIntegralTerm', +1, ...
% 'integralBounds', []);
% "kernel", [], "signOfIntegralTerm", +1, ...
% "integralBounds", []);
% 2. With the inverse backstepping kernel:
% misc.Backstepping(...
% 'kernelInverse', [], 'signOfIntegralTermInverse', +1, ...
% 'integralBounds', []);
% "kernelInverse", [], "signOfIntegralTermInverse", +1, ...
% "integralBounds", []);
% 3. With the both the backstepping kernel and its inverse:
% misc.Backstepping(...
% 'kernel', [], 'signOfIntegralTerm', +1, ...
% 'kernelInverse', [], 'signOfIntegralTermInverse', +1,
% 'integralBounds', []);
% "kernel", [], "signOfIntegralTerm", +1, ...
% "kernelInverse", [], "signOfIntegralTermInverse", +1,
% "integralBounds", []);
% In the cases 1 and 2, the other kernel (1: kernelInverse, 2. kernel)
% can be calculated easily by calling the method invert()
myParser = misc.Parser();
addParameter(myParser, 'kernel', []);
addParameter(myParser, 'kernelInverse', []);
addParameter(myParser, 'signOfIntegralTerm', []);
addParameter(myParser, 'signOfIntegralTermInverse', []);
addParameter(myParser, 'integralBounds', []);
addParameter(myParser, "kernel", []);
addParameter(myParser, "kernelInverse", []);
addParameter(myParser, "signOfIntegralTerm", []);
addParameter(myParser, "signOfIntegralTermInverse", []);
addParameter(myParser, "integralBounds", []);
myParser.parse(varargin{:});
if ~isempty(fieldnames(myParser.Unmatched))
warning('unmatched parameter in input parser');
warning("unmatched parameter in input parser");
end
assert(~isempty(myParser.Results.integralBounds), ...
'integralBounds must be specified');
"integralBounds must be specified");
if ~isempty(myParser.Results.kernel)
assert(~isempty(myParser.Results.signOfIntegralTerm), ...
'signOfIntegralTerm must be specified');
"signOfIntegralTerm must be specified");
end
if ~isempty(myParser.Results.kernelInverse)
assert(~isempty(myParser.Results.signOfIntegralTermInverse), ...
'signOfIntegralTermInverse must be specified');
"signOfIntegralTermInverse must be specified");
end
propertyNames = fieldnames(myParser.Results);
......@@ -72,18 +72,18 @@ classdef Backstepping < handle & matlab.mixin.Copyable
% xOut = xIn + obj.signOfIntegralTerm * ...
% int_integralBounds{1}^integralBounds{2}
% obj.kernel(z, zeta) * xIn(zeta) dzeta
assert(~isempty(obj.kernel), 'no kernel-property specified');
assert(size(xIn, 1) == size(obj.kernel, 2), 'xIn does not fit to kernel');
assert(~isempty(obj.kernel), "no kernel-property specified");
assert(size(xIn, 1) == size(obj.kernel, 2), "xIn does not fit to kernel");
if obj(1).signOfIntegralTerm == 1
xOut = xIn + ...
cumInt(obj(1).kernel * subs(xIn, {'z'}, {'zeta'}), ...
'zeta', obj(1).integralBounds{:});
cumInt(obj(1).kernel * subs(xIn, "z", "zeta"), ...
"zeta", obj(1).integralBounds{:});
elseif obj(1).signOfIntegralTerm == -1
xOut = xIn - ...
cumInt(obj(1).kernel * subs(xIn, {'z'}, {'zeta'}), ...
'zeta', obj(1).integralBounds{:});
cumInt(obj(1).kernel * subs(xIn, "z", "zeta"), ...
"zeta", obj(1).integralBounds{:});
else
error('signOfIntegralTerm is neither -1 nor +1')
error("signOfIntegralTerm is neither -1 nor +1")
end
end % transform
......@@ -92,18 +92,18 @@ classdef Backstepping < handle & matlab.mixin.Copyable
% xOut = xIn + obj.signOfIntegralTermInverse * ...
% int_integralBounds{1}^integralBounds{2}
% obj.kernelInverse(z, zeta) * xIn(zeta) dzeta
assert(~isempty(obj.kernelInverse), 'no kernelInverse-property specified')
assert(size(xIn, 1) == size(obj.kernel, 2), 'xIn does not fit to kernelInverse');
assert(~isempty(obj.kernelInverse), "no kernelInverse-property specified")
assert(size(xIn, 1) == size(obj.kernel, 2), "xIn does not fit to kernelInverse");
if obj(1).signOfIntegralTermInverse == 1
xOut = xIn + ...
cumInt(obj(1).kernelInverse * subs(xIn, {'z'}, {'zeta'}), ...
'zeta', obj(1).integralBounds{:});
cumInt(obj(1).kernelInverse * subs(xIn, "z", "zeta"), ...
"zeta", obj(1).integralBounds{:});
elseif obj(1).signOfIntegralTermInverse == -1
xOut = xIn - ...
cumInt(obj(1).kernelInverse * subs(xIn, {'z'}, {'zeta'}), ...
'zeta', obj(1).integralBounds{:});
cumInt(obj(1).kernelInverse * subs(xIn, "z", "zeta"), ...
"zeta", obj(1).integralBounds{:});
else
error('signOfIntegralTermInverse is neither -1 nor +1')
error("signOfIntegralTermInverse is neither -1 nor +1")
end
end % transformInverse()
......@@ -115,7 +115,7 @@ classdef Backstepping < handle & matlab.mixin.Copyable
end
xOriginal = quantity.Discrete(...
ones(101, size(obj.kernel, 1)), ...
quantity.Domain('z', linspace(0, 1, 101)));
quantity.Domain("z", linspace(0, 1, 101)));
xTilde = obj.transform(xOriginal);
xOriginalAgain = obj.transformInverse(xTilde);
myError = max(median(abs(xOriginal-xOriginalAgain)));
......@@ -156,32 +156,32 @@ classdef Backstepping < handle & matlab.mixin.Copyable
end
newKernel = quantity.Discrete(knownKernelWithSign);
if isequal(obj.integralBounds, {0, 'z'})
if isequal(obj.integralBounds, {0, "z"})
signOfReceprocityIntegral = +1;
elseif isequal(obj.integralBounds, {'z', 1})
elseif isequal(obj.integralBounds, {"z", 1})
signOfReceprocityIntegral = -1;
else
error('invalid integral bounds');
error("invalid integral bounds");
end
% successive approximation
progress = misc.ProgressBar('name', ...
progress = misc.ProgressBar("name", ...
"Successive Calculation of inverse " + ...
"of " + knownKernelWithSign(1).name + ": ", ...
'steps', 50, 'terminalValue', 50, 'printAbsolutProgress', true, 'silent', silent);
"steps", 50, "terminalValue", 50, "printAbsolutProgress", true, "silent", silent);
progress.start();
for it = 1 : progress.steps
newKernelLastIteration = newKernel;
newKernel = - knownKernelWithSign + signOfReceprocityIntegral * ...
cumInt(knownKernelWithSign.subs('zeta', 'eta') * newKernel.subs('z', 'eta'), ...
'eta', 'z', 'zeta');
cumInt(knownKernelWithSign.subs("zeta", "eta") * newKernel.subs("z", "eta"), ...
"eta", "z", "zeta");
changeInThisIteration = MAX(abs(newKernel-newKernelLastIteration));
if changeInThisIteration < tolerance
break;
end
progress.raise(it, 'addMessage', ...
['change = ', num2str(changeInThisIteration)]);
progress.raise(it, "addMessage", ...
"change = " + num2str(changeInThisIteration));
end
newKernel = newIntegralSign * newKernel;
progress.stop();
......@@ -207,27 +207,27 @@ classdef Backstepping < handle & matlab.mixin.Copyable
% boundary z=zeta.
%
% INPUT PARAMETER
% type 'kernel' or 'kernelInverse'
% type "kernel" or "kernelInverse"
% myGrid vector of discrete grid on which derivatives are
% calcualted
% OUTPUT PARAMETER
% K_dzeta spatial derivative in zeta direction as quantity.Discrete
% K_dz spatial derivative in z direction as quantity.Discrete
% set parameter domain which must be either 'zeta<=z' or 'zeta>=z'
if isequal(obj.integralBounds, {0, 'z'})
domain = 'zeta<=z';
elseif isequal(obj.integralBounds, {'z', 1})
domain = 'zeta>=z';
% set parameter domain which must be either "zeta<=z" or "zeta>=z"
if isequal(obj.integralBounds, {0, "z"})
domain = "zeta<=z";
elseif isequal(obj.integralBounds, {"z", 1})
domain = "zeta>=z";
end
% select kernel according to input
if strcmp(type, 'kernelInverse')
if strcmp(type, "kernelInverse")
thisKernel = obj.kernelInverse;
elseif strcmp(type, 'kernel')
elseif strcmp(type, "kernel")
thisKernel = obj.kernel;
else
error("type must be either 'kernel' or 'kernelInverse'");
error("type must be either kernel or kernelInverse");
end
% select default grid
......@@ -236,7 +236,7 @@ classdef Backstepping < handle & matlab.mixin.Copyable
end
spacing = myGrid(2) - myGrid(1);
assert(numeric.near(diff(myGrid), spacing, 1e-6), ...
'grid must be homogenious');
"grid must be homogenious");
% calculate gradient numerically
quantityDomain = [quantity.Domain("z", myGrid), quantity.Domain("zeta", myGrid)];
......@@ -246,11 +246,11 @@ classdef Backstepping < handle & matlab.mixin.Copyable
% create quantities as output parameters
K_dzeta = quantity.Discrete(K_dzetaMat, ...
quantityDomain , ...
'name', "d_{zeta}" + thisKernel(1).name);
"name", "d_{zeta}" + thisKernel(1).name);
K_dz = quantity.Discrete(K_dzMat, ...
quantityDomain , ...
'name', "d_z" + thisKernel(1).name);
"name", "d_z" + thisKernel(1).name);
end % gradient()
function domainSelector = get.domainSelector(obj)
......@@ -276,10 +276,10 @@ classdef Backstepping < handle & matlab.mixin.Copyable
% get domainSelector depending on obj.integralBounds
if isequal(bounds, {0, 1})
domainSelector = quantity.Discrete.ones([1, 1], domain);
elseif isequal(bounds, {0, "z"}) || isequal(bounds, {0, 'z'})
elseif isequal(bounds, {0, "z"}) || isequal(bounds, {0, "z"})
ndGrids = domain.ndgrid;
domainSelector = quantity.Discrete(double(ndGrids{1} >= ndGrids{2}), domain);
elseif isequal(bounds, {"z", 1}) || isequal(bounds, {'z', 1})
elseif isequal(bounds, {"z", 1}) || isequal(bounds, {"z", 1})
ndGrids = domain.ndgrid;
domainSelector = quantity.Discrete(double(ndGrids{1} <= ndGrids{2}), domain);
elseif isnumeric(bounds{1}) && isnumeric(bounds{2})
......
+misc/solveVolterraIntegroDifferentialIVP.m
View file @ b9fd6def
......@@ -72,7 +72,7 @@ if NameValueArgs.D.abs.MAX() >= 1e3*eps
% is used.
N = misc.solveVolterraIntegroDifferentialIVP(z, ic, ...
'M', NameValueArgs.M, 'A', NameValueArgs.A, 'B', NameValueArgs.B, ...
'C', NameValueArgs.C + cumInt(NameValueArgs.D * N.subs('z', 'zeta'), 'zeta', 0, 'z'));
'C', NameValueArgs.C + cumInt(NameValueArgs.D * N.subs("z", "zeta"), "zeta", 0, "z"));
end
else
% the integral term can be neclegted, hence the solution is easily obtained with
......@@ -87,7 +87,7 @@ else
z.grid.', ...
ic(:));
N = quantity.Discrete(reshape(myN, [z.n, n, m]), z, 'name', "N");
N = quantity.Discrete(reshape(myN, [z.n, n, m]), z, "name", "N");
end
......
+signals/PolynomialOperator.m
View file @ b9fd6def
......@@ -363,12 +363,12 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
[quantity.Domain("z", spatialDomain), quantity.Domain("zeta", spatialDomain)]);
b = quantity.Discrete(B_full, quantity.Domain("zeta", spatialDomain));
p = cumInt(f * b, 'zeta', spatialDomain(1), 'z');
p = cumInt(f * b, "zeta", spatialDomain(1), "z");
z0 = A(1).grid{1}(1);
for k = 1:N
idxk = idx + (k-1)*n;
Phi(:,:,k) = f(idxk, idx).subs('zeta', z0);
Phi(:,:,k) = f(idxk, idx).subs("zeta", z0);
Psi(:,:,k) = p(idxk,:);
end
......@@ -441,16 +441,15 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
F0 = ip.Results.F0;
end
z = F0(1).grid{1};
Phi(:,:,1) = F0.subs('zeta', 0);
Phi(:,:,1) = F0.subs("zeta", 0);
if isempty(B)
Psi(:,:,1) = quantity.Discrete.empty();
else
Psi(:,:,1) = cumInt( F0 * B(1).coefficient.subs('z', 'zeta'), ...
'zeta', z(1), 'z');
Psi(:,:,1) = cumInt( F0 * B(1).coefficient.subs("z", "zeta"), "zeta", z(1), "z");
end
On = quantity.Discrete.zeros([n n], z, 'gridName', 'z');
Om = quantity.Discrete.zeros([n m], z, 'gridName', 'z');
On = quantity.Discrete.zeros([n n], z, "gridName", "z");
Om = quantity.Discrete.zeros([n m], z, "gridName", "z");
for k = 2 : N
pbar.raise();
......@@ -468,11 +467,9 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
% compute the integration of the fundamental matrix with the temporal
% values:
% Phi_k(z) = int_0_z Phi_0(z, zeta) * M(zeta) d_zeta
Phi(:, :, k) = cumInt(F0 * sumPhi.subs('z', 'zeta') , ...
'zeta', z(1), 'z');
Phi(:, :, k) = cumInt(F0 * sumPhi.subs("z", "zeta"), "zeta", z(1), "z");
Psi(:, :, k) = cumInt(F0 * sumPsi.subs('z', 'zeta'), ...
'zeta', z(1), 'z');
Psi(:, :, k) = cumInt(F0 * sumPsi.subs("z", "zeta"), "zeta", z(1), "z");
end
%
pbar.stop();
......
+signals/WhiteGaussianNoise.m
View file @ b9fd6def
......@@ -103,7 +103,7 @@ classdef WhiteGaussianNoise
function P = cumulativeDistributionFunction(obj, varargin)
p = obj.probabilityDensityFunction(varargin{:});
P = p.cumInt('x', p(1).grid{1}(1), 'x');
P = p.cumInt("x", p(1).grid{1}(1), "x");
end
function p = probability(obj, z)
......
+unittests/+misc/testBackstepping.m
View file @ b9fd6def
......@@ -38,29 +38,29 @@ function testInvert(testCase)
% init random kernel
syms z zeta
myGrid = linspace(0, 1, 21);
myDomain = [quantity.Domain('z', myGrid), quantity.Domain('zeta', myGrid)];
myDomain = [quantity.Domain("z", myGrid), quantity.Domain("zeta", myGrid)];
kernelData = quantity.Discrete(quantity.Symbolic(...
[z^2/2, sin(z); 1/(2-z), -z] * magic(2) * [-zeta^2/2, sin(2*zeta); 1/(2-zeta), zeta], ...
myDomain, 'name', 'K'));
myDomain, "name", "K"));
% set up kernel data in misc.Backstepping objects
k(1) = misc.Backstepping('kernel', kernelData, ...
'signOfIntegralTerm', +1, 'integralBounds', {0, 'z'});
k(2) = misc.Backstepping('kernel', kernelData, ...
'signOfIntegralTerm', -1, 'integralBounds', {0, 'z'});
k(3) = misc.Backstepping('kernel', kernelData, ...
'signOfIntegralTerm', +1, 'integralBounds', {'z', 1});
k(4) = misc.Backstepping('kernel', kernelData, ...
'signOfIntegralTerm', -1, 'integralBounds', {'z', 1});
k(1) = misc.Backstepping("kernel", kernelData, ...
"signOfIntegralTerm", +1, "integralBounds", {0, "z"});
k(2) = misc.Backstepping("kernel", kernelData, ...
"signOfIntegralTerm", -1, "integralBounds", {0, "z"});
k(3) = misc.Backstepping("kernel", kernelData, ...
"signOfIntegralTerm", +1, "integralBounds", {"z", 1});
k(4) = misc.Backstepping("kernel", kernelData, ...
"signOfIntegralTerm", -1, "integralBounds", {"z", 1});
% inverse
k(5) = misc.Backstepping('kernelInverse', kernelData, ...
'signOfIntegralTermInverse', +1, 'integralBounds', {0, 'z'});
k(6) = misc.Backstepping('kernelInverse', kernelData, ...
'signOfIntegralTermInverse', -1, 'integralBounds', {0, 'z'});
k(7) = misc.Backstepping('kernelInverse', kernelData, ...
'signOfIntegralTermInverse', +1, 'integralBounds', {'z', 1});
k(8) = misc.Backstepping('kernelInverse', kernelData, ...
'signOfIntegralTermInverse', -1, 'integralBounds', {'z', 1});
k(5) = misc.Backstepping("kernelInverse", kernelData, ...
"signOfIntegralTermInverse", +1, "integralBounds", {0, "z"});
k(6) = misc.Backstepping("kernelInverse", kernelData, ...
"signOfIntegralTermInverse", -1, "integralBounds", {0, "z"});
k(7) = misc.Backstepping("kernelInverse", kernelData, ...
"signOfIntegralTermInverse", +1, "integralBounds", {"z", 1});
k(8) = misc.Backstepping("kernelInverse", kernelData, ...
"signOfIntegralTermInverse", -1, "integralBounds", {"z", 1});
% invert all
for it = 1 : numel(k)
......@@ -70,9 +70,10 @@ end
% compare to old implementation for int_0^z
myDomain = quantity.Discrete(quantity.Symbolic(zeta<z, myDomain));
zZetaDomain = [quantity.Domain('z', myGrid), quantity.Domain('zeta', myGrid)];
Kref = quantity.Discrete(...
misc.invertBacksteppingKernel(kernelData.on({myGrid, myGrid}, {'z', 'zeta'}), -1), ...
[quantity.Domain('z', myGrid), quantity.Domain('zeta', myGrid)]);
misc.invertBacksteppingKernel(kernelData.on(zZetaDomain), -1), ...
zZetaDomain);
k1Diff = myDomain * (Kref - k(1).kernelInverse);
k5Diff = myDomain * (Kref - k(5).kernel);
testCase.verifyEqual(MAX(abs(k1Diff)), 0, 'AbsTol', 2e-4);
......
+unittests/+misc/testSolveVolterraIntegroDifferentialIVP.m
View file @ b9fd6def
......@@ -23,7 +23,7 @@ classdef testSolveVolterraIntegroDifferentialIVP < matlab.unittest.TestCase
tc.C = quantity.Symbolic([1+sym("z"), 0, -1; 0, 2, -1], tc.z);
tc.D = quantity.Symbolic([1+sym("z"), sym("z")*sym("zeta"); sym("zeta"), 2], [tc.z, tc.zeta]);
tc.N = misc.solveVolterraIntegroDifferentialIVP(tc.z, tc.ic, ...
'M', tc.M, 'A', tc.A, 'B', tc.B, 'C', tc.C, 'D', tc.D);
"M", tc.M, "A", tc.A, "B", tc.B, "C", tc.C, "D", tc.D);
end
end % methods(TestClassSetup)
......@@ -31,12 +31,12 @@ classdef testSolveVolterraIntegroDifferentialIVP < matlab.unittest.TestCase
function testSolution(tc)
% verify initial condition
tc.verifyEqual(tc.ic, tc.N.at(0), 'AbsTol', 10*eps);
tc.verifyEqual(tc.ic, tc.N.at(0), "AbsTol", 10*eps);
% verify ODE
odeResiduum = tc.M * tc.N.diff("z", 1) + tc.A * tc.N + tc.N * tc.B + tc.C ...
+ cumInt(tc.D * tc.N.subs("z", "zeta"), 'zeta', 0, 'z');
tc.verifyEqual(max(odeResiduum.abs.median(), [], 'all'), 0, 'AbsTol', 1e-4);
+ cumInt(tc.D * tc.N.subs("z", "zeta"), "zeta", 0, "z");
tc.verifyEqual(max(odeResiduum.abs.median(), [], "all"), 0, "AbsTol", 1e-4);
% plot
if tc.doPlot
......@@ -53,13 +53,13 @@ classdef testSolveVolterraIntegroDifferentialIVP < matlab.unittest.TestCase
for it = 1 : numel(gridPoints)
Ncell{it} = misc.solveVolterraIntegroDifferentialIVP(...
quantity.Domain("z", linspace(0, 1, gridPoints(it))), ...
tc.ic, 'M', tc.M, 'A', tc.A, 'B', tc.B, 'C', tc.C, 'D', tc.D);
medianDifference(it) = max(median(abs(tc.N - Ncell{it})), [], 'all');
tc.ic, "M", tc.M, "A", tc.A, "B", tc.B, "C", tc.C, "D", tc.D);
medianDifference(it) = max(median(abs(tc.N - Ncell{it})), [], "all");
end % for it = 1 : numel(gridPoints)
% decay
figure()
plot(gridPoints, medianDifference, 'x');
plot(gridPoints, medianDifference, "x");
% comparison
figure();
......
+unittests/+signals/testPolynomialOperator.m
View file @ b9fd6def
......@@ -108,7 +108,7 @@ end
function testChangeDomain(testCase)
A = testCase.TestData.A;
a = testCase.TestData.a;
A0 = A.subs('z', 0);
A0 = A.subs("z", 0);
for k = 1:3
testCase.verifyEqual( ...
......
+unittests/+signals/testTimeDelayOperator.m
View file @ b9fd6def
......@@ -67,7 +67,7 @@ function testApplyTo(testCase)
z = testCase.TestData.z;
t = testCase.TestData.t;
tau = quantity.Domain('tau', linspace(-1, 3, 201));
tau = quantity.Domain("tau", linspace(-1, 3, 201));
v = testCase.TestData.delay.variable;
h = quantity.Discrete( sin(tau.grid * 2 * pi), tau );
......@@ -80,7 +80,7 @@ function testApplyTo(testCase)
d2 = testCase.TestData.delay.spatialDomain2;
H2 = d2.applyTo(h);
H2 = subs(H2, 'zeta', 0);
H2 = subs(H2, "zeta", 0);
testCase.verifyEqual(H.on(), H2.on(), 'AbsTol', 5e-3);
%% test a diagonal operator
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment