diff --git a/+export/Data.m b/+export/Data.m
index a5e3c03b17e97bf7bf80319605d7f7e229ec19d2..86723d6096ae932a29fad7ced0202d06d6386aab 100644
--- a/+export/Data.m
+++ b/+export/Data.m
@@ -3,9 +3,9 @@ classdef (Abstract) Data < handle & matlab.mixin.Heterogeneous
% Detailed explanation goes here
properties
- basepath = '.';
- foldername = '.';
- filename = 'data';
+ basepath char = '.';
+ foldername char = '.';
+ filename char = 'data';
N = 201;
end
@@ -53,7 +53,7 @@ classdef (Abstract) Data < handle & matlab.mixin.Heterogeneous
for k = 1:length(obj)
if ~obj(k).isempty
- if ~isdir(obj(k).absolutepath)
+ if ~isfolder(obj(k).absolutepath)
mkdir(obj(k).absolutepath);
end
diff --git a/+fractional/derivative.m b/+fractional/derivative.m
new file mode 100644
index 0000000000000000000000000000000000000000..d99fc000ac694d55dd10521b07af6373ccfe94cb
--- /dev/null
+++ b/+fractional/derivative.m
@@ -0,0 +1,39 @@
+function [h, t] = derivative(y, p, t0, t, optArg)
+% DERIVATIVE calculates a fractional derivative
+% [h, t] = derivative(y, p, t0, t, optArg) comput the p-th fractional derivative of a symbolic
+% expression y. The lower terminal is t0 and the symbolic independent variable is t. By the
+% name-value-pair argument, the "type" of the fractional derivative can be specified. Use "RL" for
+% the Riemann-Liouville derivative and "C" for the Caputo derivative.
+arguments
+ y sym
+ p (1,1) double {mustBePositive};
+ t0 (1,1) double;
+ t (1,1) sym;
+ optArg.type = "RL"
+end
+
+if p < 0
+ error('a must be positive! Use fintegral instead?');
+end
+
+if misc.nearInt(p)
+ h = diff(y, round(p));
+ return
+end
+
+n = ceil(p);
+tau = sym( string(t) + "_");
+
+if optArg.type == "C"
+ h = 1/gamma(n-p) * int((t-tau)^(n-p-1) * subs(diff(y, t, n), t, tau), tau, t0, t);
+ %DaY = 1/gamma(1+n-a) * ((t - t0)^(n-a) * subs( diff(y, n), t, t0) ...
+ % + int((t-tau)^(n-a) * subs(diff(y, 1+n), t, tau), tau, t0, t)); % Integration by parts
+
+elseif optArg.type == "RL"
+ h = 1 / gamma(n-p) * diff( int( (t-tau)^(n-p-1) * subs(y, t, tau), tau, t0, t), t, n);
+else
+ error("The derivative of the type " + optArg.type + " is not implemented. Only C for Caputo and RL for Riemann-Liouville are available")
+end
+
+
+end
\ No newline at end of file
diff --git a/+fractional/monomialFractionalDerivative.m b/+fractional/monomialFractionalDerivative.m
new file mode 100644
index 0000000000000000000000000000000000000000..23689944550178b97f5a17d43a20dfb103a84361
--- /dev/null
+++ b/+fractional/monomialFractionalDerivative.m
@@ -0,0 +1,25 @@
+function [y, m] = monomialFractionalDerivative(t, t0, nu, p, optArg)
+%MONOMIALFRACTIONALDERIVATIVE compute fractional derivative of a monomial
+% [y, m] = monomialFractionalDerivative(t, a, nu, p, optArg) compute the p-fractional derivative of
+% the monomial m = (t-a)^nu. At this, t must be a symbolic variable, t0 is the lower terminal as
+% double value, nu is the exponent as double, p is the fractional order as double the type of the
+% derivative can be specified by the name-value-pair "type". So far online the "RL" :
+% Riemann-Liouville derivative is implemented.
+arguments
+ t (1,1) sym; % symbolic variable for the monomial
+ t0 (1,1) double; % expansion point of the monome
+ nu (1,1) double; % order of the monome
+ p (1,1) double; % order of the fractional derivative
+ optArg.type string = "RL"; % Type of the fractional derivative
+end
+
+m = (t-t0)^nu;
+
+switch optArg.type
+
+ case "RL"
+ y = gamma( 1+nu ) / gamma( 1+nu-p ) * ( t - t0 )^(nu-p);
+ otherwise
+ error("Not jet implemented, only the RL - Rieman Liouville derivative is available")
+end
+end
\ No newline at end of file
diff --git a/+fractional/simulate.m b/+fractional/simulate.m
new file mode 100644
index 0000000000000000000000000000000000000000..c842f73f95f675911ddbafb42f3098c0acf2022f
--- /dev/null
+++ b/+fractional/simulate.m
@@ -0,0 +1,47 @@
+function [y, w, x] = simulate(stsp, u, optArgs)
+
+arguments
+ stsp ss;
+ u quantity.Discrete;
+ optArgs.timeDomain quantity.EquidistantDomain = u(1).domain;
+ optArgs.x0 = zeros(size(stsp.A, 1), 1);
+ optArgs.stepSize (1,1) double = u(1).domain.stepSize;
+ optArgs.tol = 1.0e-6;
+ optArgs.itmax = 500;
+ optArgs.spatialDomain = quantity.Domain("z", linspace(0, 1, size(stsp.A, 1)));
+ optArgs.silent (1,1) logical = false;
+end
+
+
+if ~exist('MT_FDE_PI1_Im.m', 'file')
+ error("The fractional solver MT_FDE_PI1_Im must be in the matlab path. It can be downloaded from https://de.mathworks.com/matlabcentral/fileexchange/66603-solving-multiterm-fractional-differential-equations-fde");
+end
+
+if isempty( stsp.E )
+ E = eye(size(stsp.A,1));
+else
+ E = stsp.E;
+end
+
+assert( stsp.UserData.alpha <= 1 && stsp.UserData.alpha > 0, ...
+ 'Fractional order must hold 0 < alpha <= 1!');
+
+f_fun = @(t, x) stsp.A * x + stsp.B * u.at(t);
+J_fun = @(t, x) stsp.A; % Jacobian matrix of f_fun, calculated by hand. It always has the same structure.
+
+[~, x] = FDE_PI1_Im_Ver_10(stsp.UserData.alpha, f_fun, J_fun, ...
+ optArgs.timeDomain.lower, ...
+ optArgs.timeDomain.upper, ...
+ optArgs.x0, ...
+ optArgs.stepSize, ...
+ [], ...
+ optArgs.tol, ...
+ optArgs.itmax, ....
+ optArgs.silent);
+
+w = quantity.Discrete(x, [optArgs.spatialDomain, optArgs.timeDomain]);
+x = quantity.Discrete(x', optArgs.timeDomain);
+
+y = stsp.C * x + stsp.D * u;
+
+end
\ No newline at end of file
diff --git a/+misc/+ss/combineInputSignals.m b/+misc/+ss/combineInputSignals.m
index 5dc153c73603161083e77e3cdd65a393b8bdd526..140e411c468e11125f6d38d392a7a8b10d07ec21 100644
--- a/+misc/+ss/combineInputSignals.m
+++ b/+misc/+ss/combineInputSignals.m
@@ -38,7 +38,7 @@ myParser = misc.Parser();
for it = 1 : numel(myInputNameUnique)
myParser.addParameter(myInputNameUniqueValidFieldName{it}, ...
quantity.Discrete(zeros(numel(t), myInputLength(it)), ...
- 'grid', t, 'gridName', 't', 'name', myInputNameUnique{it}), ...
+ quantity.Domain("t", t), 'name', myInputNameUnique{it}), ...
@(v) isa(v, 'quantity.Discrete') & isvector(v) & (numel(v)==myInputLength(it)));
end
diff --git a/+misc/+ss/planPolynomialTrajectory.m b/+misc/+ss/planPolynomialTrajectory.m
index 92dab9ef5a148e3210251b284c08d69b096b7b3b..1414951f805608618e5423ec0f072f3b9d170336 100644
--- a/+misc/+ss/planPolynomialTrajectory.m
+++ b/+misc/+ss/planPolynomialTrajectory.m
@@ -2,7 +2,7 @@ function [u, y, x] = planPolynomialTrajectory(sys, t, optArgs)
%FLATKERNELCOMPUTATION solves the transition x(0) -> x(T) for a given
%state space based on a flat output
%
-% [u,x,Cx] = FLATKERNELCOMPUTATION(system, xT, T) calculates
+% [u,x,Cx] = planPolynomialTrajectory(system, xT, T) calculates
% the input u, the states x and the output Cx for a given state space
% system symbolically. The calculation is based on the flat output based approach
% for trajectory planning. The system for which the
diff --git a/+misc/+ss/planTrajectory.m b/+misc/+ss/planTrajectory.m
index d252e45e0de22e6a349381cc1a526a62cef67227..1cfb386dcc68d95a2428a036ad47c7066a0c984a 100644
--- a/+misc/+ss/planTrajectory.m
+++ b/+misc/+ss/planTrajectory.m
@@ -28,18 +28,19 @@ x1 = optArgs.x1;
t = t(:);
t0 = t(1);
t1 = t(end);
+tDomain = quantity.Domain({optArgs.domainName}, t);
%% Compute the state transition matrix: Phi(t,tau) = expm(A*(t-tau) )
% Att0 = sys.A * quantity.Symbolic( sym('t') - sym('tau'), 'grid', {t, t}, 'gridName', {'t', 'tau'} );
-Phi_t0 = expm(sys.A * quantity.Discrete( t - t0, 'grid', {t}, 'gridName', {optArgs.domainName} ));
+Phi_t0 = expm(sys.A * quantity.Discrete( t - t0, tDomain));
-invPhi_t1 = expm(sys.A * quantity.Discrete( t1 - t, 'grid', {t}, 'gridName', {optArgs.domainName} ));
+invPhi_t1 = expm(sys.A * quantity.Discrete( t1 - t, tDomain));
%% compute the gramian controllability matrix
% W1 = misc.ss.gramian(sys, t, t0, optArgs.domainName);
% int_t0^t1 expm(A*(tau-t0)) * b * b^T * expm(A^T(tau-t0)) dtau
-W1 = cumInt( Phi_t0 * sys.b * sys.b' * expm(sys.A' * quantity.Discrete( t - t0, 'grid', {t}, 'gridName', {optArgs.domainName} )), ...
+W1 = cumInt( Phi_t0 * sys.b * sys.b' * expm(sys.A' * quantity.Discrete( t - t0, tDomain)), ...
Phi_t0(1).domain, t0, optArgs.domainName);
W1_t1 = W1.at(t1);
diff --git a/+misc/+ss/setPointChange.m b/+misc/+ss/setPointChange.m
new file mode 100644
index 0000000000000000000000000000000000000000..21084271049a6230626a57b23dd97987e920b388
--- /dev/null
+++ b/+misc/+ss/setPointChange.m
@@ -0,0 +1,47 @@
+function [f, a] = setPointChange(phi0, phi1, t0, t1, n, optArgs)
+%SETPOINTCHANGE polynomial fit of an initial and end value with a certain smoothness
+%
+% [f, a] = setPointChange(PHI0, PHI1, t0, t1, n, optArgs) plans a function f which connects the
+% initial value PHI0 with the end value PHI1, so that f(t0) = PHI0 and f(t1) = PHI1 as well as
+% d^k / dt^k f(t)|_{t=t0,t1} = 0, k = 1, 2, ..., n
+% The name of the variable can be set with "var" as name value pair.
+arguments
+ phi0 (1,1) double;
+ phi1 (1,1) double;
+ t0 (1,1) double;
+ t1 (1,1) double;
+ n (1,1) double;
+ optArgs.var = sym( "t" );
+end
+
+
+N = 2*n+2;
+
+M0 = zeros(n+1, N);
+M1 = zeros(n+1, N);
+
+coefficient = @(k, i, t) factorial( i + k -1 ) / factorial( i - 1) * t^(i-1);
+
+for k = 0:n
+ for i = 1:( N - k )
+ M0(k+1, i + k) = coefficient(k, i, t0);
+ M1(k+1, i + k) = coefficient(k, i, t1);
+ end
+end
+
+a = [M0; M1] \ [phi0; zeros(n,1); phi1; zeros(n,1)];
+
+if misc.issym( optArgs.var )
+ t = optArgs.var;
+else
+ t = sym( optArgs.var );
+end
+
+f = sym(0);
+
+for i = 1 : N
+ f = f + a(i) * t^(i-1);
+end
+
+end
+
diff --git a/+misc/+ss/simulationOutput2Quantity.m b/+misc/+ss/simulationOutput2Quantity.m
index bc2efcc054f2a735052bf6fa7aa0d9d7b806c026..ebcf86e877b0c150df74ee65e105e8f9bdbf13c3 100644
--- a/+misc/+ss/simulationOutput2Quantity.m
+++ b/+misc/+ss/simulationOutput2Quantity.m
@@ -110,16 +110,14 @@ for idx = 1 : numel(signalNames)
z = quantity.Domain("z", thisGrid{1});
end
quantityStruct = misc.setfield(quantityStruct, signalNames{idx}, ...
- quantity.Discrete(thisSignalSorted, ...
- 'name', signalNames{idx}, ...
- "domain", [t, z]));
+ quantity.Discrete(thisSignalSorted, [t, z], ...
+ 'name', signalNames{idx}));
else
quantityStruct = misc.setfield(quantityStruct, signalNames{idx}, ...
quantity.Discrete( ...
- simOutput(:, strcmp(outputNamesUnenumerated, signalNames{idx})), ...
- 'name', signalNames{idx}, ....
- 'domain', t));
+ simOutput(:, strcmp(outputNamesUnenumerated, signalNames{idx})), t, ...
+ 'name', signalNames{idx}));
end
end
end
diff --git a/+misc/Backstepping.m b/+misc/Backstepping.m
index 8cad2dee311ed7aa2df24f3b2842270385560e4c..64d4d3b1092a2d2856e5549adac75ab69922e9ae 100644
--- a/+misc/Backstepping.m
+++ b/+misc/Backstepping.m
@@ -115,7 +115,7 @@ classdef Backstepping < handle & matlab.mixin.Copyable
end
xOriginal = quantity.Discrete(...
ones(101, size(obj.kernel, 1)), ...
- 'grid', {linspace(0, 1, 101)}, 'gridName', {'z'});
+ quantity.Domain('z', linspace(0, 1, 101)));
xTilde = obj.transform(xOriginal);
xOriginalAgain = obj.transformInverse(xTilde);
myError = max(median(abs(xOriginal-xOriginalAgain)));
@@ -238,11 +238,11 @@ classdef Backstepping < handle & matlab.mixin.Copyable
% create quantities as output parameters
K_dzeta = quantity.Discrete(K_dzetaMat, ...
- 'domain', quantityDomain , ...
+ quantityDomain , ...
'name', "d_{zeta}" + thisKernel(1).name);
K_dz = quantity.Discrete(K_dzMat, ...
- 'domain', quantityDomain , ...
+ quantityDomain , ...
'name', "d_z" + thisKernel(1).name);
end % gradient()
@@ -269,11 +269,11 @@ classdef Backstepping < handle & matlab.mixin.Copyable
syms z zeta
if isequal(obj.integralBounds, {0, 'z'})
domainSelector = quantity.Discrete(quantity.Symbolic(zeta<=z, ...
- 'domain', obj.kernel(1).domain));
+ obj.kernel(1).domain));
elseif isequal(obj.integralBounds, {'z', 1})
domainSelector = quantity.Discrete(quantity.Symbolic(zeta>=z, ...
- 'domain', obj.kernel(1).domain));
+ obj.kernel(1).domain));
end
end % get.domainSelector()
diff --git a/+misc/Gain.m b/+misc/Gain.m
index f3f4ca4b167106a223b2d13c54f23f15d7282b42..3c74b0cf88f21967c392b09adba4ee2cdaf9bf63 100644
--- a/+misc/Gain.m
+++ b/+misc/Gain.m
@@ -256,6 +256,20 @@ classdef Gain < handle & matlab.mixin.Copyable
end
end % strrepOutputType()
+ function obj = extendOutputType(obj, position, newText)
+ % extendOutputType adds a pre- or a postfix to obj.outputType
+ assert(any(strcmp(position, {'front', 'back'})));
+ assert(ischar(newText) || isstring(newText), 'prefix must be a char-array or string');
+
+ for it = 1 : numel(obj.outputType)
+ if strcmp(position, 'front')
+ obj.outputType{it} = [newText, obj.outputType{it}];
+ elseif strcmp(position, 'back')
+ obj.outputType{it} = [obj.outputType{it}, newText];
+ end
+ end
+ end % extendOutputType()
+
function obj = strrepInputType(obj, oldText, newText)
% replace strings in type
assert(ischar(oldText) || isstring(oldText), 'oldText must be a char-array or string');
@@ -263,6 +277,18 @@ classdef Gain < handle & matlab.mixin.Copyable
obj.inputType = strrep(obj.inputType, oldText, newText);
end % strrepInputType()
+ function obj = extendInputType(obj, position, newText)
+ % extendInputType adds a pre- or a postfix to obj.inputType
+ assert(any(strcmp(position, {'front', 'back'})));
+ assert(ischar(newText) || isstring(newText), 'prefix must be a char-array or string');
+
+ if strcmp(position, 'front')
+ obj.inputType = [newText, obj.inputType];
+ elseif strcmp(position, 'back')
+ obj.inputType = [obj.inputType, newText];
+ end
+ end % extendInputType()
+
function OutputName = get.OutputName(obj)
OutputName = cell(sum(obj.lengthOutput), 1);
myCounter = 1;
diff --git a/+misc/Gains.m b/+misc/Gains.m
index e323495f329e0eb63e84fb019f8786b86ab86348..19bc477700704a9e7d5198404e1b0718b764dbc5 100644
--- a/+misc/Gains.m
+++ b/+misc/Gains.m
@@ -378,6 +378,13 @@ classdef Gains < handle & matlab.mixin.Copyable
end
end % strrepOutputType()
+ function obj = extendOutputType(obj, position, newText)
+ % extendOutputType adds a pre- or a postfix to all obj.outputType{:}
+ for it = 1 : obj.numTypes
+ obj.gain(it).extendOutputType(position, newText);
+ end
+ end % extendOutputType()
+
function obj = strrepInputType(obj, oldText, newText)
if ~isempty(obj)
% replace strings in type
@@ -387,6 +394,13 @@ classdef Gains < handle & matlab.mixin.Copyable
end
end % strrepInputType()
+ function obj = extendInputType(obj, position, newText)
+ % extendInputType adds a pre- or a postfix to all obj.inputType{:}
+ for it = 1 : obj.numTypes
+ obj.gain(it).extendInputType(position, newText);
+ end
+ end % extendInputType()
+
%% set gain value
function set.value(obj, valueTemp)
for it = 1 : obj.numTypes
diff --git a/+misc/Parser.m b/+misc/Parser.m
index b8d34560ae708cff0607fe5c6f82c3dafb8a7adc..cefb692927c27cfd3aa80cbd1c8a5c6c532d47af 100644
--- a/+misc/Parser.m
+++ b/+misc/Parser.m
@@ -15,8 +15,10 @@ classdef Parser < inputParser
obj.PartialMatching = false;
end % Parser()
- function i = isDefault(obj, name)
- i = any(ismember(obj.UsingDefaults(:), {name}));
+ function result = isDefault(obj, name)
+ % isDefault returns true if for the Parameter specified by NAME is the default value is
+ % used, and returns false otherwise.
+ result = any(ismember(obj.UsingDefaults(:), {name}));
end % isDefault()
function parse2ws(obj, varargin)
diff --git a/+misc/alln.m b/+misc/alln.m
index 485da5bee56fd6de158f6e57966fc8fb4cd5815a..2c1d184564bcaeaed4f15bdd9350280f5a8e80fb 100644
--- a/+misc/alln.m
+++ b/+misc/alln.m
@@ -23,4 +23,6 @@ import misc.alln
V = all(V(:));
+warning("depricated: use all(arg, 'all') instead of misc.alln(arg)");
+
end
\ No newline at end of file
diff --git a/+misc/binomial.m b/+misc/binomial.m
new file mode 100644
index 0000000000000000000000000000000000000000..a382b634772300767bf3ce5d5cb90a487d44d8ae
--- /dev/null
+++ b/+misc/binomial.m
@@ -0,0 +1,16 @@
+function b = binomial(p, k)
+%UNTITLED Summary of this function goes here
+% Detailed explanation goes here
+arguments
+ p (1,1) double;
+ k (1,1) double {mustBeInteger};
+end
+
+if p < 0 && k ~= 0
+ b = (-1)^k * misc.binomial( -p + k -1, k);
+else
+ b = gamma(p+1) / gamma(k+1) / gamma(p-k + 1);
+end
+
+end
+
diff --git a/+misc/functionArguments.m b/+misc/functionArguments.m
new file mode 100644
index 0000000000000000000000000000000000000000..6e7de0bc7a8d746cff61661ef924febddc2c3076
--- /dev/null
+++ b/+misc/functionArguments.m
@@ -0,0 +1,10 @@
+function names = functionArguments(func)
+%FUNCTIONARGUMENTS extracts the arguments of a function handle
+% names = misc.functionArguments(func) returns the arguments of the function handle func as the
+% string array names
+%% Example
+% names = functionArguments(@(z, zeta, t) z*zeta+t)
+
+names = split(string(regexp(func2str(func), "(?<=@\().*(?=\))", "match")), ",");
+end
+
diff --git a/+misc/indexGrid.m b/+misc/indexGrid.m
new file mode 100644
index 0000000000000000000000000000000000000000..2e5d270cb991294d2c3ef79a1924665f77819bdd
--- /dev/null
+++ b/+misc/indexGrid.m
@@ -0,0 +1,18 @@
+function myGrid = indexGrid(mySize)
+% INDEXGRID creates a cell array of grid vectors. For every element of the mySize integer array one
+% cell contain 1 : 1 : mySize(it) is created.
+%% Example
+% myGrid = misc.indexGrid([1, 2, 3]);
+% myGrid{:}
+% ans =
+% 1
+% ans =
+% 1 2
+% ans =
+% 1 2 3
+
+myGrid = cell(numel(mySize), 1);
+for it = 1 : numel(myGrid)
+ myGrid{it} = 1 : 1 : mySize(it);
+end
+end % indexGrid()
\ No newline at end of file
diff --git a/+misc/nearInt.m b/+misc/nearInt.m
new file mode 100644
index 0000000000000000000000000000000000000000..22d4d6a69eee104092e5243b95743a6612a40281
--- /dev/null
+++ b/+misc/nearInt.m
@@ -0,0 +1,7 @@
+function result = nearInt(number, tolerance)
+ if nargin == 1
+ result = numeric.near(number, round(number));
+ else
+ result = numeric.near(number, round(number), tolerance);
+ end
+end
\ No newline at end of file
diff --git a/+misc/solveVolterraIntegroDifferentialIVP.m b/+misc/solveVolterraIntegroDifferentialIVP.m
index 62773ea4a9be85384833c01a7478fc6646bbf8ce..0bec272d5f771b1d03b34a41d08419d5feeca265 100644
--- a/+misc/solveVolterraIntegroDifferentialIVP.m
+++ b/+misc/solveVolterraIntegroDifferentialIVP.m
@@ -65,7 +65,7 @@ if NameValueArgs.D.abs.MAX() >= 1e3*eps
c5a = permute(numeric.trapz_fast_nDim(z.grid(1:zIdx), integrand), [2, 3, 1]);
end % for zIdx = 2 : z.n
- N = quantity.Discrete(permute(Ndouble, [3, 1, 2]), 'domain', z, 'name', "N");
+ N = quantity.Discrete(permute(Ndouble, [3, 1, 2]), z, 'name', "N");
for it = 1 : 2
% improve by result by successive approximation for which
% Cnew = C(z) + int_0^z D(z, zeta) N(zeta) dzeta
@@ -87,7 +87,7 @@ else
z.grid.', ...
ic(:));
- N = quantity.Discrete(reshape(myN, [z.n, n, m]), 'domain', z, 'name', "N");
+ N = quantity.Discrete(reshape(myN, [z.n, n, m]), z, 'name', "N");
end
diff --git a/+mustBe/scalarOrEmpty.m b/+mustBe/scalarOrEmpty.m
new file mode 100644
index 0000000000000000000000000000000000000000..1000287d959d5bcfdaf13f2e5cdd5e8b33f5e24b
--- /dev/null
+++ b/+mustBe/scalarOrEmpty.m
@@ -0,0 +1,7 @@
+function scalarOrEmpty(A)
+%mustBe.scalarOrEmpty Validates if the argument is empty or a scalar value
+
+ if numel(A) ~= 1 && ~isempty(A)
+ error('Value assigned to Data property must be a scalar value or empty.');
+ end
+end
\ No newline at end of file
diff --git a/+numeric/femMatrices.m b/+numeric/femMatrices.m
index d84d0e9ef1bb388e3fe9e35f6f15287f2c96a499..abb7e78f278e6e6394405d53ce9c4b0ec33dd4c8 100644
--- a/+numeric/femMatrices.m
+++ b/+numeric/femMatrices.m
@@ -1,73 +1,72 @@
-function [M, K, L, P, PHI] = femMatrices(varargin)
+function [M, K, L, P, PHI, dphi] = femMatrices(spatialDomain, optArgs)
%FEMMATRICES computes the discretization matrices obtained by FE-Method
-% [M, K, L, P, PHI] = femMatrices(varargin) computes the matrices required
-% for the approximation of a system using FEM. For instance, the rhs of a pde
-% dz (alpha(z) dz x(z)) + beta(z) dz x(z) + gamma(z)
-% is approximated as K x(z)
-% The matrices that will be returned as described in the following:
+% [M, K, L, P, PHI, DPHI] = femMatrices(SPATIALDOMAIN, OPTARGS) computes the matrices required for
+% the approximation of a system using FEM. For instance, the rhs of a pde
+% PDE[ x(z) ] = dz (alpha(z) dz x(z)) + beta(z) dz x(z) + gamma(z) x(z) + b(z) u
+% y = c(z) x(z)
+% is approximated as K x(z) The matrices that will be returned are:
% The mass-matrix M is computed using
% M = int phi(z) phi(z)
% The stiffness matrix K is computed with
-% K0 = int phi(z) gamma(z) phi(z)
-% K1 = int phi(z) beta(z) (dz phi(z))
-% K2 = int - (dz phi(z)) alpha(z) (dz phi(z))
-% K = K0 + K1 + K2
+% K0 = int phi(z) gamma(z) phi(z) K1 = int phi(z) beta(z) (dz phi(z)) K2 = int - (dz phi(z))
+% alpha(z) (dz phi(z)) K = K0 + K1 - K2
% The input matrix L:
% L = int( phi(z) * b(z) )
% The output matrix P:
% P = int( c(z) * phi'(z) )
%
-% PHI is the trial function used for the discretization.
+% PHI is the trial function used for the discretization and DPHI is the derivative of the trial
+% function.
%
-% The input arguments can be given as name-value-pairs. The parameters
-% are described in the following, the default values are written in (.).
-% 'grid' (linspace(0, 1, 11)): grid of resulting finite elements
-% 'Ne'(111): the implementation uses numerical integration for each element.
-% The parameter 'Ne' defines the grid for this integration.
-% 'alpha' : the system parameter alpha, should be a quantity.Discrete.
-% 'beta' : the system parameter beta, should be a quantity.Discrete
-% 'gamma' : the system parameter gamma, should be a quantity.Discrete
-% 'b' : the input vector, should be a quantity.Discrete
-% 'c' : the output vector, should be a quantity.Discrete
-
-myParser = misc.Parser();
-myParser.addOptional('grid', linspace(0, 1, 11));
-myParser.addOptional('Ne', 111);
-myParser.addOptional('alpha', quantity.Discrete.empty(1,0));
-myParser.addOptional('beta', quantity.Discrete.empty(1,0));
-myParser.addOptional('gamma', quantity.Discrete.empty(1,0));
-myParser.addOptional('b', quantity.Discrete.empty(1,0));
-myParser.addOptional('c', quantity.Discrete.empty(0,1));
-myParser.addOptional('silent', false);
-myParser.parse(varargin{:});
-
-% set the system properties
-z = myParser.Results.grid(:);
-N = numel(z);
-z0 = z(1);
-z1 = z(end);
+% The argument SPATIALDOMAIN is mandatory and must be a quantity.Domain object. The further system
+% parameters as
+% "alpha", "beta", "gamma", "b", "c"
+% must be quantity.Discrete objects. With the parameter
+% "Ne"
+% as positive integer number one can specify the number of discretization points for one element.
+% The parameter
+% "silent"
+% is a logical flag to supress the progressbar.
+%
+
+arguments
+ spatialDomain (1,1) quantity.Domain;
+ optArgs.Ne (1,1) double {mustBeInteger} = 111;
+ optArgs.alpha quantity.Discrete {mustBe.scalarOrEmpty} = quantity.Discrete.empty(1,0);
+ optArgs.beta quantity.Discrete {mustBe.scalarOrEmpty} = quantity.Discrete.empty(1,0);
+ optArgs.gamma quantity.Discrete {mustBe.scalarOrEmpty} = quantity.Discrete.empty(1,0);
+ optArgs.b quantity.Discrete = quantity.Discrete.empty(1,0);
+ optArgs.c quantity.Discrete = quantity.Discrete.empty(1,0)';
+ optArgs.silent logical = false;
+end
+
+% set the system properties to local variables, for better readability and (maybe) faster access
+z = spatialDomain.grid;
+N = spatialDomain.n;
+z0 = spatialDomain.lower;
+z1 = spatialDomain.upper;
+
pbar = misc.ProgressBar('name', 'FEM', 'terminalValue', N-1, 'steps', N-1, ...
- 'initialStart', true, 'silent', myParser.Results.silent);
+ 'initialStart', true, 'silent', optArgs.silent);
-alf = myParser.Results.alpha;
-bta = myParser.Results.beta;
-gma = myParser.Results.gamma;
-b = myParser.Results.b;
-c = myParser.Results.c;
+alf = optArgs.alpha;
+bta = optArgs.beta;
+gma = optArgs.gamma;
+b = optArgs.b;
+c = optArgs.c;
p = size(b, 2); % number of the input values
m = size(c, 1); % number of the output values
% define the properties of one finite element
-element.N = myParser.Results.Ne;
-element.delta = (z1 - z0) / (N - 1);
-element.z = linspace(0, element.delta, element.N)';
+element = quantity.EquidistantDomain("z", 0, (z1 - z0) / (N - 1), ...
+ 'stepNumber', optArgs.Ne);
%% verify that grid is homogenious
-% Yet, this function only supports homogenious grids for most output parameters.
-% Only, the mass matrix M is calculated correct for nonhomogenious grids. The
-% following line ensures, that if the grid is not homogenious, only M is returned.
-isGridHomogenious = all(diff(z) - element.delta < 10*eps);
+% Yet, this function only supports homogenious grids for most output parameters. Only, the mass
+% matrix M is calculated correct for nonhomogenious grids. The following line ensures, that if the
+% grid is not homogenious, only M is returned.
+isGridHomogenious = all(diff(z) - element.upper < 10*eps);
if ~isGridHomogenious
assert(nargout <= 1, 'non-homogenious grids only supported for mass-matrix');
end
@@ -81,22 +80,22 @@ end
% | ------------ |
% \ deltaZ_k /
-PHI = quantity.Discrete({(element.delta - element.z) / element.delta; ...
- (element.z / element.delta) }, ...
- 'size', [2, 1], 'gridName', 'z', 'grid', element.z, 'name', 'phi');
+PHI = quantity.Discrete({(element.upper - element.grid) / element.upper; ...
+ (element.grid / element.upper) }, ...
+ element, 'size', [2, 1], 'name', 'phi');
%% compute the mass matrix
M0 = 1 / 6 * [2, 1; 1, 2];
-element.diff = diff(z);
+zDiff = diff(z);
%% compute the stiffness matrix
-% initialize the trial functions for the computation for element matrices.
-% for further comments, see the description in the loop
-pp = reshape(double(PHI * PHI'), element.N, 4);
+% initialize the trial functions for the computation for element matrices. for further comments, see
+% the description in the loop
+pp = reshape(double(PHI * PHI'), element.n, 4);
phi = double(PHI);
-dphi = [ -1 / element.delta; ...
- 1 / element.delta ];
+dphi = [ -1 / element.upper; ...
+ 1 / element.upper ];
dphidphi = reshape(dphi * dphi', 1, 4);
% initialize the element matrices
@@ -122,51 +121,48 @@ flags.c = ~isempty(c);
for k = 1:N-1
pbar.raise(k);
idxk=[k, k+1];
- zk = linspace(z(k), z(k+1), element.N)';
+ zk = linspace(z(k), z(k+1), element.n)';
if flags.gamma
% compute the elements for K0 = int phi(z) * gamma(z) * phi'(z)
- % the matrix phi(z) * phi'(z) is the same for each element. Thus,
- % this matrix is computed a priori. To simplify the multiplication
- % with the scalar value gamma(z), the matrix phi*phi' is vectorized
- % and the result of phi(z)phi'(z) alpha(z) is reshaped as a matrix.
+ % the matrix phi(z) * phi'(z) is the same for each element. Thus, this matrix is computed
+ % a priori. To simplify the multiplication with the scalar value gamma(z), the matrix
+ % phi*phi' is vectorized and the result of phi(z)phi'(z) alpha(z) is reshaped as a matrix.
k0 = misc.multArray(pp, gma.on(zk), 3, 2, 1);
K0 = reshape( numeric.trapz_fast(zk, k0'), 2, 2);
end
if flags.beta
% compute the elements for K1 = int phi(z) * beta(z) * dz phi'(z))
- % as the values dz phi(z) are constants, the multiplication can be.
- % simplified.
+ % as the values dz phi(z) are constants, the multiplication can be. simplified.
k1 = misc.multArray(phi, bta.on(zk) * dphi', 3, 3, 1);
- K1 = reshape( numeric.trapz_fast(zk, reshape(k1, element.N, 4, 1)' ), 2, 2);
+ K1 = reshape( numeric.trapz_fast(zk, reshape(k1, element.n, 4, 1)' ), 2, 2);
end
if flags.alpha
% compute the elements for K2 = int dz phi(z) * alpha(z) * dz phi'(z)
- % as the matrix dzphi*dzphi' is constnat the multiplication
- % is simplified by vectorization of dzphi*dzphi' and a reshape to
- % obtain the correct result.
+ % as the matrix dzphi*dzphi' is constnat the multiplication is simplified by vectorization
+ % of dzphi*dzphi' and a reshape to obtain the correct result.
k2 = alf.on(zk) * dphidphi;
K2 = reshape( numeric.trapz_fast(zk, k2'), 2, 2);
end
if flags.b
% compute the input matrix L = int phi(z) * b(z)
- l = reshape( misc.multArray(phi, b.on(zk), 3, 2, 1), element.N, 2*p);
+ l = reshape( misc.multArray(phi, b.on(zk), 3, 2, 1), element.n, 2*p);
L0 = reshape( numeric.trapz_fast(zk, l'), 2, p);
end
if flags.c
% compute the output matrix P = int c(z) * phi(z)'
- c0 = reshape( misc.multArray(c.on(zk), phi, 3, 3, 1), element.N, 2*m);
+ c0 = reshape( misc.multArray(c.on(zk), phi, 3, 3, 1), element.n, 2*m);
P0 = reshape( numeric.trapz_fast(zk, c0'), m, 2);
end
% assemble the matrices with the elements:
- M(idxk, idxk) = M(idxk, idxk) + M0 * element.diff(k);
+ M(idxk, idxk) = M(idxk, idxk) + M0 * zDiff(k);
- K(idxk, idxk) = K(idxk, idxk) + K0 + K1 + K2;
+ K(idxk, idxk) = K(idxk, idxk) + K0 + K1 - K2;
L(idxk, :) = L(idxk, :) + L0;
P(:, idxk) = P(:, idxk) + P0;
end
diff --git a/+quantity/BasicVariable.m b/+quantity/BasicVariable.m
index 8aa5af204359e3fd08bbfb280bb66e5cfec7ab81..d0177ea0ae018d09caddba02c6d921567726f6b8 100644
--- a/+quantity/BasicVariable.m
+++ b/+quantity/BasicVariable.m
@@ -33,8 +33,8 @@ end
methods
%% CONSTRUCTOR
-function obj = BasicVariable(valueContinuous, myDomain, varargin)
-
+function obj = BasicVariable(valueContinuous, varargin)
+ warning("DEPRICATED!: Use signals.BasicVariable.")
parentVarargin = {};
if nargin > 0 && ~isempty(varargin)
@@ -111,6 +111,7 @@ function obj = BasicVariable(valueContinuous, myDomain, varargin)
if size(obj_j.derivatives, 1) < k+1 || isempty(obj_j.derivatives(k+1,:))
% create a new object for the derivative of obj:
Dij = obj_j.copy();
+ Dij.offset = 0;
Dij.name = ['d/dt (' num2str(k) ') '];
Dij.valueDiscrete = [];
Dij.diffShift = Dij.diffShift + k;
@@ -168,7 +169,6 @@ function obj = BasicVariable(valueContinuous, myDomain, varargin)
bi.K = K;
b.derivatives(i) = bi;
end
-
end
end
@@ -211,4 +211,4 @@ methods (Static)
end
end
-end
\ No newline at end of file
+end
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 9969feb822978626b405d2a5f075d4105d64d28b..17f27dfb9877c4f758c87bf7cea4466cacc143b2 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -6,10 +6,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
valueDiscrete double;
end
- properties (Hidden, Access = protected, Dependent)
- doNotCopy;
- end
-
properties
% ID of the figure handle in which the handle is plotted
figureID double = 1;
@@ -29,9 +25,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% Grid for the evaluation of the continuous quantity. For the
% example with the function f(x,t), the grid would be
% {[<spatial domain>], [<temporal domain>]}
- % whereas <spatial domain> is the discret description of the
- % spatial domain and <temporal domain> the discrete description of
- % the temporal domain.
+ % whereas <spatial domain> is the discrete description of the
+ % spatial domain x and <temporal domain> the discrete description of
+ % the temporal domain t.
grid; % in set.grid it is ensured that, grid is a (1,:)-cell-array
end
@@ -48,9 +44,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% size(valueOriginal) == size(obj) and
% size(valueOriginal{it}) == gridSize
% Example: valueOrigin = { f(Z, T), g(Z, T) } is a cell array
- % wich contains the functions f(z,t) and g(z,t) evaluated on
- % the discrete domain (Z x T). Then, the name-value-pair
- % parameter 'domain' must be set with quantity.Domain
+ % which contains the functions f(z,t) and g(z,t) evaluated on
+ % the discrete domain (Z x T). Then, myDOmain must be quantity.Domain
% objects, according to the domains Z and T.
% OR
% 2) a double-array with
@@ -99,9 +94,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% the case if all values are empty. This is required for
% the initialization of quantity.Function and
% quantity.Symbolic objects
- assert( misc.alln( cellfun(@isempty, valueOriginal ) ) || ...
- numGridElements(myDomain) == numel(valueOriginal{1}), ...
- 'grids do not fit to valueOriginal');
+ assert(isempty(myDomain) ... % constant case
+ || all( cellfun(@isempty, valueOriginal ), 'all' ) ... % empty case
+ || isequal([myDomain.n], size(valueOriginal{1}, 1 : max(1, numel(myDomain)))) ... % usual case
+ || (isrow(valueOriginal{1}) && ... % row-vector case (including next line)
+ isequal([1, myDomain.n], size(valueOriginal{1}, 1 : max(1, numel(myDomain)+1)))), ...
+ 'grids do not fit to valueOriginal');
% allow initialization of empty objects
valueOriginalSize = size(valueOriginal);
@@ -163,9 +161,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
itIs = isempty(obj(1).domain);
end % isConstant()
- function doNotCopy = get.doNotCopy(obj)
- doNotCopy = obj.doNotCopyPropertiesName();
- end
function valueDiscrete = get.valueDiscrete(obj)
% check if the value discrete for this object
% has already been computed.
@@ -175,7 +170,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
valueDiscrete = obj.valueDiscrete;
end
-
+
%-------------------
% --- converters ---
%-------------------
@@ -189,7 +184,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
headers{i+1} = obj(i).name + "" + num2str(i);
end
exportData = export.dd(...
- 'M', [obj.grid{:}, obj.valueDiscrete], ...
+ 'M', [obj(1).grid{:}, obj.valueDiscrete], ...
'header', headers, varargin{:});
elseif obj.nargin == 2
error('Not yet implemented')
@@ -235,6 +230,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% end
[obj.name] = deal(newName);
end % setName()
+
+
+
end
methods (Access = public)
@@ -248,7 +246,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for i = 1:numel(obj)
- data = obj(i).valueDiscrete();
+ data = obj(i).valueDiscrete;
zeros = find( abs(data) <= optArg.tol);
upCrossing = find( data(1:end-1) <= 0 & data(2:end) > 0);
@@ -275,8 +273,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
data = {[ obj.valueDiscrete ], obj(1).domain.name, ...
obj(1).domain.grid, obj(1).name};
h = misc.hash(data);
-
- end
+ end % hash()
function d = compositionDomain(obj, domainName)
@@ -293,7 +290,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
function obj_hat = compose(obj, g, optionalArgs)
% COMPOSE compose two functions
- % OBJ_hat = compose(obj, G, varargin) composes the function
+ % OBJ_hat = compose(obj, G, varargin) composes the function f
% defined by OBJ with the function given by G. In particular,
% f_hat(z,t) = f( z, g(z,t) )
% if f(t) = obj, g is G and f_hat is OBJ_hat.
@@ -398,7 +395,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
'name', obj.name + "°" + g.name, ...
'size', size(obj));
- end
+ end % compose()
function value = on(obj, myDomain, gridNames)
% ON evaluation of the quantity on a certain domain.
@@ -420,7 +417,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
if nargin == 1
% case 0: no domain was specified, hence the value is requested
% on the default grid defined by obj(1).domain.
- value = obj.obj2value(obj(1).domain);
+ value = reshape(cat(numel(obj(1).domain)+1, obj(:).valueDiscrete), ...
+ [obj(1).domain.gridLength(), size(obj)]);
elseif nargin == 2 && (iscell(myDomain) || isnumeric(myDomain))
% case 1: a domain is specified by myDomain as agrid
@@ -439,8 +437,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for k = 1:length(newGrid)
myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
end
- value = reshape(obj.obj2value(myDomain), ...
- [myDomain.gridLength, size(obj)]);
+ value = obj.obj2value(myDomain);
else
% Since in the remaining cases the order of the domains is not
% neccessarily equal to the order in obj(1).domain, this is
@@ -463,7 +460,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
'The cell entries for a new grid have to be vectors')
newGrid = myDomain;
- myDomain = quantity.Domain.empty();
+ clear('myDomain');
+ myDomain(1:length(newGrid)) = quantity.Domain();
for k = 1:length(newGrid)
myDomain(k) = quantity.Domain(gridNames{k}, newGrid{k});
end
@@ -480,6 +478,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% compute the permutation index, in order to bring the
% new domain in the same order as the original one.
gridPermuteIdx = obj(1).domain.getPermutationIdx(myDomain);
+
+ assert(any(gridPermuteIdx ~= 0), "grid could not be found.")
end
% get the valueDiscrete data for this object. Apply the
% permuted myDomain. Then the obj2value will be evaluated
@@ -487,8 +487,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% the new order will be done in the next step.
originalOrderedDomain(gridPermuteIdx) = myDomain;
value = obj.obj2value(originalOrderedDomain);
- value = permute(reshape(value, [originalOrderedDomain.gridLength, size(obj)]), ...
- [gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
+ value = permute(value, [gridPermuteIdx, numel(gridPermuteIdx)+(1:ndims(obj))]);
end
end % if isempty(obj)
end % on()
@@ -497,17 +496,17 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% get the interpolant of the obj;
if isempty(obj)
value = zeros(size(obj));
- indexGrid = arrayfun(@(s)linspace(1,s,s), size(obj), 'UniformOutput', false);
+ indexGrid = misc.indexGrid(size(obj));
interpolant = numeric.interpolant(...
[indexGrid{:}], value);
else
myGrid = obj(1).grid;
value = obj.obj2value();
- indexGrid = arrayfun(@(s)linspace(1,s,s), size(obj), 'UniformOutput', false);
+ indexGrid = misc.indexGrid(size(obj));
interpolant = numeric.interpolant(...
[myGrid, indexGrid{:}], value);
end
- end
+ end % interpolant()
function assertSameGrid(a, varargin)
@@ -584,20 +583,20 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
end
- function obj = sort(obj, varargin)
+ function newObj = sort(obj, varargin)
%SORT sorts the grid of the object in a desired order
% obj = sortGrid(obj) sorts the grid in alphabetical order.
% obj = sort(obj, 'descend') sorts the grid in descending
% alphabetical order.
-
+ newObj = obj.copy();
% only sort the grids if there is something to sort
- if ~isempty(obj) && obj(1).nargin > 1
+ if ~isempty(newObj) && newObj(1).nargin > 1
- [sortedDomain, I] = obj(1).domain.sort(varargin{:});
- [obj.domain] = deal(sortedDomain);
+ [sortedDomain, I] = newObj(1).domain.sort(varargin{:});
+ [newObj.domain] = deal(sortedDomain);
- for k = 1:numel(obj)
- obj(k).valueDiscrete = permute(obj(k).valueDiscrete, I);
+ for k = 1:numel(newObj)
+ newObj(k).valueDiscrete = permute(obj(k).valueDiscrete, I);
end
end
end% sort()
@@ -849,7 +848,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
inverse = quantity.Discrete(...
repmat(obj(1).grid{obj(1).domain.gridIndex(gridName)}(:), [1, size(obj)]), ...
quantity.Domain([obj(1).name], obj.on()), ...
- 'size', size(obj), ...
'name', gridName);
end % invert()
@@ -970,6 +968,13 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
assert(numel(values) == numel(gridName2Replace), ...
'gridName2Replace and values must be of same size');
+ % set the newDomain once. If obj(1).domain is a quantity.Equidistant domain, it can
+ % not be mixed with other quantity.Domains in an array. Hence, it must be casted to
+ % a quantity.Domain. The following strange form of concatenation an empty Domain
+ % with the required domain, ensures that the result is an array of quantity.Domain
+ % objects.
+ newDomain = [quantity.Domain.empty, obj(1).domain];
+
% here substitution starts:
% The first (gridName2Replace{1}, values{1})-pair is
% replaced. If there are more cell-elements in those inputs
@@ -998,7 +1003,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% Do not use "return", since, later subs might need to be
% called recursively!
newValue = obj.on();
- newDomain = obj(1).domain;
elseif any(strcmp(values{1}, obj(1).gridName))
% if for a quantity f(z, zeta) this method is
% called with subs(f, zeta, z), then g(z) = f(z, z)
@@ -1021,7 +1025,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
newDomainForOn(all(1:1:numel(newDomainForOn) ~= domainIndices(:)))];
else
% this is the default case. just grid name is changed.
- newDomain = obj(1).domain;
+ %newDomain = obj(1).domain;
newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
quantity.Domain(values{1}, ...
obj(1).domain(obj(1).domain.gridIndex(gridName2Replace{1})).grid);
@@ -1032,7 +1036,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% if values{1} is a scalar, then obj is evaluated and
% the resulting quantity looses that spatial grid and
% gridName
- newDomain = obj(1).domain;
newDomain = newDomain(~strcmp(gridName2Replace{1}, [newDomain.name]));
% newGrid is the similar to the original grid, but the
% grid of gridName2Replace is removed.
@@ -1047,14 +1050,16 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
elseif isnumeric(values{1}) && numel(values{1}) > 1
% if values{1} is a double vector, then the grid is
% replaced.
- newDomain = obj(1).domain;
newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
quantity.Domain(gridName2Replace{1}, values{1});
newValue = obj.on(newDomain);
else
error('value must specify a gridName or a gridPoint');
- end
+ end % if ischar(values{1}) || isstring(values{1})
if isempty(newDomain)
+ % TODO@Jakob: Is it possible that this case occurrs?
+ % If I see it right, the only case where the newDomain is not set is the case
+ % where a error is thrown.
solution = newValue;
else
solution = quantity.Discrete(newValue, newDomain, ...
@@ -1075,7 +1080,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% at() evaluates the object at one point and returns it as array
% with the same size as size(obj).
value = reshape(obj.on(point), size(obj));
- end
+ end % at()
function value = atIndex(obj, varargin)
% ATINDEX returns the valueDiscrete at the requested index.
@@ -1118,7 +1123,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
value = reshape([value{:}], [outputSize, size(obj)]);
end
- end
+ end % atIndex()
function n = nargin(obj)
% FIXME: check if all funtions in this object have the same
@@ -1323,17 +1328,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
function s = struct(obj, varargin)
if (nargin == 1) && isa(obj, "quantity.Discrete")
- properties = fieldnames(obj);
+ tempProperties = fieldnames(obj);
si = num2cell( size(obj) );
s(si{:}) = struct();
for l = 1:numel(obj)
-
- doNotCopyProperties = obj(l).doNotCopy;
-
- for k = 1:length(properties)
- if ~any(strcmp(doNotCopyProperties, properties{k}))
- s(l).(properties{k}) = obj(1).(properties{k});
- end
+ for k = 1:length(tempProperties)
+ s(l).(tempProperties{k}) = obj(1).(tempProperties{k});
end
end
else
@@ -1392,8 +1392,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
gridNameNew = [gridNew.name];
gridNew = {gridNew.grid};
else
- gridNameNew = misc.ensureString(gridNameNew);
+ gridNameNew = misc.ensureString(gridNameNew);
gridNew = misc.ensureIsCell(gridNew);
+ for it = 1:numel(gridNew)
+ assert( isnumeric( [gridNew{it}] ), "The gridNew parameter must be a cell array of numeric arrays." )
+ end
end
if obj(1).isConstant
@@ -1424,7 +1427,55 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for it = 1 : numel(obj)
newObj(it).valueDiscrete = obj(it).on(newDomain);
end
- end
+ end % changeGrid()
+
+
+ function newObj = replaceGrid(obj, myNewDomain, optArgs)
+ % CHANGEGRID change the grid of the quantity.
+ % newObj = REPLACEGRID(obj, MYNEWDOMAIN, "gridName", NEWGRIDNAME)
+ % replace the grid of the obj quantity. The order of grid and
+ % gridName in the obj properties remains unchanged, only the
+ % grid points are exchanged.
+ %
+ % TODO:
+ % newObj = REPLACEGRID(obj, domain) changes the domain of the
+ % object specified by the name of DOMAIN into DOMAIN.
+ %
+ % newObj = REPLACEGRID(obj, domain, 'gridName', NEWNAME) changes the domain of the
+ % object specified by NEWNAME into DOMAIN.
+
+ arguments
+ obj
+ myNewDomain quantity.Domain
+ optArgs.gridName = [myNewDomain.name];
+ end
+
+ if isempty(obj)
+ newObj = obj.copy();
+ return;
+ end
+
+ assert( intersect(obj(1).gridName, optArgs.gridName) == optArgs.gridName );
+
+ if obj(1).isConstant
+ error("Not yet implemented")
+ else
+ gridIndexNew = obj(1).domain.gridIndex(optArgs.gridName);
+ % initialization of the newDomain array as quantity.Domain
+ % array. This is required in order to handle also
+ % quantity.EquidistantDomains:
+ newDomain(1:obj(1).nargin) = quantity.Domain();
+ newDomain(:) = obj(1).domain;
+
+ for it = 1 : length(gridIndexNew)
+ newDomain(gridIndexNew(it)) = ...
+ quantity.Domain(myNewDomain(it).name, myNewDomain(it).grid);
+ end
+ end
+
+ newObj = obj.copy();
+ [newObj.domain] = deal(newDomain);
+ end % replaceGrid
end
%% math
@@ -1556,26 +1607,50 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
error('sqrtm() is only implemented for quadratic matrices');
end
end % sqrtm()
+
function P = power(obj, p)
- % a.^p implemented by multiplication
- % #todo unittest required
- P = quantity.Discrete.zeros(size(obj), obj(1).domain);
- for k = 1:numel(obj)
- P(k) = obj(k)^p;
- end
- end
+ % a.^p elementwise power
+ P = quantity.Discrete(obj.on().^(p), obj(1).domain, ...
+ "name", obj(1).name + ".^{" + num2str(p) + "}");
+ end % power()
+
function P = mpower(a, p)
- % a^p implemented by multiplication
- assert(p==floor(p) && p > 0);
- P = a;
- for k = 1:(p-1)
- P = P * a;
+ % Matrix power a^p is matrix or scalar a to the power p.
+ if p == 0
+ P = setName(eye(size(a)) + 0*a, "I");
+ elseif p < 0
+ if numel(a) > 1
+ warning("mpower(a, p) implements a^p. " ...
+ + "For matrices a with negative exponent p, inv(a^(-p)) is returned. " ...
+ + "This represents a division from left, " ...
+ + "maybe division from right is needed in your case!")
+ end % this warning is not important in the scalar case.
+ P = inv(mpower(a, -p));
+ else % p > 0
+ assert(p==floor(p), "power p must be integer");
+ P = a;
+ for k = 1:(p-1)
+ P = P * a;
+ end
end
end % mpower()
+
function s = num2str(obj)
s = obj.name;
end % num2str()
+ function P = times(a, b)
+
+ assert( all( size(a) == size(b), 'all' ), 'the terms a and b must have the same size');
+
+ for i = 1:numel(a)
+ P(i) = a(i) * b(i);
+ end
+
+ P = reshape(P, size(a));
+
+ end
+
function P = mtimes(a, b)
% TODO rewrite the selection of the special cases! the
% if-then-cosntruct is pretty ugly!
@@ -1756,10 +1831,62 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
function y = exp(obj)
% exp() is the exponential function using obj as the exponent.
y = quantity.Discrete(exp(obj.on()), obj(1).domain, ...
- 'name', "exp(" + obj(1).name + ")", ...
- 'size', size(obj));
+ 'name', "exp(" + obj(1).name + ")");
end % exp()
+ function [V, D, W] = eig(obj)
+ % eig Eigenvalues and eigenvectors pointwise over the domain
+ % E = eig(A) produces a column vector E containing the eigenvalues of
+ % a square matrix A.
+ %
+ % [V,D] = eig(A) produces a diagonal matrix D of eigenvalues and
+ % a full matrix V whose columns are the corresponding eigenvectors
+ % so that A*V = V*D.
+ %
+ % [V,D,W] = eig(A) also produces a full matrix W whose columns are the
+ % corresponding left eigenvectors so that W'*A = D*W'.
+ assert(numel(obj(1).domain)==1, "quantity.Discrete/eig is only implemented for one domain");
+ assert((ndims(obj) <= 2) && (size(obj, 1) == size(obj, 2)), ...
+ "quantity.Discrete/eig is only implemented for quadratic quantities");
+
+ objDisc = permute(obj.on(), [2, 3, 1]);
+ if nargout == 1
+ % E = eig(A)
+ VDisc = zeros([obj(1).domain.n, size(obj, 1)]);
+ for it = 1 : obj(1).domain.n
+ VDisc(it, :) = eig(objDisc(:, :, it));
+ end % for it = 1 : obj(1).domain.n
+ else
+ % [V,D] = eig(A) or [V,D,W] = eig(A)
+ VDisc = zeros([obj(1).domain.n, size(obj)]);
+ DDisc = zeros([obj(1).domain.n, size(obj)]);
+ WDisc = zeros([obj(1).domain.n, size(obj)]);
+ for it = 1 : obj(1).domain.n
+ [VDisc(it, :, :), DDisc(it, :, :), WDisc(it, :, :)] = ...
+ eig(objDisc(:, :, it));
+ end % for it = 1 : obj(1).domain.n
+ D = quantity.Discrete(DDisc, obj(1).domain, "name", "D");
+ W = quantity.Discrete(WDisc, obj(1).domain, "name", "W");
+ end
+ V = quantity.Discrete(VDisc, obj(1).domain, "name", "V");
+ end % eig()
+
+ function y = log(obj)
+ % log Natural logarithm.
+ % log(X) is the natural logarithm of the elements of X.
+ % Complex results are produced if X is not positive.
+ y = quantity.Discrete(log(obj.on()), obj(1).domain, ...
+ 'name', "log(" + obj(1).name + ")");
+ end % log()
+
+ function y = log10(obj)
+ % log10 Common (base 10) logarithm.
+ % log10(X) is the base 10 logarithm of the elements of X.
+ % Complex results are produced if X is not positive.
+ y = quantity.Discrete(log10(obj.on()), obj(1).domain, ...
+ 'name', "log10(" + obj(1).name + ")");
+ end % log10()
+
function xNorm = l2norm(obj, integralGridName, optArg)
% calculates the l2 norm, for instance
% xNorm = sqrt(int_0^1 x.' * x dz).
@@ -1791,6 +1918,24 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
end % l2norm()
+ function xNorm = norm(obj, p)
+ % norm implements the vector norm, similar to builtin matlab function norm
+ % norm(V,P) returns the p-norm of V defined as SUM(ABS(V).^P)^(1/P).
+ %
+ % norm(X) is the same as norm(X,2). (euclidian norm)
+
+ arguments
+ obj;
+ p (1, 1) double = 2;
+ end
+
+ if p == 2
+ xNorm = sqrt(sum(abs(obj).^p));
+ else
+ xNorm = sum(abs(obj).^p).^(1/p);
+ end
+ end % norm(obj, p)
+
function xNorm = quadraticNorm(obj, varargin)
% calculates the quadratic norm, i.e,
% xNorm = sqrt(x.' * x).
@@ -1856,7 +2001,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
p = a(1).gridSize();
q = p(2);
p = p(1);
- A = a.on();
+ obj = a.on();
B = b.on();
% dimensions
@@ -1869,7 +2014,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for k = 1 : n % rows of P
for l = 1 : m % columns of P
for r = 1 : o % rows of B or columns of A
- P(:, :, k, l) = P( :, :, k, l ) + A( :, :, k, r) .* B( :, :, r, l );
+ P(:, :, k, l) = P( :, :, k, l ) + obj( :, :, k, r) .* B( :, :, r, l );
end
end
end
@@ -1960,7 +2105,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
q = p(2);
p = p(1);
- A = a.on();
+ obj = a.on();
B = repmat(b.on(), 1, 1, 1, q);
B = permute(B, [1, 4, 2, 3]);
@@ -1974,7 +2119,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
for k = 1 : n % rows of P
for l = 1 : m % columns of P
for r = 1 : o % rows of B or columns of A
- P(:, :, k, l) = P( :, :, k, l ) + A( :, :, k, r) .* B( :, :, r, l );
+ P(:, :, k, l) = P( :, :, k, l ) + obj( :, :, k, r) .* B( :, :, r, l );
end
end
end
@@ -2187,7 +2332,11 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
function [P, supremum] = relativeErrorSupremum(A, B)
-
+ %% RELATIVEERRRORSUPREMUM compute the relative error of the supremum
+ % [P, SUPREMUM] = relativeErrorSupremum(A, B) computes the supremum of the absolute
+ % error of A and B and returns the relative error of A - B with respect to this
+ % supremum. The relative error is returned as P and the supremum is returned as
+ % SUPREMUM.
assert(numel(A) == numel(B), 'Not implemented')
P = A.copy();
@@ -2356,12 +2505,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
% obj2value(obj, myDomain) returns the valueDiscrete in the
% form size(value) = [myDomain.gridLength objSize]
- v = zeros([numel(obj(1).valueDiscrete), numel(obj)]);
- for k = 1:numel(obj)
- v(:,k) = obj(k).valueDiscrete(:);
- end
-
- value = reshape(v, [obj(1).domain.gridLength(), size(obj)]);
+ value = reshape(cat(numel(obj(1).domain)+1, obj(:).valueDiscrete), ...
+ [obj(1).domain.gridLength(), size(obj)]);
if nargin >= 2 && ~isequal( myDomain, obj(1).domain )
% if a new domain is specified for the evaluation of
@@ -2369,10 +2514,10 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
if obj.isConstant
% ... duplicate the constant value on the desired
% domain
- value = repmat(v, [cellfun(@(v) numel(v), {myDomain.grid}), ones(1, length(size(obj)))]);
+ value = repmat(value(:).', [myDomain.n, ones(1, ndims(obj))]);
else
%... do an interpolation based on the old data.
- indexGrid = arrayfun(@(s) 1:s, size(obj), 'UniformOutput', false);
+ indexGrid = misc.indexGrid(size(obj));
tempInterpolant = numeric.interpolant(...
[obj(1).grid, indexGrid{:}], value);
tempGrid = {myDomain.grid};
@@ -2380,7 +2525,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
end
end
- end
+ end % obj2value()
function result = diag2vec(obj)
% This method creates a vector of quantities by selecting the
@@ -2477,35 +2622,8 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
q = reshape(mat2cell(value, s{:}), valueSize);
end % value2cell()
- function [p, q, parameters, gridSize] = inputNormalizer(a, b)
-
- parameters = [];
-
- function s = innerInputNormalizer(A, B)
- if isa(A, 'double')
- s.name = num2str(A(:)');
- s.gridName = {};
- s.grid = {};
- s.value = A;
- s.size = size(A);
- s.nargin = 0;
- elseif isa(A, 'quantity.Discrete')
- s.name = A(1).name;
- s.value = A.on();
- s.gridName = A(1).gridName;
- s.grid = A(1).grid;
- s.size = size(A);
- s.nargin = A.nargin();
- parameters = struct(A);
- gridSize = A.domain.gridLength();
- end
- end
-
- q = innerInputNormalizer(b, a);
- p = innerInputNormalizer(a, b);
-
- end
end %% (Static)
+
methods(Access = protected)
function [valDiscrete] = expandValueDiscrete(obj, newDomain)
@@ -2702,13 +2820,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
[sprintf('%ix', obj(1).domain.gridLength, size(obj)) sprintf('\b')];
end
end % getPropertyGroups
+
end % methods (Access = protected)
- methods ( Static, Access = protected )
- function doNotCopy = doNotCopyPropertiesName()
- % #DEPRECATED
- doNotCopy = {'valueContinuous', 'valueDiscrete', 'derivatives'};
- end
- end % methods ( Static, Access = protected )
-
end % classdef
diff --git a/+quantity/Domain.m b/+quantity/Domain.m
index 1b3bb92cf0144a483c3708e697489eb642faef08..187f6eac0eee8cb217e9aca186d7f80adcb059bf 100644
--- a/+quantity/Domain.m
+++ b/+quantity/Domain.m
@@ -63,9 +63,27 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
function s = string(obj)
s = obj.name;
- end
+ end
+
+ function quan = Discrete(obj)
+ % DISCRETE casts the domain into a quantity.Discrete with the values of obj.grid.
+ arguments
+ obj (1, 1) quantity.Domain;
+ end
+ quan = quantity.Discrete(obj.grid, obj, 'name', obj.name);
+ end % quantity.Discrete()
+
+ function quan = Symbolic(obj)
+ % SYMBOLIC casts the domain into a quantity.Symbolic with the grid name as symbolic
+ % value.
+ arguments
+ obj (1, 1) quantity.Domain;
+ end
+ quan = quantity.Symbolic(sym(obj.name), obj, 'name', obj.name);
+ end % Symbolic()
+
end
-
+
methods (Access = public)
function d = split(obj, splittingPoints, optArgs)
@@ -140,15 +158,29 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
end
end
- function d = copyAndRename(obj, newName)
- d = obj.copy();
+ function d = rename(obj, newName, oldName)
+ % RENAME renames the name of the domain OBJ.
+ %
+ % d = rename(obj, newName) replaces all names of obj with the names specified by the
+ % string array newName
+ %
+ % d = rename(obj, newName, oldName) replaces the names of obj specified by the string
+ % array oldName with the corresponding elements of the string array newName
- newName = misc.ensureIsCell(newName);
- assert(length(newName) == length(d), 'For the renaming of a copy of a domain, the numebr of new names must be equal the number of the domains');
- for k = 1:length(newName)
- d(k).name = newName{k};
- end
- end
+ arguments
+ obj quantity.Domain;
+ newName string;
+ oldName string = [obj.name];
+ end % arguments
+
+ assert(numel(newName) == numel(oldName), ...
+ "number of new names must be equal to number of old names.");
+ d = copy(obj);
+ for it = 1 : numel(newName)
+ d(obj.gridIndex(oldName(it))).name = newName{it};
+ end % for it = 1 : numel(newName)
+ assert(misc.isunique([d.name]), "the resulting domain must contain unique names");
+ end % rename
function nd = ndims(obj)
%NDIMS number of dimensions of the domain specified by the
@@ -159,7 +191,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
function n = numGridElements(obj)
% NUMGRIDLEMENTS returns the number of the elements of the grid
n = prod([obj.n]);
- end
+ end % numGridElements()
function s = gridLength(obj)
%GRIDLENGTH number of discretization points for each grid in the
@@ -191,6 +223,19 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
d = obj( obj.gridIndex( searchName ) );
end
+
+ function rest = remove(obj, toBeRemoved)
+ % REMOVE removes the domain specified by toBeRemoved from the quantity.Domain array obj.
+ %
+ % rest = remove(obj, toBeRemoved) removes the domain specified by toBeRemoved from the
+ % quantity.Domain array obj. toBeRemoved may be a string or a domain or a array of
+ % those.
+
+ [~, logicalIdxToBeRemoved] = obj.gridIndex(toBeRemoved);
+
+ rest = obj(~logicalIdxToBeRemoved);
+ end % remove()
+
function [idx, log] = gridIndex(obj, searchName, position)
%% GRIDINDEX returns the index of the grid
% [idx, log] = gridIndex(obj, names) searches in the name
diff --git a/+quantity/Function.m b/+quantity/Function.m
index 9cb0182d035e2221c16054aec3396b992c3c0579..a0b4c1d435e9681a08b1b7307b1a379edc776e85 100644
--- a/+quantity/Function.m
+++ b/+quantity/Function.m
@@ -66,7 +66,7 @@ classdef Function < quantity.Discrete
myParser.addParameter('name', obj(1).name);
myParser.parse(varargin{:});
Q = quantity.Discrete(obj.on(), myParser.Results.domain, 'name', myParser.Results.name);
- end
+ end % quantity.Discrete()
%-----------------------------
% --- overloaded functions ---
@@ -91,14 +91,29 @@ classdef Function < quantity.Discrete
end
value = reshape( cell2mat(value), [ gridLength(myDomain), size(obj)]);
end
- end
+ end % obj2value()
+
+ function value = at(obj, point)
+ % at() evaluates the object at one point and returns it as array
+ % with the same size as size(obj).
+ value = zeros(size(obj));
+ if ~iscell(point)
+ for it = 1 : numel(value)
+ value(it) = obj(it).valueContinuous(point);
+ end
+ else
+ for it = 1 : numel(value)
+ value(it) = obj(it).valueContinuous(point{:});
+ end
+ end
+ end % at()
function n = nargin(obj)
% FIXME: check if all funtions in this object have the same
% number of input values.
n = numel(obj(1).grid);
- end
- end
+ end % nargin()
+ end % methods
% -------------
% --- Mathe ---
diff --git a/+quantity/Symbolic.m b/+quantity/Symbolic.m
index 380276572c3ff07f8b47bbc1895551f814a12004..e5b4c95e1375008fe5f9f9480178231074130a94 100644
--- a/+quantity/Symbolic.m
+++ b/+quantity/Symbolic.m
@@ -4,9 +4,7 @@ classdef Symbolic < quantity.Function
valueSymbolic sym;
variable sym;
end
- properties (Constant)
- defaultSymVar = sym('z', 'real');
- end
+
properties
symbolicEvaluation = false;
end
@@ -30,7 +28,7 @@ classdef Symbolic < quantity.Function
end
variableParser = misc.Parser();
- variableParser.addParameter('variable', quantity.Symbolic.getVariable(valueContinuous));
+ variableParser.addParameter('variable', []);
variableParser.addParameter('symbolicEvaluation', false);
variableParser.parse(varargin{:});
@@ -43,7 +41,6 @@ classdef Symbolic < quantity.Function
fun = cell(S{:});
symb = cell(S{:});
-
for k = 1:numel(valueContinuous)
if iscell(valueContinuous)
@@ -256,13 +253,13 @@ classdef Symbolic < quantity.Function
c = cat@quantity.Discrete(dim, objCell{:});
end % cat()
- function solution = subs(obj, gridName, values)
- %% SUBS substitues gridName and values
+ function solution = subs(obj, gridName2Replace, values)
+ % SUBS substitues gridName and values
% solution = subs(obj, DOMAIN) can be used to change the grid
% of this quantity. The domains with the same domain name as
% DOMAIN will be replaced with the corresponding grid in
% DOMAIN.
- %
+ %J
% solution = subs(obj, DOMAINNAME, GRID) can be used to
% change the grid of this quantity. The grid of the domains with the
% DOMAINNAME will be replaced by the values specified in
@@ -275,124 +272,82 @@ classdef Symbolic < quantity.Function
% quantity.Domain objects. The NEWDOMAIN must be an
% object-array of quantity.Domain.
- if isa(gridName, 'quantity.Domain')
+ if isa(gridName2Replace, 'quantity.Domain')
if nargin == 2
- gridName = {gridName.name};
- values = {gridName.grid};
+ gridName2Replace = {gridName2Replace.name};
+ values = {gridName2Replace.grid};
else
- gridName = {gridName.name};
+ gridName2Replace = {gridName2Replace.name};
end
end
- gridName = misc.ensureString(gridName);
+ gridName2Replace = misc.ensureString(gridName2Replace);
if nargin == 3 && isa(values, 'quantity.Domain')
% replacement of the grid AND the gridName
% 1) replace the grid
- solution = obj.changeGrid( {values.grid}, gridName );
+ solution = obj.changeGrid( {values.grid}, gridName2Replace );
% 2) replace the name
- solution = solution.subs( gridName, {values.name} );
+ solution = solution.subs( gridName2Replace, {values.name} );
return
- else
- values = misc.ensureIsCell(values);
- for it = 1 : numel(values)
- % for later strcmp usage, values must be a cell-array with only
- % numeric or char elements.
- if isstring(values{it})
- values{it} = char(values{it});
- end
- end % for it = 1 : numel(values)
end
- isNumericValue = cellfun(@isnumeric, values);
- if any((cellfun(@(v) numel(v(:)), values)>1) & isNumericValue)
- error('only implemented for one value per grid');
- end
- numericValues = values(isNumericValue);
- if ~isempty(numericValues) && numel(obj(1).gridName) == numel(numericValues)
- % if all grids are evaluated, solution is a double array.
- solution = reshape(obj.on(numericValues), size(obj));
- else
- % evaluate numeric values
- subsGrid = obj(1).grid;
- selectRemainingGrid = false(1, numel(obj(1).grid));
- for currentGridName = gridName(isNumericValue)
- selectGrid = strcmp(obj(1).gridName, currentGridName);
- subsGrid{selectGrid} = values{strcmp(gridName, currentGridName)};
- selectRemainingGrid = selectRemainingGrid | selectGrid;
- end
- newGrid = obj(1).grid(~selectRemainingGrid);
- newGridName = obj(1).gridName(~selectRemainingGrid);
- for it = 1 : numel(values)
- if ~isNumericValue(it) && ~isempty(obj(1).gridName(~selectRemainingGrid))
- % check if there is a symbolic value and if this value exists in the object
- if isstring( values{it} ) || ischar( values{it} )
- newGridName( strcmp(obj(1).gridName(~selectRemainingGrid), gridName(it)) ) ...
- = values{it};
- elseif isa(values{it}, 'sym')
- gridNameTemp = symvar(values{it});
- assert(numel(gridNameTemp) == 1, 'replacing one gridName with 2 gridName is not supported');
- newGridName{strcmp(obj(1).gridName(~selectRemainingGrid), gridName(it))} ...
- = char(gridNameTemp);
-
- else
- error(['value{', num2str(it), '} is of unsupported type']);
- end
- end
+ values = misc.ensureIsCell(values);
+ isNumericValue = false(numel(values), 1);
+ for it = 1 : numel(values)
+ % check which values are substitutions with numerical values in the domain
+ if isnumeric(values{it})
+ isNumericValue(it) = true;
end
-
- symbolicSolution = subs(obj.sym(), cellstr(gridName), values);
+ end % for it = 1 : numel(values)
+
+ % new domain is the old one without the domains that are evaluated at specific numeric
+ % values. Furthermore, the other domains are renamed.
+ newDomain = obj(1).domain.remove(gridName2Replace(isNumericValue));
+
+ if any(~isNumericValue)
+ % rename domains that are affected by symbolic substitution
+ gridName2ReplaceSymbolicOld = gridName2Replace(~isNumericValue);
+ gridName2ReplaceSymbolicNew = string(values(~isNumericValue));
- if isempty(symvar(symbolicSolution(:)))
- % take the new values that are not numeric as the new
- % variables:
- newGridName = unique(newGridName, 'stable');
- newGrid = cell(1, numel(newGridName));
- charValues = cell(0);
- for it = 1 : numel(values)
- if isstring(values{it}) || ischar(values{it})
- charValues{end+1} = char(values{it});
- end
- end
- for it = 1 : numel(newGridName)
- if ~isempty(charValues) && any(contains(charValues, newGridName{it}))
- newGrid{it} = obj.gridOf(gridName( strcmp(newGridName{it}, values) ));
+ resultingName = [newDomain.name];
+ for it = 1 : numel(gridName2ReplaceSymbolicOld)
+ resultingName(strcmp([newDomain.name], gridName2ReplaceSymbolicOld(it))) ...
+ = gridName2ReplaceSymbolicNew(it);
+ end % for it = 1 : numel(gridName2ReplaceSymbolicOld)
+ if misc.isunique(resultingName)
+ newDomain = newDomain.rename(gridName2ReplaceSymbolicNew, gridName2ReplaceSymbolicOld);
+ else % complicated case, for instance f(x, y).subs("x", "y") = f(y, y)
+ % for this, join is used.
+ for it = 1 : numel(gridName2ReplaceSymbolicNew)
+ if any(strcmp([newDomain.name], gridName2ReplaceSymbolicNew(it)))
+ % remove domain to be removed and join it with the single renamed domain.
+ newDomain = join(newDomain.remove(gridName2ReplaceSymbolicOld(it)), ...
+ newDomain.find(gridName2ReplaceSymbolicOld(it)).rename(...
+ gridName2ReplaceSymbolicNew(it)));
else
- newGrid{it} = obj.gridOf(newGridName{it});
- end
- end
- solution = quantity.Symbolic(double(symbolicSolution), ...
- 'grid', newGrid, 'gridName', newGridName, 'name', obj(1).name);
- else
- % before creating a new quantity, it is checked that
- % newGridName is unique. If there are non-unique
- % gridName, multiple are removed and the finest grid
- % from newGrid is taken.
- uniqueGridName = unique(newGridName, 'stable');
- if numel(newGridName) == numel(uniqueGridName)
- solution = quantity.Symbolic(symbolicSolution, ...
- 'grid', newGrid, 'gridName', newGridName, 'name', obj(1).name);
- else
- uniqueGrid = cell(1, numel(uniqueGridName));
- for it = 1 : numel(uniqueGrid)
- gridCandidatesTemp = newGrid(...
- strcmp(uniqueGridName{it}, newGridName));
- for jt = 1 : numel(gridCandidatesTemp)
- if ~isempty(uniqueGrid{it})
- assert(uniqueGrid{it}(1) == gridCandidatesTemp{jt}(1) ...
- && uniqueGrid{it}(end) == gridCandidatesTemp{jt}(end), ...
- 'grids must have same domain');
- end
- if isempty(uniqueGrid{it}) || ...
- numel(gridCandidatesTemp{jt}) > numel(uniqueGrid{it})
- uniqueGrid{it} = gridCandidatesTemp{jt};
- end
- end
+ newDomain = newDomain.rename(...
+ gridName2ReplaceSymbolicNew(it), gridName2ReplaceSymbolicOld(it));
end
- solution = quantity.Symbolic(symbolicSolution, ...
- 'grid', uniqueGrid, 'gridName', uniqueGridName, 'name', obj(1).name);
- end
+ end % for it = 1 : numel(gridName2ReplaceSymbolicNew)
+
end
+
+ % convert elements of values to symbolic expressions
+ for it = 1 : numel(values)
+ if ~isnumeric(values{it})
+ values{it} = sym(values{it});
+ end
+ end % for it = 1 : numel(values)
+ end
+
+ if numel(newDomain) > 0
+ solution = quantity.Symbolic(...
+ subs(obj.sym(), sym(gridName2Replace), values), ...
+ newDomain, ...
+ 'name', obj(1).name);
+ else
+ solution = double(subs(obj.sym(), sym(gridName2Replace), values));
end
end % subs()
@@ -488,8 +443,8 @@ classdef Symbolic < quantity.Function
obj(it).variable, varNew(s)), varNew(0)==ic);
end
solution = quantity.Symbolic(symbolicSolution, ...
- 'gridName', {myParser.Results.newGridName, 'ic'}, ...
- 'grid', {myParser.Results.variableGrid, myParser.Results.initialValueGrid}, ...
+ [quantity.Domain(myParser.Results.newGridName, myParser.Results.variableGrid), ...
+ quantity.Domain('ic', myParser.Results.initialValueGrid)], ...
'name', "solve(" + obj(1).name + ")");
end % solveDVariableEqualQuantity()
@@ -540,6 +495,19 @@ classdef Symbolic < quantity.Function
P = prod(quantity.Discrete(obj), dim);
end
end % prod()
+
+ function P = power(obj, p)
+ % a.^p elementwise power
+ P = quantity.Symbolic(power(obj.sym(), p), obj(1).domain, ...
+ "name", obj(1).name + ".^{" + num2str(p) + "}");
+ end % power()
+
+ function P = mpower(obj, p)
+ % Matrix power a^p is matrix or scalar a to the power p.
+ P = quantity.Symbolic(mpower(obj.sym(), p), obj(1).domain, ...
+ "name", obj(1).name + "^{" + num2str(p) + "}");
+ end % mpower()
+
function d = det(obj)
% det(X) returns the the determinant of the squre matrix X
if isempty(obj)
@@ -710,6 +678,22 @@ classdef Symbolic < quantity.Function
'name', "expm(" + obj(1).name + ")");
end % expm()
+ function y = log(obj)
+ % log Natural logarithm.
+ % log(X) is the natural logarithm of the elements of X.
+ % Complex results are produced if X is not positive.
+ y = quantity.Symbolic(log(obj.sym()), obj(1).domain, ...
+ 'name', "log(" + obj(1).name + ")");
+ end % log()
+
+ function y = log10(obj)
+ % log10 Common (base 10) logarithm.
+ % log10(X) is the base 10 logarithm of the elements of X.
+ % Complex results are produced if X is not positive.
+ y = quantity.Symbolic(log10(obj.sym()), obj(1).domain, ...
+ 'name', "log10(" + obj(1).name + ")");
+ end % log10()
+
function y = ctranspose(obj)
% ctranspose() or ' is the complex conjugate tranpose
y = quantity.Symbolic(conj(obj.sym().'), obj(1).domain, ...
@@ -846,65 +830,31 @@ classdef Symbolic < quantity.Function
P = quantity.Symbolic(zeros(valueSize), varargin{:});
end
- function [p, q, parameters] = inputNormalizer(a, b)
-
- parameters = [];
-
- function s = innerInputNormalizer(a)
- if isa(a, 'double')
- s.name = num2str(a(:)');
- s.value = a;
- % s.variable = [];
- elseif isa(a, 'sym')
- s.name = char(a);
- s.value = a;
- % s.variable = [];
- elseif isa(a, 'quantity.Symbolic')
- s.name = a(1).name;
- s.value = a.sym();
- % s.variable = a.variable;
- parameters = struct(a);
- parameters = rmfield(parameters, 'gridName');
- end
- end
-
- p = innerInputNormalizer(a);
- q = innerInputNormalizer(b);
-
- end
- end
+ end % methods (Static)
methods (Static, Access = protected)
function var = getVariable(symbolicFunction, nVar)
+ % nVar is needed for overloaded functionality in PolynomialOperator
if misc.issym(symbolicFunction(:).')
var = symvar(symbolicFunction(:).');
- else
- var = [];
- end
-
- if isempty(var)
- % var = quantity.Symbolic.defaultSymVar;
- else
-
- if var(1) ~= quantity.Symbolic.defaultSymVar
+ if ~isempty(var) && ~strcmp(string(var(1)), "z")
% fast solution to order the spatial and temporal variable
% in the common order: (z, t)
var = flip(var);
end
+ else
+ var = [];
end
- end
- function doNotCopy = doNotCopyPropertiesName()
- parentalProperties = quantity.Discrete.doNotCopyPropertiesName;
- doNotCopy = {parentalProperties{:}, ...
- 'valueSymbolic', 'isConstant'};
- end
- function f = setValueContinuous(f, symVar)
+ end % getVariable()
+
+ function [f, argList] = setValueContinuous(f, symVar)
% converts symVar into a matlabFunction.
+ argList = string(symVar);
if isa(f, 'sym')
if isempty(symvar(f))
fDouble = double(f);
- f = @(varargin) fDouble + quantity.Function.zero(varargin{:});
+ f = eval("@(" + argList + ") fDouble + quantity.Function.zero(varargin{:})");
else
if iscell( symVar )
symVar = [symVar{:}];
@@ -914,5 +864,5 @@ classdef Symbolic < quantity.Function
end
end
end % setValueContinuous()
- end
-end
+ end % methods (Static, Access = protected)
+end % classdef
diff --git a/+signals/+faultmodels/TaylorPolynomial.m b/+signals/+faultmodels/TaylorPolynomial.m
index d451e9b9e72eb5b671f1eb77a1955beb8431c005..a98619af394b66dd586fa76b9d344db04da496a4 100644
--- a/+signals/+faultmodels/TaylorPolynomial.m
+++ b/+signals/+faultmodels/TaylorPolynomial.m
@@ -79,8 +79,7 @@ classdef TaylorPolynomial < quantity.Discrete
occurrence);
obj@quantity.Discrete( val, ...
- 'gridName', {'t'}, ...
- 'grid', {myGlobalGrid}, ...
+ quantity.Domain("t", myGlobalGrid), ...
'name', ['f_' num2str(index)]);
% put all the values to local variables
@@ -128,12 +127,11 @@ classdef TaylorPolynomial < quantity.Discrete
end
function monomials = getMonomials(degree, t, tStar, myLocalGrid, T)
monomials = quantity.Symbolic.zeros([degree, 1], ...
- 'domain', quantity.Domain.empty(), ...
+ quantity.Domain.empty(), ...
'name', 'varphi_k');
for k = 1:degree
monomials(k, 1) = quantity.Symbolic((t-tStar)^(k-1), ...
- 'gridName', cellfun(@char, {t, tStar}, 'UniformOutput', false), ...
- 'grid', {myLocalGrid, T}, ...
+ [quantity.Domain(char(t), myLocalGrid), quantity.Domain(char(tStar), T)], ...
'name', 'varphi_k');
end
end
diff --git a/+signals/BasicVariable.m b/+signals/BasicVariable.m
index de37d3850c89f4013f76549accff985d302920b5..f259d2625d13e4c91fc2302db8ba635022b0f43d 100644
--- a/+signals/BasicVariable.m
+++ b/+signals/BasicVariable.m
@@ -1,39 +1,239 @@
classdef BasicVariable < handle
- %BasicVariable Summary of this class goes here
- % Detailed explanation goes here
+ %BasicVariable Class to handle quantities with a priori known derivatives
+ % This class is used to simplify the usage of quantities for which a lot of computations must be
+ % done on its derivatives. So a BasicVariable consists of a basic function and its derivatives,
+ % which are precomputed and safed. So computations can be performed without recalculation of the
+ % derivatives.
properties ( SetAccess = protected )
- derivatives cell;
- fun quantity.Discrete;
+ derivatives cell;
+ fun (1,1) quantity.Discrete;
end
- methods
+ properties ( Dependent )
+ highestDerivative; % order of the highest computed derivative
+ T; % the transition time T
+ dt; % step size
+ end
+
+ methods
function obj = BasicVariable(fun, derivatives)
obj.fun = fun;
obj.derivatives = derivatives;
end
+ function h = get.highestDerivative(obj)
+ h = numel(obj(1).derivatives);
+ end
+ function T = get.T(obj)
+ T = obj.domain.upper;
+ end
+ function dt = get.dt(obj)
+ if isa(obj.domain, "quantity.EquidistantDomain")
+ dt = obj.domain.stepSize;
+ else
+ delta_t = diff(obj.domain.grid);
+ assert( all( diff( delta_t ) < 1e-15 ), "Grid is not equidistant spaced" );
+ dt = delta_t(1);
+ end
+ end
end
methods ( Access = public)
- function D = diff(obj, k, ~)
- % DIFF get the k-th derivative from saved derivatives
+ function D = diff(obj, domain, n)
+ % DIFF get the k-th derivative
+ % D = diff(obj, domain, n) tries to find the n-th derivative of obj.fun in the saved
+ % derivatives. If it does not exist, it tries to compute it with the evaluateFunction.
+ % If n is a non integer number, the Riemann-Liouville fractional derivative is
+ % computed. For this, a formula is used which requires the next integer + 1 derivative
+ % of the function to avoid the singularity under the integral. Actually, the domain is
+ % clear, but it is used hear to be consistent with the quantity.Discrete/diff.
arguments
obj
- k double {mustBeInteger} = 0;
- ~;
+ domain (1,1) quantity.Domain;
+ n (1,1) double = 0;
+ end
+
+ for l = 1:numel(obj)
+
+ assert( domain.isequal( obj(l).domain ),...
+ "The derivative wrt. this domain is not possible.")
+
+ if n == 0
+ D(l) = obj(l).fun;
+
+ elseif ~numeric.near(round(n), n)
+ % fractional derivative:
+ D(l) = obj(l).fDiff(n);
+
+ elseif n > length(obj(l).derivatives)
+ % For some functions their maybe a possibility to compute new derivatives on the
+ % fly. For this, the evaluateFunction(k) can be used, where k is the order of
+ % the derivative.
+ D(l) = obj(l).evaluateFunction(obj.domain, n);
+ obj(l).derivatives{end+1} = D(l);
+ else
+ % use the pre-computed derivatives
+ D(l) = obj(l).derivatives{n};
+ end
+ end
+ D = reshape(D, size(obj));
+ end
+
+ function D = diffs(obj, k)
+ % DIFFS compute multiple derivatives
+ % D = diffs(obj, k) computes the k-th derivatives of OBJ. At this, k can be a vector
+ % of derivative orders.
+ arguments
+ obj,
+ k (:,1) double {mustBeInteger, mustBeNonnegative} = 0:numel(obj(1).derivatives);
end
+ D_ = arrayfun( @(i) obj.diff(obj.domain, i), k, 'UniformOutput', false);
- if k == 0
- D = obj.fun;
- elseif k > length(obj.derivatives)
- error(['Requested derivative ' num2str(k) ' is not available']);
+ for i = 1:numel(k)
+ D(i,:) = D_{i};
+ end
+ end
+
+
+
+ function h = productRule(obj, phi, n)
+ % LEIBNIZ compute the k-derivative of the product of the basic variable and phi
+ % h = productRule(obj, phi, n) computes the k-th derivative of obj.fun*phi, i.e.,
+ % h = d_t^n ( obj.fun(t) * phi )
+ % using the product rule. This is usefull if the derivatives of the basic variable obj
+ % are known exactly and the derivatives of phi, but not the multiplication
+ % obj.diff(0)*phi.
+ t = phi.domain;
+ h = quantity.Discrete.zeros( size(obj), t );
+
+ if numeric.near(round(n), n)
+ % for integer order derivatives
+ for k = 0:n
+ h = h + misc.binomial(n, k) * obj.diff(t, n-k) * phi.diff(t, k);
+ end
else
- D = obj.derivatives{k};
+ % for fractional order derivatives: [see Podlubny]
+ % p \~~ oo / p \ (k) p-k
+ % D (phi(t) f(t)) = > | | | phi (t) D f(t)
+ % a t /__ k=0 \ k / a t
+ h = misc.binomial( n, 0 ) * phi * obj.fDiff(n);
+ h0 = h;
+ k = 1;
+ % do the iteration on the loop until the difference is very small.
+ while max( abs( h0 ) ) / max( abs( h ) ) > 1e-9
+ h0 = misc.binomial(n,k) * phi.diff(t, k) * obj.fDiff( n - k);
+ h = h + h0;
+ k = k+1;
+ end
end
end
+
function d = domain(obj)
- d = obj.fun(1).domain;
+ d = obj(1).fun(1).domain;
end
- end
-end
\ No newline at end of file
+
+ function plot(obj, varargin)
+ F = quantity.Discrete(obj);
+ F.plot(varargin{:});
+ end
+
+ function F = quantity.Discrete(obj)
+ F (1,:) = [obj.fun];
+
+ for i = 1:max([obj.highestDerivative])
+
+ for j = 1:numel(obj)
+ F_i(j) = obj(j).derivatives{i};
+ end
+
+ F(i+1,:) = reshape(F_i, size(obj));
+ end
+ end
+
+ function c = mtimes(a, b)
+ arguments
+ a double;
+ b signals.BasicVariable
+ end
+
+ assert( size(a,2) == size(b,1), "dimensions of the terms for multiplication do not match")
+ F = quantity.Discrete(b);
+
+ % do the multiplication for each derivative:
+ for i = 1:max([b.highestDerivative])+1
+ tmp(i,:) = a * shiftdim(F(i,:));
+ end
+
+ % restore the result as signals.BasicVariable
+ for k = 1:numel(b)
+ for l = 2:size(tmp,1)
+ derivs_{l-1} = tmp(l,k);
+ end
+ c(k) = signals.BasicVariable( tmp(1,k), derivs_);
+ end
+ c = reshape(c, size(b));
+ end
+
+ end % methods (Access = public)
+
+ methods (Access = protected)
+ function v = evaluateFunction(obj, domain, k)
+ error('conI:signals:BasicVariable:derivative', ['Requested derivative ' num2str(k) ' is not available']);
+ end
+
+ function D = fDiff(obj, p)
+ % FDIFF compute the p-fractional derivative
+ % D = fDiff(obj, p) computes the fractional derivative of order p. To be
+ % concrete the Riemann-Liouville fractional derivative with initial point t0 = lower
+ % end of the considered time domain is computed.
+ % Because the basisc variable must be C^infinity function, we use the form
+ % of the fractional derivative described in eq. (2.80) in the book [Podlubny: Fractional
+ % differential equations]. Then, integration by parts is applied to eliminate the
+ % singularity:
+ %
+
+ n = ceil(p);
+ tDomain = obj(1).domain;
+ t0 = tDomain.lower;
+ t = Discrete( tDomain );
+ tauDomain = tDomain.rename( tDomain.name + "_Hat" );
+ tau = Discrete( tauDomain );
+ if p > 0
+
+ % TODO Try to implement the polynomial part as function or as symbolic
+ D = quantity.Discrete.zeros(size(obj), tDomain);
+ for k = 0:n
+ f_dk = obj.diff(tDomain, k);
+
+ if f_dk.at(t0) ~= 0
+ D = D + f_dk.at(t0) * (t - t0).^(k-p) / gamma( k-p + 1);
+ end
+
+ end
+ assert( ~ any( isinf( obj.diff(tDomain, n+1).on() ) ), "conI:BasicVariable:fDiff", ...
+ "The function has a singularity! Hence, the fractional derivative can not be computed")
+ D = D + cumInt( ( t - tau ).^(n-p) * subs(obj.diff(tDomain, n+1), tDomain.name, tauDomain.name), ...
+ tauDomain, tauDomain.lower, tDomain.name) / gamma( n+1-p );
+
+ elseif p < 0 && p > -1
+ alpha = -p;
+ D = (t - t0).^(alpha) / gamma( alpha + 1) * obj.diff(tDomain, 0).at(t0) + ...
+ cumInt( (t-tau).^alpha * subs(obj.diff(tDomain, 1), tDomain.name, tauDomain.name), ...
+ tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha + 1 );
+
+ elseif p <= -1
+ alpha = -p;
+ D = cumInt( (t-tau).^(alpha-1) * subs( obj.diff(tDomain, 0), tDomain.name, tauDomain.name), ...
+ tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha );
+
+ elseif p == 0
+ D = [obj.fun()];
+ else
+ error("This should not happen")
+ end
+ end
+
+ end
+
+end
diff --git a/+signals/GevreyFunction.m b/+signals/GevreyFunction.m
index bc1b546485e4678d3e5941ff3c0679e8943b1bd0..49ae6323fd7f7d41f40fcd8537a727b79019d96e 100644
--- a/+signals/GevreyFunction.m
+++ b/+signals/GevreyFunction.m
@@ -1,47 +1,91 @@
-classdef GevreyFunction < quantity.BasicVariable
- %UNTITLED Summary of this class goes here
- % Detailed explanation goes here
+classdef GevreyFunction < signals.BasicVariable
+ %GEVREYFUNCTION -
+ % Class to generate a gevrey function and its derivative. The general function has the form
+ %
+ % / 0 : t <= 0
+ % |
+ % f(t) = | f0 + K * h(t) : 0 < t < T
+ % |
+ % \ 0 : t >= T
+ %
+ % with the additional property int_0^T h(t) = 1.
properties
% Gevrey-order
- order double = -1;
+ order (1,1) double = 1.99;
+ % the constant gain
+ gain double = 1;
+ % the offset:
+ offset (1,1) double = 0;
+ % the number of derivatives to be considered
+ numDiff (1,1) double = 5;
+ % underlying gevery function
+ g;
+ % derivative shift
+ diffShift (1,1) double = 0;
end
properties (Dependent)
+ % the exponent of the exponential function
sigma;
end
methods
- function obj = GevreyFunction(varargin)
+ function obj = GevreyFunction(optArgs)
+ % GEVREYFUNCTION initialization of a gevrey function as signals.BasicVariable
+ % obj = GevreyFunction(optArgs) creates a gevrey function object. The properties can
+ % be set by name-value-pairs and are "order", "gain", "offset", "numDiff",
+ % "timeDomain", "g", "diffShift".
+ arguments
+ optArgs.order (1,1) double = 1.99;
+ optArgs.gain double = 1;
+ optArgs.offset (1,1) double = 0;
+ optArgs.numDiff (1,1) double = 5;
+ optArgs.timeDomain (1,1) quantity.Domain = quantity.Domain.defaultDomain(100, "t");
+ optArgs.g = @signals.gevrey.bump.g;
+ optArgs.diffShift (1,1) double = 0;
+ end
- if nargin == 0
- g = [];
- else
- g = @signals.gevrey.bump.g;
+ % generate the underlying gevrey functions:
+ f_derivatives = cell(optArgs.numDiff + 1, 1);
+ for k = 0:optArgs.numDiff
+ f_derivatives{k+1} = ...
+ signals.GevreyFunction.evaluateFunction_p(optArgs, optArgs.timeDomain, k);
end
- obj@quantity.BasicVariable(g, varargin{:});
- if nargin > 0
- prsr = misc.Parser();
- prsr.addParameter('order', 1.5 );
- prsr.parse(varargin{:});
- obj.order = prsr.Results.order;
+
+ obj@signals.BasicVariable(f_derivatives{1}, f_derivatives(2:end));
+
+ % set the parameters to the object properties:
+ for k = fieldnames(rmfield(optArgs, "timeDomain"))'
+ obj.(k{:}) = optArgs.(k{:});
end
+
end
-
-
function s = get.sigma(obj)
s = 1 ./ ( obj.order - 1);
end
function set.sigma(obj, s)
obj.order = 1 + 1/s;
end
-
+ end
+
+ methods (Static, Access = private)
+ function v = evaluateFunction_p(optArgs, domain, k)
+
+ derivativeOrder = k + optArgs.diffShift;
+ % The offset is only added to the zero-order derivative.
+ offset_k = optArgs.offset * 0^derivativeOrder;
+
+ v = quantity.Function(@(t) offset_k + optArgs.gain * ...
+ optArgs.g(t, derivativeOrder, ...
+ domain.upper, 1 ./ ( optArgs.order - 1)), domain);
+ end
end
methods (Access = protected)
- function v = evaluateFunction(obj, grid)
- v = obj.K * obj.valueContinuous(grid, obj.diffShift, obj.T, obj.sigma);
+ function v = evaluateFunction(obj, domain, k)
+ v = signals.GevreyFunction.evaluateFunction_p(obj, domain, k);
end
end
end
diff --git a/+signals/PolyBump.m b/+signals/PolyBump.m
index a7087598e4a1072348bd9401b90290040eb5eef6..e841e98328bb737c5950b04ed5f99faa92183c33 100644
--- a/+signals/PolyBump.m
+++ b/+signals/PolyBump.m
@@ -40,7 +40,7 @@ classdef PolyBump < quantity.Function
fun = signals.PolyBump.polyBumpFunction(a, b, c, T0, T1);
- obj@quantity.Function(fun, 'domain', t);
+ obj@quantity.Function(fun, t);
obj.a = a;
obj.b = b;
diff --git a/+signals/PolynomialOperator.m b/+signals/PolynomialOperator.m
index 43c45824587e28b3f5acaffb8015e013b800bbea..1558cb401cdee65e5d03ded4a173bc97c8b27d16 100644
--- a/+signals/PolynomialOperator.m
+++ b/+signals/PolynomialOperator.m
@@ -69,13 +69,13 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
for k = 1 : numel(A)
if isnumeric(A{k})
obj(k).coefficient = ...
- quantity.Discrete(A{k}, 'domain', prsr.Results.domain);
+ quantity.Discrete(A{k}, prsr.Results.domain);
else
obj(k).coefficient = A{k}; %#ok<AGROW>
end
end
elseif isnumeric(A)
- obj(1).coefficient = quantity.Discrete(A, 'domain', prsr.Results.domain);
+ obj(1).coefficient = quantity.Discrete(A, prsr.Results.domain);
end
[obj.s] = deal(s);
@@ -200,7 +200,7 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
c = cell(1, size(C, 3));
for k = 1:size(C, 3)
c{k} = quantity.Discrete(double(C(:,:,k)), ...
- 'gridName', {}, 'grid', {}, 'name', "adj(" + obj(1).coefficient(1).name + ")");
+ quantity.Domain.empty(), 'name', "adj(" + obj(1).coefficient(1).name + ")");
end
C = signals.PolynomialOperator(c);
@@ -214,7 +214,7 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
c = cell(1, size(C, 3));
for k = 1:size(C, 3)
c{k} = quantity.Discrete(double(C(:,:,k)), ...
- 'domain', quantity.Domain.empty(), ...
+ quantity.Domain.empty(), ...
'name', "det(" + A(1).coefficient(1).name + ")");
end
@@ -363,10 +363,9 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
end
f = quantity.Discrete(...
misc.fundamentalMatrix.odeSolver_par( A_full, spatialDomain ), ...
- 'grid', {spatialDomain, spatialDomain}, 'gridName', {'z', 'zeta'});
+ [quantity.Domain("z", spatialDomain), quantity.Domain("zeta", spatialDomain)]);
- b = quantity.Discrete(B_full, 'grid', spatialDomain, ...
- 'gridName', 'zeta');
+ b = quantity.Discrete(B_full, quantity.Domain("zeta", spatialDomain));
p = cumInt(f * b, 'zeta', spatialDomain(1), 'z');
z0 = A(1).grid{1}(1);
@@ -489,6 +488,7 @@ classdef PolynomialOperator < handle & matlab.mixin.Copyable
end
function newOperator = changeGrid(obj, gridNew, gridNameNew)
+ warning("DEPRICATED: Use subs instead")
newOperator = signals.PolynomialOperator(obj.M.changeGrid(gridNew, gridNameNew));
end
diff --git a/+signals/RandomStairs.m b/+signals/RandomStairs.m
index 9dab4c4b4601c8bde95133ac715258754830da49..2d593f8d0d3fee7e9a6bbd874499ef7f85de9b9d 100644
--- a/+signals/RandomStairs.m
+++ b/+signals/RandomStairs.m
@@ -98,11 +98,9 @@ classdef RandomStairs
end
if optArg.smooth
- q = quantity.Discrete(obj.smoothSignal(timeDomain.grid), ...
- 'domain', timeDomain, 'name', obj.name);
+ q = quantity.Discrete(obj.smoothSignal(timeDomain.grid), timeDomain, 'name', obj.name);
else
- q = quantity.Discrete(obj.signal(timeDomain.grid), ...
- 'domain', timeDomain, 'name', obj.name);
+ q = quantity.Discrete(obj.signal(timeDomain.grid), timeDomain, 'name', obj.name);
end
end
diff --git a/+signals/SignalModel.m b/+signals/SignalModel.m
index b46488429ab24461989808b18e626cb8ba1ebe09..8e8916bac3e4ff454b74f9995c0a1e1aac03511c 100644
--- a/+signals/SignalModel.m
+++ b/+signals/SignalModel.m
@@ -128,8 +128,8 @@ classdef SignalModel
result.v( t >= obj(i).occurrence, : ) = sol.v;
end
- r = [r; quantity.Discrete(result.r, 'domain', optArg.time, 'size', [obj(i).n_p, 1])];
- v = [v; quantity.Discrete(result.v, 'domain', optArg.time, 'size', [obj(i).n_v, 1])];
+ r = [r; quantity.Discrete(result.r, optArg.time, 'size', [obj(i).n_p, 1])];
+ v = [v; quantity.Discrete(result.v, optArg.time, 'size', [obj(i).n_v, 1])];
end
diff --git a/+signals/TimeDelayOperator.m b/+signals/TimeDelayOperator.m
index 449f5312cbf9e8d4e463e8ac2d56e4d8e30289d8..ac7a86bb9654a3e9eb4721a4cdf766ef4f1800ff 100644
--- a/+signals/TimeDelayOperator.m
+++ b/+signals/TimeDelayOperator.m
@@ -83,7 +83,7 @@ classdef TimeDelayOperator < handle & matlab.mixin.Copyable
end
function t = t(obj)
- t = quantity.Symbolic(sym(obj.timeDomain.name), 'domain', obj.timeDomain);
+ t = quantity.Symbolic(sym(obj.timeDomain.name), obj.timeDomain);
end
function d = diag(obj)
@@ -129,7 +129,7 @@ classdef TimeDelayOperator < handle & matlab.mixin.Copyable
end
qinf = quantity.Discrete(infValue, ...
- 'domain', namedArgs.spatialDomain);
+ namedArgs.spatialDomain);
o = signals.TimeDelayOperator(qinf, timeDomain, 'isZero', true);
end
end
diff --git a/+signals/WhiteGaussianNoise.asv b/+signals/WhiteGaussianNoise.asv
deleted file mode 100644
index 20d0bce1383cb17e84238449fe095dc29ae93a73..0000000000000000000000000000000000000000
--- a/+signals/WhiteGaussianNoise.asv
+++ /dev/null
@@ -1,78 +0,0 @@
-classdef WhiteGaussianNoise
- %WhiteGaussianNoise
-
- properties
- sigma;
- mue = 0;
- end
-
- properties (Dependent)
- variance;
- end
-
-
- methods
- function obj = WhiteGaussianNoise(sigma)
- obj.sigma = sigma;
- end
-
- function s = get.variance(obj)
- s = obj.sigma^2;
- end
-
- function o = noise(obj, t)
- o = randn(size(t)) * sqrt(obj.variance);
- end
-
- function plotProbabilityDistributionFunction(obj, noise)
- histogram(noise);
- end
-
- function p = probabilityDensityFunction(obj, varargin)
-
- prsr = misc.Parser();
- prsr.addParameter('grid', linspace(- 5* obj.sigma, 5 * obj.sigma)');
- prsr.parse(varargin{:});
-
- x = sym('x', 'real');
- p = quantity.Symbolic( 1 / sqrt( 2 * pi ) / obj.sigma * ...
- exp( - (x - obj.mue)^2 / 2 / obj.sigma^2 ), ...
- 'grid', prsr.Results.grid, ...
- 'variable', x);
- end
-
- function P = cumulativeDistributionFunction(obj, varargin)
- p = obj.probabilityDensityFunction(varargin{:});
- P = p.cumInt('x', p.grid{1}(1), 'x');
- end
-
- function p = probability(obj, z)
- % computes the probability, that a value is inside the domain
- % abs(x) < z
-
- P = obj.cumulativeDistributionFunction;
-
- p = 2 * P.at(abs(z - obj.mue)) - 1;
-
- end
- function
- %CALCULATETREHSHOLDFORPROBABILITY computes the threshold for fault
-%detection for a given probability detection
-%
-% fB = calculateProbabilityforDetection(X, MUE, SIGMA), calculate the
-% treshold fB for a correct detection of the fault given the probability for
-% detection X. Based on a normal distributed random variable with probability
-% function P for a mean value MUE and a standard deviation SIGMA.
-% The probability X is computed by
-% P(|f(t)| <= fB) = P(-fB < f(t) <= fB) - P( |f(t)| <= fB)
-% = F(fB) - F(-fB)
-% F(fB) is the cumulative distribution function for a normal distribution
-% with mean mue and sigma X = F(fB) = 0.5(1+erf((fB-mu)/(sqrt(2)*sigma)))
-% Solving for fB with the inverse erf function and erf(-x) = -erf(x) yields
-% fB = sqrt(2) * sigma * erfinv(X) + mue
-
-fB = sqrt(2)*sigma*erfinv(X) + mue;
-
- end
-end
-
diff --git a/+signals/WhiteGaussianNoise.m b/+signals/WhiteGaussianNoise.m
index dce235f8fb6cff495ac0c5747d4afa0af21d72de..0c33d8c8691ad95689f8139abf8cd49c12c54a5f 100644
--- a/+signals/WhiteGaussianNoise.m
+++ b/+signals/WhiteGaussianNoise.m
@@ -89,16 +89,14 @@ classdef WhiteGaussianNoise
prsr.parse(varargin{:});
x = sym('x', 'real');
+ xDomain = quantity.Domain("x", prsr.Results.grid);
- p = quantity.Symbolic(zeros(size(obj)), ...
- 'grid', prsr.Results.grid, 'variable', x);
+ p = quantity.Symbolic(zeros(size(obj)), xDomain);
for i = 1:numel(obj)
p(i) = quantity.Symbolic( 1 / sqrt( 2 * pi ) / obj(i).sigma * ...
- exp( - (x - obj(i).mue)^2 / 2 / obj(i).sigma^2 ), ...
- 'grid', prsr.Results.grid, ...
- 'variable', x);
+ exp( - (x - obj(i).mue)^2 / 2 / obj(i).sigma^2 ), xDomain);
end
end
diff --git a/+signals/smoothStep.m b/+signals/smoothStep.m
index f649d53a14f4358470664b80468f72d5cb795ed1..c2e7addd6bf0c766d2887649dc567850a131f008 100644
--- a/+signals/smoothStep.m
+++ b/+signals/smoothStep.m
@@ -10,6 +10,9 @@ arguments
optArg.s0 = 0;
optArg.s1 = 1;
end
+%
+
+warning("DEPRICATED: Use misc.ss.setPointChange instead");
assert( all( size(optArg.s0) == [deg, 1]))
assert( all( size(optArg.s1) == [deg, 1]))
@@ -19,6 +22,6 @@ o = zeros(1, deg);
s = spline([optArg.domain.lower, optArg.domain.upper], ...
[o optArg.s0.', optArg.s1.' o], optArg.domain.grid);
-S = quantity.Discrete(s, 'domain', optArg.domain);
+S = quantity.Discrete(s, optArg.domain);
end
\ No newline at end of file
diff --git a/+unittests/+misc/+ss/testSs.m b/+unittests/+misc/+ss/testSs.m
index 633f9f1c54c49613edd77aab148f58751c548f28..924a32d552143ac58021438e36d088370bb8cabc 100644
--- a/+unittests/+misc/+ss/testSs.m
+++ b/+unittests/+misc/+ss/testSs.m
@@ -4,12 +4,37 @@ function [tests] = testSs()
tests = functiontests(localfunctions());
end
+function testSetPointChange(testCase)
+
+time = quantity.Domain('t', linspace(log(2), 17 / 8));
+n = 5;
+
+phi0 = pi;
+phi1 = exp(1);
+
+p = quantity.Symbolic( misc.ss.setPointChange(phi0, phi1, time.lower, time.upper, n), time);
+
+testCase.verifyEqual(p.at(time.lower), phi0, 'AbsTol', 1e-12);
+testCase.verifyEqual(p.at(time.upper), phi1, 'AbsTol', 1e-10);
+
+for k = 1:n
+ p_k = diff(p, 't', k);
+ testCase.verifyEqual(p_k.at(time.lower), 0, 'AbsTol', 1e-7);
+ testCase.verifyEqual(p_k.at(time.upper), 0, 'AbsTol', 1e-7);
+end
+
+p_np1 = diff(p, 't', n + 1);
+testCase.verifyFalse( abs( p_np1.at(time.lower) ) < 1e-6);
+testCase.verifyFalse( abs( p_np1.at(time.upper) ) < 1e-6);
+
+end % testSetPointChange
+
function testRelativeDegree(tc)
mySs = ss([0, 1; 0, 0], [0, 1; 1, 0], [1, 0; 0, 1], []);
tc.verifyEqual(misc.ss.relativeDegree(mySs), [2, 1; 1, NaN]);
tc.verifyEqual(misc.ss.relativeDegree(mySs, "min"), 1);
tc.verifyEqual(misc.ss.relativeDegree(mySs, "max"), 2);
-end
+end % testRelativeDegree
function planPolynomialTrajectory(tc)
%%
@@ -139,19 +164,19 @@ end
function testCombineInputSignals(tc)
myStateSpace = ss(-1, [1, 2, 3], 1, [], 'InputName', {'control', 'disturbance', 'my.fault'});
t = quantity.Domain('t', linspace(0, 4, 201));
-disturbanceSignal = quantity.Symbolic(sin(sym('t')), 'domain', t, 'name', 'disturbance');
-faultSignal = quantity.Symbolic(sin(sym('t')), 'domain', t, 'name', 'my.fault');
+disturbanceSignal = quantity.Symbolic(sin(sym('t')), t, 'name', 'disturbance');
+faultSignal = quantity.Symbolic(sin(sym('t')), t, 'name', 'my.fault');
% case 1: only disturbance
u = misc.ss.combineInputSignals(myStateSpace, t.grid, 'disturbance', disturbanceSignal);
-y = quantity.Discrete(lsim(myStateSpace, u.on(), t.grid), 'domain', t, 'name', 'y');
+y = quantity.Discrete(lsim(myStateSpace, u.on(), t.grid), t, 'name', 'y');
odeResiduum = - y.diff("t") - y + 2*disturbanceSignal;
tc.verifyEqual(odeResiduum.abs.median(), 0, 'AbsTol', 5e-4);
% case 2: fault and disturbance
u2 = misc.ss.combineInputSignals(myStateSpace, t.grid, 'disturbance', disturbanceSignal, ...
'my.fault', faultSignal);
-y2 = quantity.Discrete(lsim(myStateSpace, u2.on(), t.grid), 'domain', t, 'name', 'y');
+y2 = quantity.Discrete(lsim(myStateSpace, u2.on(), t.grid), t, 'name', 'y');
odeResiduum = - y2.diff("t") - y2 + 2*disturbanceSignal + 3*faultSignal;
tc.verifyEqual(odeResiduum.abs.median(), 0, 'AbsTol', 5e-4);
end
diff --git a/+unittests/+misc/testBackstepping.m b/+unittests/+misc/testBackstepping.m
index 08033a9a0cf41af5f8d207288c7137e3be851773..5c549fb39e763cbbdbf976280ede8dd800f6e7f1 100644
--- a/+unittests/+misc/testBackstepping.m
+++ b/+unittests/+misc/testBackstepping.m
@@ -10,17 +10,17 @@ myGrid = linspace(0, 1, 31);
myDomain = [quantity.Domain('z', myGrid), quantity.Domain('zeta', myGrid)];
kernelData = quantity.Symbolic(...
[z^2/2, sin(z); 1/(2-z), -z] * magic(2) * [-zeta^2/2, sin(2*zeta); 1/(2-zeta), zeta], ...
- 'domain', myDomain, 'name', 'K');
+ myDomain, 'name', 'K');
% calculate derivatives symbolic on raw data
domainSelector = quantity.Discrete(quantity.Symbolic(zeta<=z, ...
- 'domain', myDomain, 'name', ''));
+ myDomain, 'name', ''));
kernelData_dz = domainSelector * quantity.Symbolic(...
diff([z^2/2, sin(z); 1/(2-z), -z] * magic(2) * [-zeta^2/2, sin(2*zeta); 1/(2-zeta), zeta], z), ...
- 'domain', myDomain, 'name', 'K');
+ myDomain, 'name', 'K');
kernelData_dzeta = domainSelector * quantity.Symbolic(...
diff([z^2/2, sin(z); 1/(2-z), -z] * magic(2) * [-zeta^2/2, sin(2*zeta); 1/(2-zeta), zeta], zeta), ...
- 'domain', myDomain, 'name', 'K');
+ myDomain, 'name', 'K');
% set up kernel data in misc.Backstepping objects
k = misc.Backstepping('kernelInverse', kernelData, ...
@@ -40,7 +40,7 @@ myGrid = linspace(0, 1, 31);
myDomain = [quantity.Domain('z', myGrid), quantity.Domain('zeta', myGrid)];
kernelData = quantity.Symbolic(...
[z^2/2, sin(z); 1/(2-z), -z] * magic(2) * [-zeta^2/2, sin(2*zeta); 1/(2-zeta), zeta], ...
- 'domain', myDomain, 'name', 'K');
+ myDomain, 'name', 'K');
% set up kernel data in misc.Backstepping objects
k(1) = misc.Backstepping('kernel', kernelData, ...
@@ -68,10 +68,10 @@ for it = 1 : numel(k)
end
% compare to old implementation for int_0^z
-myDomain = quantity.Discrete(quantity.Symbolic(zeta<z, 'domain', myDomain));
+myDomain = quantity.Discrete(quantity.Symbolic(zeta<z, myDomain));
Kref = quantity.Discrete(...
misc.invertBacksteppingKernel(kernelData.on({myGrid, myGrid}, {'z', 'zeta'}), -1), ...
- 'domain', [quantity.Domain('z', myGrid), quantity.Domain('zeta', myGrid)]);
+ [quantity.Domain('z', myGrid), quantity.Domain('zeta', myGrid)]);
k1Diff = myDomain * (Kref - k(1).kernelInverse);
k5Diff = myDomain * (Kref - k(5).kernel);
testCase.verifyEqual(MAX(abs(k1Diff)), 0, 'AbsTol', 1e-5);
diff --git a/+unittests/+misc/testMisc.m b/+unittests/+misc/testMisc.m
index ab189e71d150885cf2d2aaf3f06ea57bbe778f75..6151c86686d862a980206705bab936e1a7797866 100644
--- a/+unittests/+misc/testMisc.m
+++ b/+unittests/+misc/testMisc.m
@@ -4,6 +4,35 @@ function [tests] = testMisc()
tests = functiontests(localfunctions());
end
+function testFunctionArguments(tc)
+ names = misc.functionArguments(@(z, zeta, t) z*zeta+t);
+ tc.verifyEqual(names, ["z"; "zeta"; "t"]);
+end % testFunctionArguments()
+
+function testBinomial(testCase)
+
+for k = 1:3
+ for n = k:k+3
+ % test the integer order
+ b1 = nchoosek(n, k);
+ b2 = misc.binomial(n, k);
+ testCase.verifyEqual(b1, b2);
+ end
+end
+
+% test the non integer case - example from wikipedia
+testCase.verifyEqual( misc.binomial(2.5, 2), 1.875)
+testCase.verifyEqual( misc.binomial(-1, 4), (-1)^4)
+
+for k = -1:3
+ % test non integer over zero
+ testCase.verifyEqual( misc.binomial( k*0.3, 0), 1);
+ % test non integer over one
+ testCase.verifyEqual( misc.binomial( k*0.3, 1), k * 0.3, 'AbsTol', 1e-16);
+end
+
+end
+
function testUnique4cells(tc)
Moriginal = {zeros(2), 0, rand(3), [0, 1; 0, 0], "asdf", 0, 'asdf', [0, 1; 0, 0], eps, eps};
Mexpected = {zeros(2), 0, Moriginal{3}, [0, 1; 0, 0], "asdf", eps};
diff --git a/+unittests/+misc/testMultArray.m b/+unittests/+misc/testMultArray.m
index 810df6bf86f8565408a42dd742d78e1fae7dc2ab..b896e3ffb25d1b41a6e12ad489d7cb2c1db73c6c 100644
--- a/+unittests/+misc/testMultArray.m
+++ b/+unittests/+misc/testMultArray.m
@@ -121,11 +121,11 @@ end
function test3variables(testCase)
-z = linspace(0, pi, 7);
-t = linspace(0, pi, 5);
-v = 0:pi/2:3*pi/2;
-[z, t] = ndgrid(z, t);
-[v, t2] = ndgrid(v, t(1,:));
+zDom = quantity.Domain("z", linspace(0, pi, 7));
+tDom = quantity.Domain("t", linspace(0, pi, 5));
+vDom = quantity.Domain("v", 0:pi/2:3*pi/2);
+[z, t] = ndgrid(zDom.grid, tDom.grid);
+[v, t2] = ndgrid(vDom.grid, t(1,:));
syms V T Z
bb = [cos(Z), 0; 0, sin(T)];
@@ -154,12 +154,10 @@ bc = misc.multArray(b, c, 4, 3, 2);
bc = permute(bc, [ 1 2 4 3 5 ]);
testCase.verifyTrue(numeric.near(bc, bbcc));
-
-B = quantity.Discrete(b, 'size', [ 2, 2 ], 'gridName', {'z', 't'}, 'grid', {z(:,1), t(1,:)});
-C = quantity.Discrete(c, 'size', [ 2, 3 ], 'gridName', {'v', 't'}, 'grid', {v(:,1), t2(1,:)});
+B = quantity.Discrete(b, [zDom, tDom]);
+C = quantity.Discrete(c, [vDom, tDom]);
BC = B*C;
-bcon = BC.on();
testCase.verifyTrue(numeric.near(bc, BC.on()));
diff --git a/+unittests/+misc/testSolveVolterraIntegroDifferentialIVP.m b/+unittests/+misc/testSolveVolterraIntegroDifferentialIVP.m
index 7bde4c0b73965f8bb8ce38f083c7b1aedb5ae35d..dc29fa361353c1a9f4b164d72378571edf76d7ed 100644
--- a/+unittests/+misc/testSolveVolterraIntegroDifferentialIVP.m
+++ b/+unittests/+misc/testSolveVolterraIntegroDifferentialIVP.m
@@ -17,11 +17,11 @@ classdef testSolveVolterraIntegroDifferentialIVP < matlab.unittest.TestCase
methods(TestClassSetup)
function tc = initData(tc)
tc.zeta = quantity.Domain("zeta", tc.z.grid);
- tc.M = quantity.Symbolic([1+sym("z"), sin(sym("z")); 0, 2], 'domain', tc.z);
- tc.A = quantity.Symbolic([1+sym("z"), 0; -sym("z")^2, 2], 'domain', tc.z);
- tc.B = quantity.Symbolic([1+sym("z"), 0, 0.1; 0.2, 0, 2; -1, 0.5, 0.5], 'domain', tc.z);
- tc.C = quantity.Symbolic([1+sym("z"), 0, -1; 0, 2, -1], 'domain', tc.z);
- tc.D = quantity.Symbolic([1+sym("z"), sym("z")*sym("zeta"); sym("zeta"), 2], 'domain', [tc.z, tc.zeta]);
+ tc.M = quantity.Symbolic([1+sym("z"), sin(sym("z")); 0, 2], tc.z);
+ tc.A = quantity.Symbolic([1+sym("z"), 0; -sym("z")^2, 2], tc.z);
+ tc.B = quantity.Symbolic([1+sym("z"), 0, 0.1; 0.2, 0, 2; -1, 0.5, 0.5], tc.z);
+ 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);
end
diff --git a/+unittests/+numeric/femMatricesTest.m b/+unittests/+numeric/femMatricesTest.m
index 7ac8fe8232f3f0c4c805a25ee5de2b4384f6bd77..447e93b279720b6f6eb608229aa0424fc220cff7 100644
--- a/+unittests/+numeric/femMatricesTest.m
+++ b/+unittests/+numeric/femMatricesTest.m
@@ -3,11 +3,46 @@ tests = functiontests(localfunctions);
end
function mTest(testCase)
-[M] = numeric.femMatrices('grid', linspace(0,1,3));
-M_ = [1/6, 1/12 0; ...
- 1/12, 1/3, 1/12; ...
- 0, 1/12, 1/6];
+spatialDomain = quantity.EquidistantDomain("z", 0, 1, 'stepNumber', 3);
+I = quantity.Discrete.ones(1, spatialDomain);
+[M, D2, L, P] = numeric.femMatrices(spatialDomain, ...
+ 'alpha', quantity.Discrete.ones(1, spatialDomain), ...
+ 'b', I, 'c', I);
+
+[~, D1] = numeric.femMatrices(spatialDomain, ...
+ 'beta', quantity.Discrete.ones(1, spatialDomain));
+
+[~, D0] = numeric.femMatrices(spatialDomain, ...
+ 'gamma', quantity.Discrete.ones(1, spatialDomain));
+
+
+
+M_ = spatialDomain.stepSize / 6 * ...
+ [2, 1, 0; ...
+ 1, 4, 1; ...
+ 0, 1, 2];
+
+D2_ = 1 / spatialDomain.stepSize * ...
+ [-1 1 0; ...
+ 1 -2 1; ...
+ 0 1 -1];
+
+D1_ = 1 / 2 * ...
+ [-1 1 0; ...
+ -1 0 1; ...
+ 0 -1 1];
+
+D0_ = spatialDomain.stepSize / 6 * ...
+ [2, 1, 0; ...
+ 1, 4, 1; ...
+ 0, 1, 2];
+
testCase.verifyEqual( M, M_, 'AbsTol', 10*eps );
+testCase.verifyEqual( D2, D2_, 'AbsTol', 10*eps);
+testCase.verifyEqual( D1, D1_, 'AbsTol', 10*eps);
+testCase.verifyEqual( D0, D0_, 'AbsTol', 5e-5);
+testCase.verifyEqual( L, M_ * I.on(), 'AbsTol', 10*eps);
+testCase.verifyEqual( P, I.on()' * M_, 'AbsTol', 10*eps);
end
\ No newline at end of file
diff --git a/+unittests/+quantity/testBasicVariable.m b/+unittests/+quantity/testBasicVariable.m
deleted file mode 100644
index daeca9c876dfcac1170fd05d59ee63556dd164e7..0000000000000000000000000000000000000000
--- a/+unittests/+quantity/testBasicVariable.m
+++ /dev/null
@@ -1,25 +0,0 @@
-function [tests] = testGevrey()
-%TESTGEVREY Summary of this function goes here
-% Detailed explanation goes here
-tests = functiontests(localfunctions);
-end
-
-function testBasicVariableInit(testCase)
-
-%%
-t = linspace(0, 1, 11)';
-g1 = signals.GevreyFunction('order', 1.8, 't', t);
-g2 = signals.GevreyFunction('order', 1.5, 't', t);
-
-G = [g1; g2];
-N_diff = 2;
-derivatives = G.diff(0:N_diff);
-
-GD = zeros(length(t), N_diff+1, 2);
-for k = 1:N_diff+1
- GD(:, k, 1) = signals.gevrey.bump.g(t, k-1, t(end), g1.sigma);
- GD(:, k, 2) = signals.gevrey.bump.g(t, k-1, t(end), g2.sigma);
-end
-
-testCase.verifyEqual(derivatives.on(), GD);
-end
\ No newline at end of file
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index becc80275197e369379b46853205d52b5d5af274..29945255d0d0da590bf163f05230491119109922 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -4,6 +4,56 @@ function [tests] = testDiscrete()
tests = functiontests(localfunctions);
end
+function testEigScalar(tc)
+z = quantity.Domain("z", linspace(-1, 1, 11));
+[Vz, Dz, Wz] = z.Discrete.eig();
+tc.verifyEqual(z.grid, Dz.on());
+tc.verifyEqual(ones(z.n, 1), Vz.on());
+tc.verifyEqual(ones(z.n, 1), Wz.on());
+end % testEigScalar()
+
+function testEigMat(tc)
+z = quantity.Domain("z", linspace(-1, 1, 11));
+zSym = sym("z");
+quan = quantity.Discrete(quantity.Symbolic([1+zSym^2, 2; -zSym, zSym], z));
+% all output arguments
+[Vq, Dq, Wq] = quan.eig();
+tc.verifyEqual(diag(eig(quan.at(0.4))), Dq.at(0.4), "AbsTol", 1e-15);
+tc.verifyEqual(on(quan * Vq), on(Vq * Dq), "AbsTol", 1e-15);
+tc.verifyEqual(on(Wq' * quan), on(Dq * Wq'), "AbsTol", 1e-15);
+% 2 output arguments
+[Vq, Dq] = quan.eig();
+tc.verifyEqual(diag(eig(quan.at(0.4))), Dq.at(0.4), "AbsTol", 1e-15);
+tc.verifyEqual(on(quan * Vq), on(Vq * Dq), "AbsTol", 1e-15);
+% 1 output argument
+Eq = quan.eig();
+tc.verifyEqual(eig(quan.at(0.4)), Eq.at(0.4), "AbsTol", 1e-15);
+tc.verifyEqual(on(Eq), on(Dq.diag2vec()), "AbsTol", 1e-15);
+end % testEigMat()
+
+function testPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Discrete(sin(z.grid * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
+AA = A.^2;
+
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A.^0), ones(z.n, 1));
+tc.verifyEqual( on(A.^3), on(AA * A), 'AbsTol', 1e-15);
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Discrete + zeta.Discrete, A + 0*zeta.Discrete; ...
+ z.Discrete+zeta.Discrete, A * subs(A, "z", "zeta")];
+BBB = B.^3;
+BBB_alt = B .* B .* B;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6})).^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B.^0), ones([B(1).domain(1).n, B(1).domain(2).n, 2, 2]), ...
+ 'AbsTol', 1e-15);
+end % testPower()
+
function setupOnce(testCase)
syms z zeta t sy
testCase.TestData.sym.z = z;
@@ -12,6 +62,64 @@ testCase.TestData.sym.t = t;
testCase.TestData.sym.sy = sy;
end
+function testNorm(tc)
+z = quantity.Domain("z", linspace(0, 1, 11));
+quan = Discrete(z) * [-1; 2; 3];
+tc.verifyEqual(quan.norm.on(), on(z.Discrete*sqrt(1+4+9)), 'AbsTol', 1e-15);
+tc.verifyEqual(quan.norm.on(), on(quan.norm(2)), 'AbsTol', 1e-15);
+end
+
+function testLog(tc)
+z = quantity.Domain("z", linspace(0, 1, 11));
+quan = Discrete(z) * ones(2, 2) + [1, 1; 0, 0]+0.1*ones(2,2);
+
+tc.verifyEqual(log(quan.on()), on(log(quan)), 'AbsTol', 1e-15);
+tc.verifyEqual(on(quan), on(log(exp(quan))), 'AbsTol', 1e-15);
+end % testLog()
+
+function testLog10(tc)
+z = quantity.Domain("z", linspace(0, 1, 11));
+quan = Discrete(z) * ones(2, 2) + [1, 1; 0, 0]+0.1*ones(2,2);
+
+tc.verifyEqual(log10(quan.on()), on(log10(quan)), 'AbsTol', 1e-15);
+end % testLog10()
+
+function powerTest(testCase)
+t = quantity.Domain("t", linspace(0,1));
+
+p(1,1) = Discrete(t);
+p(2,1) = Discrete(t)^2;
+
+p2 = p.^2;
+
+testCase.verifyEqual(p2(1).on(), t.grid.^2)
+testCase.verifyEqual(p2(2).on(), t.grid.^4, 'AbsTol', 1e-15)
+
+p25 = p.^2.5;
+
+testCase.verifyEqual(p25(1).on(), t.grid.^2.5);
+testCase.verifyEqual(p25(2).on(), t.grid.^5, 'AbsTol', 1e-15);
+
+end
+
+function timesTest(testCase)
+
+ t = quantity.Domain.defaultDomain(13, "t");
+
+ a(1) = quantity.Discrete(t.grid, t);
+ a(2) = quantity.Discrete(t.grid.^2, t);
+ a(3) = quantity.Discrete(t.grid.^3, t);
+
+ p = a .* a;
+
+ p_(1) = quantity.Discrete(t.grid.^2, t);
+ p_(2) = quantity.Discrete(t.grid.^4, t);
+ p_(3) = quantity.Discrete(t.grid.^6, t);
+
+ testCase.verifyEqual( p.on(), p_.on(),'AbsTol', 1e-15)
+
+end
+
function testFindZeros(testCase)
x = quantity.Domain('x', linspace(0,1));
@@ -480,14 +588,14 @@ end
function testSqrt(testCase)
% 1 spatial variable
quanScalar1d = quantity.Discrete((linspace(0,2).^2).', quantity.Domain("z", linspace(0, 1)), ...
- 'size', [1, 1], 'name', "s1d");
+ 'name', "s1d");
quanScalarRoot = quanScalar1d.sqrt();
testCase.verifyEqual(quanScalarRoot.on(), linspace(0, 2).');
% 2 spatial variables
quanScalar2d = quantity.Discrete((linspace(0, 2, 21).' + linspace(0, 2, 41)).^2, ...
[quantity.Domain("z", linspace(0, 1, 21)), quantity.Domain("zeta", linspace(0, 1, 41))], ...
- 'size', [1, 1], 'name', 's2d');
+ 'name', 's2d');
quanScalar2dRoot = quanScalar2d.sqrt();
testCase.verifyEqual(quanScalar2dRoot.on(), (linspace(0, 2, 21).' + linspace(0, 2, 41)));
@@ -504,10 +612,10 @@ end
function testSqrtm(testCase)
quanScalar1d = quantity.Discrete((linspace(0,2).^2).', quantity.Domain("z", linspace(0, 1)), ...
- 'size', [1, 1], 'name', 's1d');
+ 'name', 's1d');
quanScalar2d = quantity.Discrete((linspace(0, 2, 21).' + linspace(0, 2, 41)).^2, ...
[quantity.Domain("z", linspace(0, 1, 21)), quantity.Domain("zeta", linspace(0, 1, 41))], ...
- 'size', [1, 1], 'name', 's2d');
+ 'name', 's2d');
% diagonal matrix - 1 spatial variable
quanDiag = [quanScalar1d, 0*quanScalar1d; 0*quanScalar1d, (2)*quanScalar1d];
quanDiagRoot = quanDiag.sqrt(); % only works for diagonal matrices
@@ -593,11 +701,16 @@ end
function testSolveDVariableEqualQuantityConstant(testCase)
%% simple constant case
-quan = quantity.Discrete(ones(3, 1), quantity.Domain("z", linspace(0, 1, 3)), ...
- 'size', [1, 1], 'name', 'a');
+quan = quantity.Discrete(ones(3, 1), quantity.Domain("z", linspace(0, 1, 3)));
solution = quan.solveDVariableEqualQuantity();
[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 3), linspace(0, 1, 3));
testCase.verifyEqual(solution.on(), referenceResult1 + referenceResult2, 'AbsTol', 10*eps);
+% %% 2x2 constant case
+%Lambda = quantity.Discrete(quantity.Symbolic([1+sym("z")^2; sym("z")-2], quantity.Domain("z", linspace(0, 1, 101))));
+%Lambda = quantity.Discrete(quantity.Symbolic([2; -1], quantity.Domain("z", linspace(0, 1, 11))));
+%solution = Lambda.solveDVariableEqualQuantity();
+%[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 3), linspace(0, 1, 3));
+%testCase.verifyEqual(solution.on(), referenceResult1 + referenceResult2, 'AbsTol', 10*eps);
end
function testSolveDVariableEqualQuantityNegative(testCase)
@@ -624,8 +737,8 @@ function testSolveDVariableEqualQuantityComparedToSym(testCase)
z = testCase.TestData.sym.z;
Z = quantity.Domain("z", linspace(0, 1, 5));
assume(z>0 & z<1);
-quanBSym = quantity.Symbolic(1+z, Z, 'name', 'bSym', 'gridName', {'z'});
-quanBDiscrete = quantity.Discrete(quanBSym.on(), Z, 'name', 'bDiscrete', 'size', size(quanBSym));
+quanBSym = quantity.Symbolic(1+z, Z);
+quanBDiscrete = quantity.Discrete(quanBSym.on(), Z);
solutionBSym = quanBSym.solveDVariableEqualQuantity();
solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
%solutionBSym.plot(); solutionBDiscrete.plot();
@@ -655,9 +768,9 @@ end
function testMtimesDifferentGridLength(testCase)
%%
a = quantity.Discrete(ones(11, 1), quantity.Domain("z", linspace(0, 1, 11)), ...
- 'size', [1, 1], 'name', 'a');
+ 'name', 'a');
b = quantity.Discrete(ones(5, 1), quantity.Domain("z", linspace(0, 1, 5)), ...
- 'size', [1, 1], 'name', 'b');
+ 'name', 'b');
ab = a*b;
%
z = testCase.TestData.sym.z;
@@ -907,7 +1020,6 @@ vInv = v.inv();
%% matrix example
t = quantity.Domain("t", linspace(0, pi));
-
[T, ~] = ndgrid(t.grid, s.grid);
K = quantity.Discrete(permute(repmat(2*eye(2), [1, 1, size(T)]), [3, 4, 1, 2]), [t, s], ...
@@ -1116,7 +1228,7 @@ bc = permute(bc, [ 1 2 4 3 5 ]);
B = quantity.Discrete(b, [Z, T], 'size', [ 2, 2 ]);
C = quantity.Discrete(c, [V, T], 'size', [2, 3]);
BC = B*C;
-testCase.verifyTrue(numeric.near(bc, BC.on));
+testCase.verifyTrue(numeric.near(bc, BC.on()));
%%
z = linspace(0,1).';
@@ -1148,13 +1260,13 @@ function testOnConstant(testCase)
A = [0 1; 2 3];
const = quantity.Discrete(A, quantity.Domain.empty);
-C = const.on([1:5]);
+C = const.on(1:1:5);
for k = 1:numel(A)
testCase.verifyTrue( all( A(k) == C(:,k) ));
end
-end
+end % testOnConstant()
function testIsConstant(testCase)
@@ -1191,21 +1303,28 @@ verifyTrue(testCase, numeric.near(CA.on(), ca, 1e-12));
end
-function testMPower(testCase)
-
-%%
-z = linspace(0,1).';
-a = sin(z * pi);
-
-A = quantity.Discrete({a}, quantity.Domain("z", z));
-
-aa = sin(z * pi) .* sin(z * pi);
+function testMPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Discrete(sin(z.grid * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
AA = A^2;
-verifyTrue(testCase, numeric.near(aa, AA.on()));
-
-
-end
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A^0), ones(z.n, 1));
+tc.verifyEqual( on(A^3), on(AA * A));
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Discrete + zeta.Discrete, A + 0*zeta.Discrete; ...
+ z.Discrete+zeta.Discrete, A * subs(A, "z", "zeta")];
+BBB = B^3;
+BBB_alt = B * B * B;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6}))^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B^0), permute(repmat(eye(2), [1, 1, B(1).domain(1).n, B(1).domain(2).n]), [3, 4, 1, 2]), ...
+ 'AbsTol', 1e-15);
+end % testMPower()
function testPlus(testCase)
%%
@@ -1223,9 +1342,9 @@ bzZeta = cos(zGrid * pi) + cos(zetaGrid * pi);
Z = quantity.Domain("z", z);
Zeta = quantity.Domain("zeta", zeta);
-AB = quantity.Discrete({a, b}, Z.copyAndRename("blub"));
+AB = quantity.Discrete({a, b}, Z.rename("blub"));
A = quantity.Discrete({a}, Z);
-ALr = quantity.Discrete({aLr}, Zeta.copyAndRename("z"));
+ALr = quantity.Discrete({aLr}, Zeta.rename("z"));
B = quantity.Discrete({b}, Z);
C = quantity.Discrete({C}, Zeta);
AZZETA = quantity.Discrete({azZeta}, [Z, Zeta]);
@@ -1292,7 +1411,7 @@ testCase.verifyEqual( on( TZX + XTZ ), on( 2 * TZX ), 'AbsTol', eps )
end
-function testInit(testCase)
+function testInit(tc)
%%
z = linspace(0,1).';
t = linspace(0,1,101);
@@ -1306,20 +1425,18 @@ a = quantity.Discrete(v, [Z, T]);
b = quantity.Discrete(V, [Z, T]);
%
-verifyTrue(testCase, misc.alln(a.on() == b.on()));
+verifyTrue(tc, all(a.on() == b.on(), 'all'));
%%
-z = linspace(0,1,31).';
-zeta = linspace(0,1,51);
-
+zeta = T.grid;
[zGrid, zetaGrid] = ndgrid(z, zeta);
-
Z = quantity.Domain("z", z);
Zeta = quantity.Domain("zeta", zeta);
-c = quantity.Discrete(zGrid, [Z, Zeta], 'name', 'a');
-d = quantity.Discrete(zetaGrid', [Z, Zeta], 'name', 'b');
-%TODO write test case, the last line should not work
+c = quantity.Discrete(zGrid, [Z, Zeta], 'name', 'a'); % this should work
+tc.verifyError(...
+ @() quantity.Discrete(zetaGrid.', [Z, Zeta], 'name', 'd'), ...% this should not work
+ '');
end
@@ -1369,7 +1486,7 @@ v = rand([100, 42, 2, 3]);
Z = quantity.Domain('z', linspace(0,1,100));
T = quantity.Domain('t', linspace(0,1,42));
V = quantity.Discrete( quantity.Discrete.value2cell(v, [100, 42], [2, 3]), [Z, T]);
-verifyTrue(testCase, misc.alln(v == V.on()));
+verifyTrue(testCase, all(v == V.on(), 'all'));
end
function testMtimesDifferentGrid(testCase)
@@ -1508,6 +1625,13 @@ quan4dArbitrary = quan4d.subs({'z', 'zeta', 'eta', 'beta'}, {'zeta', 'beta', 'z'
testCase.verifyEqual(reshape(quan4d.on({1, 1, 1, 1}), size(quan4d)), ...
reshape(quan4dArbitrary.on({1, 1, 1}), size(quan4dAllEta)));
+% test a substution where a quantity.EquidistantDomain is used:
+e = quantity.EquidistantDomain("t", 0, 1, "stepSize", 0.1);
+q = Discrete(e);
+p = q.subs(e, "z");
+testCase.verifyEqual( p.domain.name, "z")
+testCase.verifyEqual( on( q.subs(e, z.grid) ), q.on(z.grid) )
+
end
diff --git a/+unittests/+quantity/testDomain.m b/+unittests/+quantity/testDomain.m
index 11ee1963e5db1e7dbf3a4d522e3f796a5a34a554..197e367c4f5bed581279357e15f2590ad1487bf7 100644
--- a/+unittests/+quantity/testDomain.m
+++ b/+unittests/+quantity/testDomain.m
@@ -18,9 +18,34 @@ tc.TestData.d = d;
end
-function testParser(testCase)
-
-%profile on
+function testDiscrete(tc)
+d = quantity.Domain("a", linspace(0, 1, 3));
+tc.verifyEqual(d.Discrete.on(), d.grid);
+end % testDiscrete
+
+function testSymbolic(tc)
+d = quantity.Domain("a", linspace(0, 1, 3));
+tc.verifyEqual(d.Symbolic.on(), d.grid);
+end % testDiscrete
+
+function testRename(tc)
+d = [quantity.Domain("a", linspace(0, 1, 3)), ...
+ quantity.Domain("b", linspace(0, 2, 5)), ...
+ quantity.Domain("ab", linspace(0, pi, 4))];
+xyz = d.rename(["x", "y", "z"]);
+xyz2 = d.rename(["z", "y", "x"], ["ab", "b", "a"]);
+tc.verifyEqual([xyz.name], ["x", "y", "z"]);
+tc.verifyEqual([xyz.name], [xyz2.name]);
+
+cba = d.rename(["c", "b", "a"], ["a", "b", "ab"]);
+tc.verifyEqual([cba.name], ["c", "b", "a"]);
+tc.verifyEqual([cba.n], [d.n]);
+
+zZetaDomain = [quantity.Domain("z", linspace(0, 1, 3)), quantity.Domain("zeta", linspace(0, 1, 3))];
+tc.verifyEqual([zZetaDomain.rename("a", "z").name], ["a", "zeta"])
+end % testRename()
+
+function testParser(tc)%profile on
d = quantity.Domain.parser('blub', 1, 'grid', 1:3, 'gridName', 'test');
a = quantity.Domain.parser('domain', d);
b = quantity.Domain.parser('blabla', 0, 'domain', a);
@@ -29,47 +54,46 @@ b = quantity.Domain.parser('blabla', 0, 'domain', a, 'blub', 'a');
%profile viewer
-testCase.verifyEqual(d, a);
-testCase.verifyEqual(b, a);
+tc.verifyEqual(d, a);
+tc.verifyEqual(b, a);
end
-function testSplit(testCase)
+function testSplit(tc)
a = quantity.Domain('a', -5:15 );
b = a.split([-1, 0, 4]);
-testCase.verifyEqual(b(1).grid', -5:-1)
-testCase.verifyEqual(b(2).grid', -1:0)
-testCase.verifyEqual(b(3).grid', 0:4)
-testCase.verifyEqual(b(4).grid', 4:15)
+tc.verifyEqual(b(1).grid', -5:-1)
+tc.verifyEqual(b(2).grid', -1:0)
+tc.verifyEqual(b(3).grid', 0:4)
+tc.verifyEqual(b(4).grid', 4:15)
-testCase.verifyWarning(@() a.split([0.1, 2.6]), 'DOMAIN:split')
-testCase.verifyWarningFree(@() a.split([0.1, 2.6], 'AbsTol', 0.4), 'DOMAIN:split')
+tc.verifyWarning(@() a.split([0.1, 2.6]), 'DOMAIN:split')
+tc.verifyWarningFree(@() a.split([0.1, 2.6], 'AbsTol', 0.4), 'DOMAIN:split')
c = a.split([0.1, 2.6], 'warning', false);
-testCase.verifyEqual(c(1).grid', -5:0)
-testCase.verifyEqual(c(2).grid', 0:3)
-testCase.verifyEqual(c(3).grid', 3:15)
+tc.verifyEqual(c(1).grid', -5:0)
+tc.verifyEqual(c(2).grid', 0:3)
+tc.verifyEqual(c(3).grid', 3:15)
end
-function testContains(testCase)
+function testContains(tc)
a = quantity.Domain('a', 0:10);
b = quantity.Domain('b', 3:5);
-testCase.verifyTrue( a.contains(b) )
-testCase.verifyTrue( a.contains(1) )
-testCase.verifyFalse( a.contains(11) )
-testCase.verifyFalse( a.contains(1.5) )
+tc.verifyTrue( a.contains(b) )
+tc.verifyTrue( a.contains(1) )
+tc.verifyFalse( a.contains(11) )
+tc.verifyFalse( a.contains(1.5) )
A = [quantity.Domain('a1', 0:10), quantity.Domain('a2', 0:10)];
-testCase.verifyTrue( A.contains(5) )
-
+tc.verifyTrue( A.contains(5) )
end
@@ -98,6 +122,15 @@ tc.verifyEqual( d.find([a.name b.name]), [a b])
end
+function testRemove(tc)
+d = tc.TestData.d;
+dWithoutA = d.remove("a");
+tc.verifyEqual([dWithoutA.name], ["b" "c"]);
+
+a = d.remove(dWithoutA);
+tc.verifyEqual(a.name, "a");
+end
+
function testGet(tc)
z = linspace(0, pi, 3);
diff --git a/+unittests/+quantity/testEquidistantDomain.m b/+unittests/+quantity/testEquidistantDomain.m
index cb0b05ee38a4aae9cd926b904aaf25938a8e84eb..981f819a43ce6a288a1b7aad998889d01d3afa6f 100644
--- a/+unittests/+quantity/testEquidistantDomain.m
+++ b/+unittests/+quantity/testEquidistantDomain.m
@@ -40,4 +40,17 @@ EE(1:2) = quantity.Domain();
EE(1) = e;
EE(2) = d;
+end
+
+function testSubs(testCase)
+
+f = quantity.Discrete((1:5)', quantity.EquidistantDomain("t", 0, 1, "stepNumber", 5));
+testCase.verifyEqual( f.subs("t", 1), f.on(1))
+testCase.verifyEqual( f.changeGrid(quantity.Domain("t", linspace(0,1))).on(), linspace(1,5)', 'AbsTol', 1e-15);
+
+f = quantity.Symbolic(sym("t"), quantity.EquidistantDomain("t", 0, 1, "stepNumber", 5));
+testCase.verifyEqual( f.subs("t", 1), f.on(1))
+testCase.verifyEqual( f.changeGrid(quantity.Domain("t", linspace(0,1))).on(), linspace(0,1)', 'AbsTol', 1e-15);
+
+
end
\ No newline at end of file
diff --git a/+unittests/+quantity/testFunction.m b/+unittests/+quantity/testFunction.m
index c476cf957282700ae272f3ea18e0c3b373f7d073..966ae75b9b5dfa0fdd9e2c8dbb716820002a1b4c 100644
--- a/+unittests/+quantity/testFunction.m
+++ b/+unittests/+quantity/testFunction.m
@@ -4,6 +4,17 @@ function [tests] = testFunction()
tests = functiontests(localfunctions);
end
+function testAt(tc)
+z = quantity.Domain("z", linspace(0, 1, 11));
+f = quantity.Function({@(z)sin(z), @(z)1+z; @(z)cos(z), @(z)z.^2}, z);
+fZZeta = quantity.Function({@(z,zeta)sin(zeta)+z, @(z,zeta)1+zeta-z; @(z,zeta)cos(z), @(z,zeta)zeta.^2}, ...
+ [z, z.rename("zeta")]);
+fDisc = quantity.Discrete(f);
+tc.verifyLessThan(abs(f.at({0.2}) - fDisc.at({0.2})), 10*eps);
+fZZetaDisc = quantity.Discrete(fZZeta);
+tc.verifyLessThan(abs(fZZeta.at({0.2, 0.4}) - fZZetaDisc.at({0.2, 0.4})), 10*eps);
+end % testAt()
+
function testInt(testCase)
z = quantity.Domain("z", linspace(0,1));
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index f5d7ca4f8f812e85c7381b4651d254931bc236b9..c378c5dc5bd39a892feddaa5d84d9ce75a3b2c02 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -4,6 +4,99 @@ function [tests ] = testSymbolic()
tests = functiontests(localfunctions());
end
+function testValueContinuousBug(tc)
+z = quantity.Domain("z", linspace(0, 1, 7));
+zeta = quantity.Domain("zeta", linspace(0, 1, 5));
+A = quantity.Symbolic(sin(sym("z") * pi), z);
+B1 = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
+ z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
+B2 = [zeta.Symbolic + 0*z.Symbolic, A + 0*zeta.Symbolic; ...
+ z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
+tc.verifyEqual(B1.on([z, zeta]), B2.on([z, zeta]));
+
+C1 = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic];
+c1 = zeta.Symbolic + 0*z.Symbolic;
+c2 = A + 0*zeta.Symbolic;
+C2 = [c1, c2];
+tc.verifyEqual(C1.on([z, zeta]), C2.on([z, zeta]));
+
+% TODO: fix this error! See issue #41
+% this problem is maybe caused as input arguments of B2(1).valueContinuous is somehow ordered
+% differently compared to B2(2:4).valueContinuous or all B(:).valueContinuous,
+% i.e. B2(1).valueContinuous is (zeta, z) instead of (z, zeta).
+end % testValueContinuousBug()
+
+function testPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Symbolic(sin(sym("z") * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
+AA = A.^2;
+
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A.^0), ones(z.n, 1));
+tc.verifyEqual( on(A.^3), on(AA * A), 'AbsTol', 1e-15);
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
+ z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
+BBB = B.^3;
+BBB_alt = B .* B .* B;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6})).^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B.^0), ones([B(1).domain(1).n, B(1).domain(2).n, 2, 2]), ...
+ 'AbsTol', 1e-15);
+end % testPower()
+
+function testMPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Symbolic(sin(sym("z") * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
+AA = A^2;
+
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A^0), ones(z.n, 1));
+tc.verifyEqual( on(A^3), on(AA * A));
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
+ z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
+B_discrete = quantity.Discrete(B);
+BBB = B^3;
+BBB_alt = B * B * B;
+BBB_discrete = B_discrete^3;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6}))^3, 'AbsTol', 1e-15);
+tc.verifyEqual(BBB_discrete.at({0.5, 0.6}), (B.at({0.5, 0.6}))^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B^0), permute(repmat(eye(2), [1, 1, B(1).domain(1).n, B(1).domain(2).n]), [3, 4, 1, 2]), ...
+ 'AbsTol', 1e-15);
+end % testMPower()
+
+function testNorm(tc)
+z = quantity.Domain("z", linspace(0, 1, 11));
+quan = Symbolic(z) * [-1; 2; 3];
+tc.verifyEqual(quan.norm.on(), on(z.Discrete*sqrt(1+4+9)), 'AbsTol', 1e-15);
+tc.verifyEqual(quan.norm.on(), on(quan.norm(2)), 'AbsTol', 1e-15);
+end % testNorm()
+
+function testLog(tc)
+z = quantity.Domain("z", linspace(0, 1, 11));
+quan = Symbolic(z) * ones(2, 2) + [1, 1; 0, 0]+0.1*ones(2,2);
+
+tc.verifyEqual(log(quan.on()), on(log(quan)), 'AbsTol', 1e-15);
+tc.verifyEqual(on(quan), on(log(exp(quan))), 'AbsTol', 1e-15);
+end % testLog()
+
+function testLog10(tc)
+z = quantity.Domain("z", linspace(0, 1, 11));
+quan = Symbolic(z) * ones(2, 2) + [1, 1; 0, 0]+0.1*ones(2,2);
+
+tc.verifyEqual(log10(quan.on()), on(log10(quan)), 'AbsTol', 1e-15);
+end % testLog10()
+
function testDet(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
zSymbolic = quantity.Symbolic(sym("z"), z);
@@ -230,10 +323,10 @@ syms z zeta eta e a
myGrid = linspace(0, 1, 7);
myDomain(1) = quantity.Domain("z", myGrid);
-myDomain(2) = myDomain(1).copyAndRename("zeta");
-myDomain(3) = myDomain(1).copyAndRename("eta");
-myDomain(4) = myDomain(1).copyAndRename("e");
-myDomain(5) = myDomain(1).copyAndRename("a");
+myDomain(2) = myDomain(1).rename("zeta");
+myDomain(3) = myDomain(1).rename("eta");
+myDomain(4) = myDomain(1).rename("e");
+myDomain(5) = myDomain(1).rename("a");
obj = quantity.Symbolic([z*zeta, eta*z; e*a, a], myDomain);
@@ -470,7 +563,7 @@ function testZero(testCase)
x = quantity.Symbolic(0, quantity.Domain.empty());
verifyEqual(testCase, x.on(), 0);
-verifyTrue(testCase, misc.alln(on(x) * magic(5) == zeros(5)));
+verifyTrue(testCase, all(on(x) * magic(5) == zeros(5), 'all'));
end
@@ -494,7 +587,7 @@ d = [quantity.Domain('z', linspace(0, 1, 5)), ...
x = quantity.Symbolic(sin(z * t * pi), d);
o = x / x;
-verifyTrue(testCase, misc.alln(o.on == 1));
+verifyTrue(testCase, all(o.on() == 1, 'all'));
% test constant values
invsin = 1 / ( x + 8);
@@ -627,7 +720,7 @@ b = quantity.Symbolic(v, myDom2);
c = quantity.Symbolic(v, myDom2);
%%
-verifyTrue(testCase, misc.alln(b.on() == c.on()));
+verifyTrue(testCase, all(b.on() == c.on(), 'all'));
end
function testPlus(testCase)
@@ -868,13 +961,23 @@ solution = f.solveAlgebraic([1, 1], f(1).gridName{1});
testCase.verifyEqual(solution, 0.5);
end
+function testSubs3(tc)
+ zZeta = [quantity.Domain("z", linspace(0, 1, 11)), quantity.Domain("zeta", linspace(0, 1, 11))];
+ Id = quantity.Symbolic([sym("z"), 1; 2 sym("zeta")], zZeta);
+ Id_z1 = quantity.Symbolic([1, 1; 2 sym("zeta")], zZeta.find("zeta"));
+ Id_zetaZ = quantity.Symbolic( [sym("zeta"), 1; 2 sym("z")], zZeta([2 1]));
+
+ tc.verifyEqual(Id.subs("z", 1).on(), Id_z1.on(), 'AbsTol', 10*eps);
+ tc.verifyEqual( Id.subs({"z", "zeta"}, {"zeta", "z"}).on(), Id_zetaZ.on())
+end % testSubs3()
+
function testSubs2(tc)
z = quantity.Domain("z", linspace(0, 1, 11));
Id = quantity.Symbolic(eye(2), z);
tc.verifyEqual(Id.on(), Id.subs("z", "zeta").on(), 'AbsTol', 10*eps);
end % testSubs2()
-function testSubs(testCase)
+function testSubs(tc)
%%
% init
syms x y z w
@@ -901,9 +1004,9 @@ Z0 = quantity.Domain('z', 0);
fOfz = f.subs([X, Y], [Z Z]);
fOfyx = f.subs({'x', 'y'}, {'w', 'x'});
-testCase.verifyEqual(Ff1, fOf1);
-testCase.verifyEqual(Ffz.sym(), fOfz.sym());
-testCase.verifyEqual(Fyx.sym(), fOfyx.sym());
+tc.verifyEqual(Ff1, fOf1);
+tc.verifyEqual(Ffz.sym(), fOfz.sym());
+tc.verifyEqual(Fyx.sym(), fOfyx.sym());
%%
syms z zeta
@@ -911,28 +1014,15 @@ assume(z>0 & z<1); assume(zeta>0 & zeta<1);
F = quantity.Symbolic(zeros(1), ...
[quantity.Domain.defaultDomain(3, "z"), quantity.Domain.defaultDomain(3, "zeta")]);
Feta = F.subs('z', 'eta');
-testCase.verifyEqual(Feta(1).gridName, ["eta", "zeta"]);
+tc.verifyEqual(Feta(1).gridName, ["eta", "zeta"]);
%%
fMessy = quantity.Symbolic(x^1*y^2*z^4*w^5, [X, Y, Z, W]);
fMessyA = fMessy.subs({'x', 'y', 'z', 'w'}, {'a', 'a', 'a', 'a'});
fMessyYY = fMessy.subs({'w', 'x', 'z'}, {'xi', 'y', 'y'});
-testCase.verifyEqual(fMessyA.gridName, "a");
-testCase.verifyEqual(fMessyA.grid, {linspace(0, 1, 11)'});
-testCase.verifyEqual(fMessyYY.gridName, ["y", "xi"]);
-testCase.verifyEqual(fMessyYY.grid, {linspace(0, 1, 11)', linspace(0, 1, 3)'});
-% % sub multiple numerics -> not implemented yet
-% f11 = f.subs('x', [1; 2]);
-% f1 = f.subs('x', 1);
-% f2 = f.subs('x', 2);
-% testCase.verifyEqual(reshape(f11(1,:,:).sym(), size(f)), fOf1);
-% testCase.verifyEqual(reshape(f11(2,:,:).sym(), size(f)), f2.sym());
-%
-% % sub multiple multiple variables
-% f1232 = f.subs({'x', 'y'}, {[1; 2], [3; 2]});
-% f13 = f.subs({'x', 'y'}, {1, 3});
-% f22 = f.subs({'x', 'y'}, {2, 2});
-% testCase.verifyEqual(reshape(f1232(1,:,:), size(f)), f22);
-% testCase.verifyEqual(reshape(f1232(2,:,:), size(f)), f13);
-end
+tc.verifyEqual(fMessyA.gridName, "a");
+tc.verifyEqual(fMessyA.grid, {linspace(0, 1, 11)'});
+tc.verifyEqual(fMessyYY.gridName, ["y", "xi"]);
+tc.verifyEqual(fMessyYY.grid, {linspace(0, 1, 11)', linspace(0, 1, 3)'});
+end % testSubs()
diff --git a/+unittests/+signals/signalsTest.m b/+unittests/+signals/signalsTest.m
index 7fbced81d578400c443fb413aa5808a168ee142b..6980acce85b50c7b60eafcf1d046b302ef2effab 100644
--- a/+unittests/+signals/signalsTest.m
+++ b/+unittests/+signals/signalsTest.m
@@ -21,10 +21,10 @@ d = 1;
sig = 3 / d^2 * z^2 - 2 / d^3 * z^3;
Z = quantity.Domain('z', linspace(0, 1, 101));
-SI = quantity.Symbolic(sig, 'domain', Z);
+SI = quantity.Symbolic(sig, Z);
% compute step by function
-si = signals.smoothStep(1, 'domain', Z);
+si = quantity.Symbolic( misc.ss.setPointChange(0, 1, 0, 1, 1, "var", Z.name), Z );
%
testCase.verifyEqual(si.on(), SI.on(), 'AbsTol', 6e-16);
@@ -36,14 +36,212 @@ function testGevreyBump(testCase)
% signals.
g = signals.GevreyFunction('diffShift', 1);
-b = signals.gevrey.bump.dgdt_1(g.T, g.sigma, g.grid{1}) * 1 / signals.gevrey.bump.dgdt_0(g.T, g.sigma, g.T);
+b = signals.gevrey.bump.dgdt_1(g.T, g.sigma, g.domain.grid) * 1 / signals.gevrey.bump.dgdt_0(g.T, g.sigma, g.T);
%% Test scaling of Gevrey function
-verifyEqual(testCase, b, g.on(), 'AbsTol', 1e-13);
+verifyEqual(testCase, b, g.fun.on(), 'AbsTol', 1e-13);
% test the derivatives of the bump
g = signals.GevreyFunction('diffShift', 0);
-g.diff(0:5);
+
+orders = 1:5;
+derivatives = arrayfun(@(i) g.diff(g.domain, i), orders);
+verifyEqual(testCase, [derivatives.at(0), derivatives.at(g.T)], zeros( 1, 2*length(orders)), 'AbsTol', 1e-13)
+
+end
+
+function testBasicVariableInit(testCase)
+
+%%
+t = quantity.Domain("t", linspace(0, 1, 11)');
+g1 = signals.GevreyFunction('order', 1.8, 'timeDomain', t);
+g2 = signals.GevreyFunction('order', 1.5, 'timeDomain', t);
+
+G = [g1; g2];
+N_diff = 2;
+derivatives = G.diffs(0:N_diff);
+
+GD = zeros(t.n, N_diff+1, 2);
+for k = 1:N_diff+1
+ GD(:, k, 1) = signals.gevrey.bump.g(t.grid, k-1, t.grid(end), g1.sigma);
+ GD(:, k, 2) = signals.gevrey.bump.g(t.grid, k-1, t.grid(end), g2.sigma);
+end
+
+testCase.verifyEqual(derivatives.on(), GD);
+end
+
+function basicVariable_Diff_test(tc)
+ t = sym('t');
+ T = quantity.Domain('t', linspace(0,1));
+ fun = t^2 * (1-t)^2;
+
+ F = quantity.Symbolic(fun, T);
+ d1F = F.diff("t", 1);
+ d2F = F.diff("t", 2);
+
+ b(1,:) = signals.BasicVariable(F, {d1F, d2F});
+ b(2,:) = signals.BasicVariable(F, {d1F, d2F});
+
+ tc.verifyEqual(b.diff(T, 0), [F; F]);
+ tc.verifyEqual(b.diff(T, 1).on(), on(d1F * ones(2,1)));
+ tc.verifyEqual(b.diff(T, 2).on(), on(d2F * ones(2,1)));
+ tc.verifyError(@() b.diff(T, 3), 'conI:signals:BasicVariable:derivative')
+end
+
+function testGevreyFunction(testCase)
+
+N = 3;
+tau = quantity.Domain.defaultDomain(42, "tau");
+g(1,1) = signals.GevreyFunction('order', 1.5, 'gain', -5, 'offset', 2.5, 'numDiff', N, ...
+ 'timeDomain', tau, 'diffShift', 0);
+
+g(2,1) = signals.GevreyFunction('order', 1.9, 'gain', 5, 'offset',- 2.5, 'numDiff', N, ...
+ 'timeDomain', tau, 'diffShift', 0);
+
+testCase.verifyEqual(g.diff(tau, 0).at(0), [2.5; -2.5])
+testCase.verifyEqual(g.diff(tau, 0).at(1), [-2.5; 2.5])
+
+% test all the precomputed derivatives plus one higher derivative
+for k = 1:N+1
+ testCase.verifyEqual(g.diff(tau, k).at(0), zeros(2,1))
+ testCase.verifyEqual(g.diff(tau, k).at(1), zeros(2,1))
+end
+
+end
+
+function testGevreyFunctionTimes(testCase)
+
+N = 3;
+tau = quantity.Domain.defaultDomain(42, "tau");
+g(1,1) = signals.GevreyFunction('order', 1.5, 'gain', -5, 'offset', 2.5, 'numDiff', N, ...
+ 'timeDomain', tau, 'diffShift', 0);
+
+g(2,1) = signals.GevreyFunction('order', 1.9, 'gain', 5, 'offset',- 2.5, 'numDiff', N, ...
+ 'timeDomain', tau, 'diffShift', 0);
+
+A = [1 2; 3 4];
+
+b = A*g;
+testCase.verifyEqual(b.diff(tau, 0).at(0), A * [2.5; -2.5])
+testCase.verifyEqual(b(1).highestDerivative, g(1).highestDerivative)
+
+end
+
+function testBasicVariableFractionalDerivative(testCase)
+
+ t = sym('t');
+ t0 = 0;
+ timeDomain = quantity.Domain('t', linspace(t0, t0+1, 1e2));
+
+ %% verify singularity case
+ beta = 1.5;
+ F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
+ b(1,:) = signals.BasicVariable(F, {F.diff("t", 1), F.diff("t", 2), F.diff("t", 3)});
+ p = 0.3;
+ testCase.verifyError( @() b.diff(timeDomain, p), "conI:BasicVariable:fDiff")
+
+
+ %% verify regular case
+ for beta = [1, 2, 3]
+ F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
+ b(1,:) = signals.BasicVariable(F, {F.diff("t", 1), F.diff("t", 2), F.diff("t", 3)});
+
+ alpha = 0.3;
+ for k = 1:5
+ p = k*alpha;
+ % compute the fractional derivative with the formula for the polynomial (see Podlubny: 2.3.4)
+ f_diff(k,1) = quantity.Symbolic( fractional.monomialFractionalDerivative(t, t0, beta, p), timeDomain);
+
+ % compute the fractional derivative numerically
+ b_diff(k,1) = b.diff(timeDomain, p);
+ end
+ % plot(f_diff - b_diff)
+
+ testCase.verifyEqual( b_diff.on(), f_diff.on(), 'AbsTol', 2e-2)
+ end
+
+ %% verify fractional integration
+ testCase.verifyEqual( b.diff(timeDomain, 0).on(), F.on() )
+
+ for k = 1:5
+ alpha = 0.3 * k;
+ J = int( (t- sym("tau")).^(alpha -1) * subs(F.sym, "t", "tau"), "tau", t0, t) / gamma(alpha);
+ J_num = b.diff( timeDomain, -alpha );
+
+ J = quantity.Symbolic(J, timeDomain);
+
+ testCase.verifyEqual( J.on() , J_num.on, 'AbsTol', 3e-3 );
+ end
+
+end
+
+function testBasicVariableFractionalProductRule(testCase)
+
+t = quantity.Domain("t", linspace(0,1));
+% choose a polynomial basis
+N = 5;
+
+
+beta = 1;
+h(1,1) = quantity.Symbolic(sym("t").^beta, t);
+poly(1,1) = quantity.Symbolic( sym("t"), t);
+for k = 2:N
+ h(k,1) = h(k-1).diff();
+ poly(k,1) = quantity.Symbolic( sym("t")^(k-1) / factorial(k-1), t);
+end
+
+b = signals.BasicVariable( h(1), num2cell(h(2:end)) );
+
+for k = 1:3
+
+ f = sym(b.diff(t, 0)) * sym(poly(k));
+
+ alpha = 0.3 * k;
+ n = ceil(alpha);
+
+ frac_f = quantity.Symbolic(fractional.derivative(f, alpha, t.lower, t.name), t);
+ frac_p = b.productRule(poly(k), alpha);
+
+ testCase.verifyEqual(frac_f.on(), frac_p.on(), 'AbsTol', 4e-3)
+end
+
+end
+
+
+function test_BasicVariable_ProductRule(testCase)
+
+t = quantity.Domain("t", linspace(0,1,11));
+phi = quantity.Symbolic( sin( sym("t") * pi), t);
+bfun = quantity.Symbolic( cos(sym("t") * pi), t);
+b = signals.BasicVariable( bfun, { bfun.diff(t, 1), bfun.diff(t, 2), bfun.diff(t, 3) });
+
+for k = 0:3
+ h1 = b.productRule(phi, k);
+ h2 = diff( phi * bfun, t, k);
+ testCase.verifyEqual( h1.on, h2.on, 'AbsTol', 1e-13 );
+end
+
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-end
\ No newline at end of file
diff --git a/+unittests/+signals/testBasicVariable.m b/+unittests/+signals/testBasicVariable.m
deleted file mode 100644
index 0efc4a94cf92f3dd7477404a33fd24d9bc656697..0000000000000000000000000000000000000000
--- a/+unittests/+signals/testBasicVariable.m
+++ /dev/null
@@ -1,26 +0,0 @@
-function [tests] = testBasicVariable()
-%testQuantity Summary of this function goes here
-% Detailed explanation goes here
-tests = functiontests(localfunctions);
-end
-
-function testDiff(tc)
- t = sym('t');
- T = quantity.Domain('t', linspace(0,1));
- fun = t^2 * (1-t)^2;
-
- F = quantity.Symbolic([fun; fun], 'domain', T);
- d1F = F.diff("t", 1);
- d2F = F.diff("t", 2);
-
-
- b = signals.BasicVariable(F, {d1F, d2F});
-
- tc.verifyEqual(b.diff(0), F);
- tc.verifyEqual(b.diff(1), d1F);
- tc.verifyEqual(b.diff(2), d2F);
-end
-
-function testApplyPolynomialOperator(testCase)
- %op =
-end
\ No newline at end of file
diff --git a/+unittests/+signals/testPolynomialOperator.m b/+unittests/+signals/testPolynomialOperator.m
index 32dcc78f9da73f0de1da19ebc9c0c40e81dd5ae9..b83c50ce0f237e40201dfaf4178abf965b78ee5b 100644
--- a/+unittests/+signals/testPolynomialOperator.m
+++ b/+unittests/+signals/testPolynomialOperator.m
@@ -8,7 +8,7 @@ function setupOnce(testCase)
%%
z = sym('z', 'real');
s = sym('s');
-Z = linspace(0, 1, 11)';
+Z = quantity.Domain("z", linspace(0, 1, 11));
% init with symbolic values
A = cat(3, [ 0, 0, 0, 0; ...
@@ -27,7 +27,7 @@ A = cat(3, [ 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, 'gridName', 'z');
+ a{k} = quantity.Symbolic(A(:,:,k), Z);
end
A = signals.PolynomialOperator(a, 's', s);
@@ -49,7 +49,7 @@ function testApplyTo(testCase)
t = sym('t', 'real');
sol.y = quantity.Symbolic( t.^3 .* (1 - t).^3, ...
- 'grid', linspace(0, 1, 1e2), 'gridName', 't');
+ quantity.Domain("t", linspace(0, 1, 1e2)));
sol.u = detsIA.applyTo(sol.y);
sol.x = adjB.applyTo(sol.y);
@@ -71,9 +71,9 @@ function testFundamentalMatrixSpaceDependent(testCase)
%
a = 20;
z = sym('z', 'real');
- Z = linspace(0, 1, 5)';
+ Z = quantity.Domain("z", linspace(0, 1, 5));
A = signals.PolynomialOperator(...
- {quantity.Symbolic([1, z; 0, a], 'gridName', 'z', 'grid', Z)});
+ {quantity.Symbolic([1, z; 0, a], Z)});
F = A.stateTransitionMatrix();
@@ -84,7 +84,7 @@ function testFundamentalMatrixSpaceDependent(testCase)
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, 'gridName', 'z', 'name', 'Phi_A');
+ PhiExact = quantity.Symbolic([exp(z), f; 0, exp(a*z)], Z, 'name', 'Phi_A');
testCase.verifyEqual(double(F), PhiExact.on(), 'RelTol', 1e-6);
@@ -92,11 +92,11 @@ end
%
function testStateTransitionMatrix(testCase)
N = 2;
-Z = linspace(0,1,11)';
-A0 = quantity.Symbolic([0 1; 1 0], 'grid', Z, 'gridName', 'z');
-A1 = quantity.Symbolic([0 1; 0 0], 'grid', Z, 'gridName', 'z');
+z = quantity.Domain("z", linspace(0,1,11));
+A0 = quantity.Symbolic([0 1; 1 0], z);
+A1 = quantity.Symbolic([0 1; 0 0], z);
A = signals.PolynomialOperator({A0, A1});
-B = signals.PolynomialOperator(quantity.Symbolic([-1 -1; 0 0], 'grid', Z, 'gridName', 'z'));
+B = signals.PolynomialOperator(quantity.Symbolic([-1 -1; 0 0], z));
[Phi1, Psi1] = A.stateTransitionMatrix(N, B);
[Phi2, Psi2] = A.stateTransitionMatrixByOdeSystem(N, B);
@@ -108,9 +108,7 @@ end
function testChangeGrid(testCase)
A = testCase.TestData.A;
a = testCase.TestData.a;
- A0 = A.changeGrid({0}, {'z'});
-
- testCase.verifyEqual(cellfun(@(v) numel(v), A0(1).grid), 1);
+ A0 = A.subs('z', 0);
for k = 1:3
testCase.verifyEqual( ...
@@ -123,11 +121,11 @@ function testMTimes(testCase)
zGrid = linspace(0, pi, 3)';
z = quantity.Domain("z", zGrid);
A0 = quantity.Discrete(reshape((repmat([1 2; 3 4], length(zGrid), 1)), length(zGrid), 2, 2),...
- 'domain', z, 'name', 'A');
+ z, 'name', 'A');
A1 = quantity.Discrete(reshape((repmat([5 6; 7 8], length(zGrid), 1)), length(zGrid), 2, 2), ...
- 'domain', z, 'name', 'A');
+ z, 'name', 'A');
A2 = quantity.Discrete(reshape((repmat([9 10; 11 12], length(zGrid), 1)), length(zGrid), 2, 2), ...
- 'domain', z, 'name', 'A');
+ z, 'name', 'A');
A = signals.PolynomialOperator({A0, A1, A2});
B = signals.PolynomialOperator({A0});
@@ -157,13 +155,13 @@ function testMTimes(testCase)
end
function testSum(testCase)
-z = linspace(0, pi, 5)';
- A0 = quantity.Discrete(reshape((repmat([1 2; 3 4], length(z), 1)), length(z), 2, 2),...
- 'gridName', 'z', 'grid', z, 'name', 'A');
- A1 = quantity.Discrete(reshape((repmat([5 6; 7 8], length(z), 1)), length(z), 2, 2), ...
- 'gridName', 'z', 'grid', z, 'name', 'A');
- A2 = quantity.Discrete(reshape((repmat([9 10; 11 12], length(z), 1)), length(z), 2, 2), ...
- 'gridName', 'z', 'grid', z, 'name', 'A');
+ z = quantity.Domain("z", linspace(0, pi, 5));
+ A0 = quantity.Discrete(reshape((repmat([1 2; 3 4], z.n, 1)), z.n, 2, 2),...
+ z, 'name', 'A');
+ A1 = quantity.Discrete(reshape((repmat([5 6; 7 8], z.n, 1)), z.n, 2, 2), ...
+ z, 'name', 'A');
+ A2 = quantity.Discrete(reshape((repmat([9 10; 11 12], z.n, 1)), z.n, 2, 2), ...
+ z, 'name', 'A');
A = signals.PolynomialOperator({A0, A1, A2});
B = signals.PolynomialOperator({A0});
@@ -215,24 +213,12 @@ function testUminus(testCase)
%%
syms z s
Z = quantity.Domain('z', linspace(0,1,11));
-A0 = quantity.Symbolic( [z, z.^2; z.^3, z.^4], 'domain', Z);
-A1 = quantity.Symbolic( [sin(z), cos(z); tan(z), sinh(z)], 'domain', Z);
+A0 = quantity.Symbolic( [z, z.^2; z.^3, z.^4], Z);
+A1 = quantity.Symbolic( [sin(z), cos(z); tan(z), sinh(z)], Z);
o = signals.PolynomialOperator({A0, A1});
mo = -o;
%%
verifyEqual(testCase, o(1).coefficient.on , -mo(1).coefficient.on);
-end
-
-
-
-
-
-
-
-
-
-
-
-
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/+unittests/+signals/testTimeDelayOperator.m b/+unittests/+signals/testTimeDelayOperator.m
index 63def61783f0e5bccd055cf631ed5909dbb6aed5..83026e6aeb6c4f2c6f0019dc9512ddd0ad0d9478 100644
--- a/+unittests/+signals/testTimeDelayOperator.m
+++ b/+unittests/+signals/testTimeDelayOperator.m
@@ -9,9 +9,9 @@ function setupOnce(testCase)
z = quantity.Domain('z', linspace(0,1));
zeta = quantity.Domain('zeta', linspace(0,1));
t = quantity.Domain('t', linspace(0, 2, 31));
-v = quantity.Discrete( z.grid, 'domain', z );
-c_zeta = quantity.Discrete( zeta.grid * 0, 'domain', zeta );
-c = quantity.Discrete( 1, 'domain', quantity.Domain.empty);
+v = quantity.Discrete( z.grid, z );
+c_zeta = quantity.Discrete( zeta.grid * 0, zeta );
+c = quantity.Discrete( 1, quantity.Domain.empty);
testCase.TestData.z = z;
testCase.TestData.t = t;
@@ -70,9 +70,9 @@ function testApplyTo(testCase)
tau = quantity.Domain('tau', linspace(-1, 3, 201));
v = testCase.TestData.delay.variable;
- h = quantity.Discrete( sin(tau.grid * 2 * pi), 'domain', tau );
+ h = quantity.Discrete( sin(tau.grid * 2 * pi), tau );
H = quantity.Discrete( sin(z.grid * 2 * pi + t.grid' * 2 * pi), ...
- 'domain', [z, t]);
+ [z, t]);
H_ = v.applyTo(h);
%% test a scalar operator
@@ -85,9 +85,9 @@ function testApplyTo(testCase)
%% test a diagonal operator
h = quantity.Discrete( {sin(tau.grid * 2 * pi); ...
- sin(tau.grid * 2 * pi)}, 'domain', tau );
+ sin(tau.grid * 2 * pi)}, tau );
H = quantity.Discrete( {sin(z.grid * 0 + 2 * pi + t.grid' * 2 * pi); ...
- sin(z.grid * 2 * pi + t.grid' * 2 * pi)}, 'domain', [z, t]);
+ sin(z.grid * 2 * pi + t.grid' * 2 * pi)}, [z, t]);
D = testCase.TestData.delay.diagonal;
H_ = D.applyTo(h);
@@ -99,11 +99,9 @@ function testApplyTo(testCase)
% g(z) = z
% D[h] = z + sin( ( t + z ) * 2 * pi )
- h = quantity.Discrete( z.grid + sin(tau.grid' * 2 * pi), ...
- 'domain', [z, tau]);
+ h = quantity.Discrete( z.grid + sin(tau.grid' * 2 * pi), [z, tau]);
- Dh = quantity.Discrete( z.grid + sin( ( z.grid + t.grid') * 2 * pi ), ...
- 'domain', [z, t]);
+ Dh = quantity.Discrete( z.grid + sin( ( z.grid + t.grid') * 2 * pi ), [z, t]);
D = testCase.TestData.delay.variable;
diff --git a/+unittests/testFractional.m b/+unittests/testFractional.m
new file mode 100644
index 0000000000000000000000000000000000000000..1ebc3e1da83040c3fa8d2aeeed35406fa52c0474
--- /dev/null
+++ b/+unittests/testFractional.m
@@ -0,0 +1,36 @@
+function [tests] = testFractional()
+%UNTITLED Summary of this function goes here
+% Detailed explanation goes here
+tests = functiontests(localfunctions);
+end
+
+function testDerivative(testCase)
+t = sym("t");
+t0 = 0.5;
+
+ for p = [0.1 0.4 0.7 0.9]
+
+ % define the analytic solution for the fractional derivative
+ [y, m] = fractional.monomialFractionalDerivative( t, t0, 1, p);
+
+ % compute the fractional derivative with the function
+ f = fractional.derivative(m, p, t0, t);
+
+ % evaluate both and compare them
+ dt = 0.1;
+ % skip the first value because there is a singularity that can not be computed by the
+ % default symbolic evaluation.
+ % TODO@ferdinand: allow the quantity.Symbolic to use the limit function for the evaluation
+ % at certain points.
+ tDomain = quantity.EquidistantDomain("t", t0+dt, t0+1, "stepSize", dt);
+
+ Y = quantity.Symbolic( y, tDomain);
+ F = quantity.Symbolic( f, tDomain);
+
+ testCase.verifyEqual( Y.on(), F.on(), 'AbsTol', 1e-15);
+
+ h = fractional.derivative(m, p, t0, t, 'type', 'C');
+ H = quantity.Symbolic( h, tDomain);
+ testCase.verifyEqual( H.on(), Y.on(), 'AbsTol', 1e-15);
+ end
+end
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 43f5a61208957257007dc077b3ccbdcf34925299..2809e14d921afd710c454a6fd6274640916379cb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@
***.aux
***.log
.PowerFolder/
+*.asv