Skip to content

GitLab

Commit f69f944e authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

Merge branch 'master' of https://gitlab.cs.fau.de/lrt_infinite_dimensional_systems/coni

parents fcfe4b03 58354f82
Hide whitespace changes
Inline Side-by-side
+misc/+ss/planPolynomialTrajectory.m
View file @ f69f944e
......@@ -30,17 +30,16 @@ function [u, y, x] = planPolynomialTrajectory(sys, t, optArgs)
% T: the time for which the transition should occur
arguments
sys ss;
t double;
t quantity.Domain;
optArgs.x0 double = zeros(size(sys.A, 1), 1);
optArgs.x1 double = zeros(size(sys.A, 1), 1);
optArgs.domainName string = "t";
end
x0 = optArgs.x0;
x1 = optArgs.x1;
t0 = t(1);
t1 = t(end);
t0 = t.lower;
t1 = t.upper;
%% Definitions
A = sys.A;
......@@ -55,7 +54,7 @@ n_u = size(B,2);
d = sym('d',[n_u,(2*n)]);
% variables
s = sym('s');
tau = sym( optArgs.domainName );
t_sym = sym( t.name );
%% Computation of the PI matrix
% The PI matrix is a matrix that contains each coefficient of the PI_s
% matrix, the matrix PI_s is given as PI_s = adj(sI - A)*B
......@@ -79,7 +78,7 @@ for i = 1:n_u
temp = size(coef,2);
PI_i(j,1:size(coef,2)) = coef;
end
y0(i) = poly2sym(d(i,:),tau);
y0(i) = poly2sym(d(i,:),t_sym);
PI(i) = {PI_i};
end
PI = cell2mat(PI);
......@@ -102,20 +101,20 @@ c1 = pinv(PI,1e-12)*x1;
Y = sym('y',[n,n_u]);
Y(1,:) = y0;
for j = 1:(n-1)
Y(j+1,:) = diff(y0,tau,j);
Y(j+1,:) = diff(y0,t_sym,j);
end
% The actual computation of the coefficients d done for each element of
% y_schlange.
i_temp = 0;
for i = 1:n_u
EQN1 = subs(Y(:,i),tau,t0) == c0((i_temp + 1):(i_temp + n));
EQN1 = subs(Y(:,i),t_sym,t0) == c0((i_temp + 1):(i_temp + n));
EQN1 = EQN1.';
EQN2 = subs(Y(:,i),tau,t1) == c1((i_temp + 1):(i_temp + n));
EQN2 = subs(Y(:,i),t_sym,t1) == c1((i_temp + 1):(i_temp + n));
EQN2 = EQN2.';
[Ai,Bi] = equationsToMatrix([EQN1, EQN2], d(i,:));
coefpoly = linsolve(Ai,Bi);
y_schlange(i) = poly2sym(coefpoly,tau);
y_schlange(i) = poly2sym(coefpoly,t_sym);
i_temp = i_temp + n;
end
%% Computation of the input u, the states x and the output y = Cx
......@@ -133,19 +132,18 @@ end
Y_schlange = sym('y',[n,n_u]);
Y_schlange(1,:) = y_schlange;
for j = 1:(n-1)
Y_schlange(j+1,:) = diff(y_schlange,tau,j);
Y_schlange(j+1,:) = diff(y_schlange,t_sym,j);
end
% The computation of the input. Y_schlange must be extended to have one
% more row that containts the n-th derivative of y_schlange
u = flip(charpoly(A))* [Y_schlange;diff(y_schlange,tau,n)];
u = quantity.Symbolic( u.', 'grid', {t}, 'gridName', optArgs.domainName);
u = flip(charpoly(A))* [Y_schlange;diff(y_schlange,t_sym,n)];
u = quantity.Symbolic( u.', t);
% Y_schlange is now reshaped so that it corresponds to the size of the PI
% matrix and the calculation of the states follows directly.
Y_schlange = reshape(Y_schlange,n_u*n,1);
x = quantity.Symbolic( PI*Y_schlange, 'grid', {t}, 'gridName', optArgs.domainName);
y = quantity.Symbolic( sys.C*x.sym, ...
'grid', {t}, 'gridName', optArgs.domainName) + sys.D * u;
x = quantity.Symbolic( PI*Y_schlange, t);
y = quantity.Symbolic( sys.C*x.sym, t) + sys.D * u;
end
+misc/+ss/planTrajectory.m
View file @ f69f944e
function [u, y, x] = planTrajectory(sys, t, optArgs)
function [u, y, x, W] = planTrajectory(sys, t, optArgs)
% PLANTRAJECTORY Computes a trajectory for a lumped-parameter systems
% [u, y, x] = planTrajectory(sys, t, varargin) computes a transition
% [u, y, x] = planTrajectory(sys, t, varargin) computes a transition
% u(t) : x(t0) = x0 --> x(t1) = x1
% for a state-space system sys of the form
% d/dt x = A x + bu.
......@@ -15,40 +15,71 @@ function [u, y, x] = planTrajectory(sys, t, optArgs)
%
arguments
sys ss;
t double;
t quantity.Domain;
optArgs.x0 double = zeros(size(sys.A, 1), 1);
optArgs.x1 double = zeros(size(sys.A, 1), 1);
optArgs.domainName string = "t";
optArgs.w quantity.Discrete;
optArgs.method = "Chen3";
optArgs.weight = quantity.Symbolic.ones(1, t);
end
x0 = optArgs.x0;
x1 = optArgs.x1;
% prepare time vectors
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, tDomain));
invPhi_t1 = expm(sys.A * quantity.Discrete( t1 - t, tDomain));
if optArgs.method == "Bernstein"
t0 = t.lower;
t1 = t.upper;
% Compute the state transition matrix: Phi(t,tau) = expm(A*(t-tau) )
Phi_t0 = expm(sys.A * (t.Discrete - t.lower));
invPhi_t1 = expm(sys.A * (t.upper - t.Discrete));
% 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
W = cumInt( Phi_t0 * sys.b * sys.b' * expm(sys.A' * (t.Symbolic - t.lower)), t, t0, t.name);
W1_t1 = W.at(t1);
%Formel aus dem Bernstein:
u = sys.b.' * invPhi_t1.' / W1_t1 * (x1 - Phi_t0.at(t1) * x0);
%Berechnung der
x = Phi_t0 * x0 + W * invPhi_t1' / W1_t1 * (x1 - Phi_t0.at(t1) * x0);
end
%% compute the gramian controllability matrix
% W1 = misc.ss.gramian(sys, t, t0, optArgs.domainName);
%W2 = expm( sys.A * (t.Discrete - t1 ) ) * cumInt( invPhi_t1 * sys.B * sys.B' * invPhi_t1', t, t0, t.name);
%% Chen^3
if optArgs.method == "Chen3"
Phi_t_t0 = expm(sys.A * ( t.Symbolic - t.lower ));
Phi_t1_t = expm(sys.A * (t.upper - t.Symbolic ) );
Phi_t_t1 = expm(sys.A * (t.Symbolic - t.upper ) );
Wt0t = cumInt( Phi_t1_t * sys.B * sys.B' * Phi_t1_t' * optArgs.weight, t, t.lower, t.name );
Wt0t1 = Wt0t.at(t.upper);
u = - sys.B' * Phi_t1_t' / Wt0t1 * ( Phi_t_t0.at(t.upper) * x0 - x1 ) * optArgs.weight;
x = Phi_t_t0 * x0 - Phi_t_t1 * Wt0t / Wt0t1 * ( Phi_t_t0.at(t.upper) * x0 - x1 );
end
% 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, tDomain)), ...
Phi_t0(1).domain, t0, optArgs.domainName);
%% Chen^1
if optArgs.method == "Chen1"
%% compute the gramian controllability matrix
Phi_t_t0 = expm(sys.A * ( t.Symbolic - t.lower ));
Phi_t0_t = expm(sys.A * ( t.lower - t.Symbolic)); % = inv( Phi_t_t0 )
% int Phi(t,t0) * b * b^T * Phi(t,t0)' dt
W = cumInt( Phi_t0_t * sys.b * sys.b' * Phi_t0_t' * optArgs.weight, t, t.lower, t.name);
% solution for the control (see Chen (5-15))
% u(t) = - B' * Phi'(t0, t) * W(t0, t1)^-1 * ( x0 - Phi(t0, t1) * x1)
u = - sys.b.' * Phi_t0_t.' / W.at(t.upper) * (x0 - Phi_t0_t.at(t.upper) * x1) * optArgs.weight;
% solution for the state:
% x(t) = Phi(t,t0) * ( x0 - W(t0, t) * W^-1(t0,t1) * (x0 - Phi(t0, t1*x1) );
x = Phi_t_t0 * ( x0 - W / W.at(t.upper) * (x0 - Phi_t0_t.at(t.upper) * x1) );
%
else
error("Method " + optArgs.method + " not yet implemented")
end
W1_t1 = W1.at(t1);
% Formel aus dem Bernstein:
u = sys.b.' * invPhi_t1.' / W1_t1 * (x1 - Phi_t0.at(t1) * x0);
% Berechnung der
x = Phi_t0 * x0 + W1 * invPhi_t1' / W1_t1 * (x1 - Phi_t0.at(t1) * x0);
%
y = sys.C * x + sys.D * u;
%% Alternative Lösung nach Chen:
......@@ -56,6 +87,6 @@ y = sys.C * x + sys.D * u;
% W0 = invPhi_t0 * W1 * invPhi_t0';
% u1 = -sys.B' * invPhi_t0' / W0.at(t1) * (x0 - Phi_t1.at(t0) * x1);
% xW0 = Phi_t0 * ( x0 - W0 / W0.at(t1) * (x0 - Phi_t1.at(t0) * x1));
%
%
% % simulation results have shown, does the average of both soultions leads to a better than each.
% u = (u0 + u1) / 2;
+misc/binomial.m
View file @ f69f944e
function b = binomial(p, k)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
% BINOMIAL computation of the binomial coefficient
% b = binomial(p, k) computes the binomial coefficient:
% / p \ p!^p
% | | = -----------
% \ k / (p-k)! k!
% to use this for fractional calculus, the p can be an integer. Then, the faculty is replaced by the
% gamma function.
arguments
p (1,1) double;
k (1,1) double {mustBeInteger};
......
+numeric/femMatrices.m
View file @ f69f944e
function [M, K, L, P, PHI, dphi] = femMatrices(spatialDomain, optArgs)
function [M, K, L, P, PHI, dphi, Dz] = femMatrices(spatialDomain, optArgs)
%FEMMATRICES computes the discretization matrices obtained by FE-Method
% [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
......@@ -8,12 +8,16 @@ function [M, K, L, P, PHI, dphi] = femMatrices(spatialDomain, optArgs)
% 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) )
% Derivative Matrix
% Dz = int( phi(z) * dz phi'(z) dz
%
% PHI is the trial function used for the discretization and DPHI is the derivative of the trial
% function.
......@@ -102,12 +106,14 @@ dphidphi = reshape(dphi * dphi', 1, 4);
K0 = zeros(2,2);
K1 = zeros(2,2);
K2 = zeros(2,2);
Dz_k = zeros(2,2);
L0 = zeros(2,0);
P0 = zeros(0,2);
% initialize the assembly matrices:
M = zeros(N, N);
K = zeros(N, N);
Dz = zeros(N, N);
L = zeros(N, p);
P = zeros(m, N);
......@@ -141,7 +147,7 @@ for k = 1:N-1
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
% as the matrix dzphi*dzphi' is constant 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);
......@@ -159,12 +165,19 @@ for k = 1:N-1
P0 = reshape( numeric.trapz_fast(zk, c0'), m, 2);
end
if nargout >= 7
Dz_ = misc.multArray(phi, repmat( dphi', [optArgs.Ne,1]), 3, 3, 1);
Dz_k = reshape( numeric.trapz_fast(zk, reshape(Dz_, element.n, 4, 1)' ), 2, 2);
end
% assemble the matrices with the elements:
M(idxk, idxk) = M(idxk, idxk) + M0 * zDiff(k);
K(idxk, idxk) = K(idxk, idxk) + K0 + K1 - K2;
L(idxk, :) = L(idxk, :) + L0;
P(:, idxk) = P(:, idxk) + P0;
Dz(idxk, idxk) = Dz(idxk, idxk) + Dz_k;
end
......
+quantity/Discrete.m
View file @ f69f944e
......@@ -962,22 +962,23 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
[newValues, oldDomainName] = quantity.Discrete.subsParser(oldDomainName, varargin);
solution = copy(obj);
for it = 1 : numel(oldDomainName)
if isstring(newValues{it})
% if newValues{it} is a string, it is converted in a domain, to use the same
% implementation as for quantity.Domain-objects in newValues.
newValues{it} = solution(1).domain.find(oldDomainName(it)).rename(newValues{it}, oldDomainName(it));
end
if isnumeric(newValues{it})
solution = solution.subsNumeric(oldDomainName(it), newValues{it});
else
solution = solution.subsDomain(oldDomainName(it), newValues{it});
end % if-else
end % for it = 1 : numel(oldDomainName)
if ~isempty(solution)
for it = 1 : numel(oldDomainName)
if isstring(newValues{it})
% if newValues{it} is a string, it is converted in a domain, to use the same
% implementation as for quantity.Domain-objects in newValues.
newValues{it} = solution(1).domain.find(oldDomainName(it)).rename(newValues{it}, oldDomainName(it));
end
if isnumeric(newValues{it})
solution = solution.subsNumeric(oldDomainName(it), newValues{it});
else
solution = solution.subsDomain(oldDomainName(it), newValues{it});
end % if-else
end % for it = 1 : numel(oldDomainName)
end % if - isempty
end % subs()
function solution = subsNumeric(obj, oldDomainName, value)
......@@ -1118,7 +1119,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete ...
p.addParameter('titleWithIndex', true');
p.addParameter('hold', false);
p.addParameter('export', false);
p.addParameter('currentFigure', false);
p.addParameter('currentFigure', quantity.Settings.instance().plot.currentFigure);
p.addParameter('exportOptions', ...
{'height', [num2str(0.25*size(obj, 1)), '\textwidth'], ...
'width', '0.8\textwidth', 'externalData', false, ...
......
+quantity/Settings.m
View file @ f69f944e
......@@ -3,7 +3,7 @@ classdef Settings < handle
%
properties (Access = public)
plot = struct('dock', false, 'showTitle', true);
plot = struct('dock', false, 'showTitle', true, 'currentFigure', false);
progressBar = struct('show', true);
end
......
+signals/BasicVariable.m
View file @ f69f944e
......@@ -132,32 +132,41 @@ classdef BasicVariable < handle
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 );
function h = productRule(a, b, n)
% PRODUCTRULE compute the k-derivative of the product of the basic variable a and b
% h = productRule(a, b, n) computes the k-th derivative
% h = d_t^n ( a * b )
% using the product rule. This is usefull if the derivatives of the basic variable a
% are known exactly and the derivatives of b, but not the multiplication a*b
t = b.domain;
assert( size(a,2) == size(b,1) || all( size(b) == 1) || all( size(a) == 1) )
if size(a,2) == size(b,1)
h = quantity.Discrete.zeros( [size(a, 1), size(b,2)], t );
elseif all( size(a) == 1)
h = quantity.Discrete.zeros( size(b), t );
elseif all( size(b) == 1 )
h = quantity.Discrete.zeros( size(a), t );
else
error("The product rule for quantities of the size " + size(a) + " and " + size(b) + " is not possible.")
end
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);
h = h + misc.binomial(n, k) * a.diff(t, n-k) * b.diff(t, k);
end
else
% 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);
h = misc.binomial( n, 0 ) * b * a.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);
h0 = misc.binomial(n,k) * b.diff(t, k) * a.fDiff( n - k);
h = h + h0;
k = k+1;
end
......@@ -208,31 +217,60 @@ classdef BasicVariable < handle
% todo: verify that the pre computed values in b are used.
arguments
a double;
b signals.BasicVariable
a
b
end
assert( size(a,2) == size(b,1), "dimensions of the terms for multiplication do not match")
assert( numel(size(a)) <= 2 )
assert( numel(size(b)) <= 2 )
degree = max([b.highestDerivative]);
F = quantity.Discrete(b);
% do the multiplication for each derivative:
for i = 1:degree+1
tmp(i,:) = a * shiftdim(F(i,:));
if isa( a, "signals.BasicVariable")
I = a.domain;
elseif isa( b, "signals.BasicVariable")
I = b.domain;
end
if size(a,2) == size(b,1)
O = quantity.Discrete.zeros( [size(a, 1), size(b,2)], I );
elseif all( size(a) == 1)
O = quantity.Discrete.zeros( size(b), I );
elseif all( size(b) == 1 )
O = quantity.Discrete.zeros( size(a), I );
else
error("The product rule for quantities of the size " + size(a) + " and " + size(b) + " is not possible.")
end
% restore the result as signals.BasicVariable
for k = 1: (size(a,1) * size(b,2))
for l = 2:size(tmp,1)
derivs_{l-1} = tmp(l,k);
if isa( a, "double") && isa( b, "signals.BasicVariable" )
% assert( size(a,2) == size(b,1), "dimensions of the terms for multiplication do not match")
assert( numel(size(a)) <= 2 )
assert( numel(size(b)) <= 2 )
degree = max([b.highestDerivative]);
% do the multiplication for each derivative:
for i = 1:degree+1
diffs{i} = a * b.diff(I, i-1);
end
c(k) = signals.BasicVariable( tmp(1,k), derivs_);
elseif isa( a, "signals.BasicVariable") && isa( b, "double" )
c = ( b.' * a.' ).';
return;
elseif isa( a, "signals.BasicVariable" ) && isa( b, "signals.BasicVariable" )
assert( a(1).domain.isequal(b(1).domain) , "Both BasicVariables must have the same domain");
order = max( [a.highestDerivative], [b.highestDerivative] );
for n = 0 : order
diffs{n+1} = O;
for k = 0 : n
diffs{n+1} = diffs{n+1} + nchoosek(n, k) * a.diff(I, k) * b.diff(I, n-k);
end
end
else
error("The multiplication of " + class(a) + " and " + class(b) + " is not yet implemented");
end
c = reshape(c, [size(a,1), size(b,2)]);
c = signals.BasicVariable(diffs{1}, diffs(2:end));
c = reshape(c, size(O) );
end
end % methods (Access = public)
......
+signals/PolynomialOperator.m
View file @ f69f944e
......@@ -426,8 +426,9 @@ classdef (InferiorClasses = {?quantity.Discrete, ?quantity.Function, ?quantity.S
Psi = signals.PolynomialOperator(Psi);
end
function [Phi, Psi, F0] = stateTransitionMatrix(obj, optArgs)
function [Phi, F0] = stateTransitionMatrix(obj, optArgs)
% STATETRANSITIONMATRIX computation of the
% TODO
% state-transition matrix
% [Phi, Psi] = stateTransitioMatrix(obj, varargin) computes
% the state-transition matrix of the boundary value problem
......@@ -441,25 +442,18 @@ classdef (InferiorClasses = {?quantity.Discrete, ?quantity.Function, ?quantity.S
arguments
obj signals.PolynomialOperator;
optArgs.N double = 1;
optArgs.B signals.PolynomialOperator = signals.PolynomialOperator.empty(1, 0);
optArgs.F0;
end
n = size(obj(1).coefficient, 2);
B = optArgs.B;
N = optArgs.N;
if isempty(B)
m = 0;
else
m = size(B(1).coefficient, 2);
end
nA = length(obj); % As order of a polynom is zero based
% initialization of the ProgressBar
counter = N - 1;
pbar = misc.ProgressBar(...
'name', 'Trajectory Planning (Woittennek)', ...
'name', 'Computation of the state transition matrix', ...
'terminalValue', counter, ...
'steps', counter, ...
'initialStart', true);
......@@ -492,22 +486,97 @@ classdef (InferiorClasses = {?quantity.Discrete, ?quantity.Function, ?quantity.S
end
z = F0(1).domain(1);
Phi(:,:,1) = F0.subs("zeta", 0);
On = quantity.Discrete.zeros([n n], z, "gridName", "z");
for k = 2 : N
pbar.raise();
% compute the temporal matrices:
sumPhi = On;
for l = 2 : min(nA, k)
sumPhi = sumPhi + obj(l).coefficient * Phi(:, :, k-l+1);
end
% compute the integration of the fundamental matrix with the temporal
% values:
% Phi_k(z) = int_0_z Phi_0(z, zeta) * M(zeta) d_zeta
Phi(:, :, k) = cumInt(F0 * sumPhi.subs("z", "zeta"), "zeta", z(1).lower, "z");
end
%
pbar.stop();
% TODO: warning schreiben, dass überprüft ab welcher Ordnung
% die Matrizen nur noch numerisches Rauschen enthalten!
Phi = signals.PolynomialOperator(Phi);
end
function [Psi] = inputTransitionMatrix(obj, B, Phi0, optArgs)
% STATETRANSITIONMATRIX computation of the
% TODO
% state-transition matrix
% [Phi, Psi] = stateTransitioMatrix(obj, varargin) computes
% the state-transition matrix of the boundary value problem
% dz x(z,s) = A(z,s) x(xi,s) + B(z,s) * u(s)
% so that Phi is the solution of
% dz Phi(z,xi,s) = A(z,s) Phi(z,xi,s)
% and Psi is given by
% Psi(z,xi,s) = int_xi_z Phi(z,zeta,s) B(zeta,s) d zeta
% if B(z,s) is defined as signals.PolynomialOperator object by
% name-value-pair.
arguments
obj signals.PolynomialOperator;
B signals.PolynomialOperator
Phi0 quantity.Discrete;
optArgs.N double = 1;
optArgs.F0;
end
n = size(obj(1).coefficient, 2);
N = optArgs.N;
if isempty(B)
Psi(:,:,1) = quantity.Discrete.empty();
m = 0;
else
Psi(:,:,1) = cumInt( F0 * B(1).coefficient.subs("z", "zeta"), "zeta", z.lower, "z");
m = size(B(1).coefficient, 2);
end
nA = length(obj); % As order of a polynom is zero based
On = quantity.Discrete.zeros([n n], z);
Om = quantity.Discrete.zeros([n m], z);
% initialization of the ProgressBar
counter = N - 1;
pbar = misc.ProgressBar(...
'name', 'Computation of the state transition input matrix', ...
'terminalValue', counter, ...
'steps', counter, ...
'initialStart', true);
%% computation of transition matrix as power series
% The fundamental matrix is considered in the laplace space as the state
% transition matrix of the ODE with complex variable s:
% dz w(z,s) = A(z, s) w(z,s) * B(z) v(s),
% Hence, fundamental matrix is a solution of the ODE:
<