Commit a2365056 authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files
parents 6e6785f9 8ded2d23
...@@ -20,13 +20,11 @@ p.parse(varargin{:}); ...@@ -20,13 +20,11 @@ p.parse(varargin{:});
x0 = p.Results.x0; x0 = p.Results.x0;
x1 = p.Results.x1; x1 = p.Results.x1;
% prepare time vectors
t = t(:);
t0 = t(1); t0 = t(1);
t1 = t(end); t1 = t(end);
if ~numeric.near(t0, 0)
warning('Not yet tested for non zero initial conditions');
end
%% Compute the state transition matrix: Phi(t,tau) = expm(A*(t-tau) ) %% 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'} ); % 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', {'t'} )); Phi_t0 = expm(sys.A * quantity.Discrete( t - t0, 'grid', {t}, 'gridName', {'t'} ));
......
function xFolded = fold(x,z0,dim,nFolded) function xFolded = fold(x,z0,dim,varargin)
% misc.fold fold discretized vector or array around folding point % misc.fold fold discretized vector or array around folding point
% %
% xFolded = FOLD(x,z0) folds the vector x which is discretized over linspace(0,1,ndisc) at the % xFolded = FOLD(x,z0) folds the vector x which is discretized over linspace(0,1,ndisc) at the
...@@ -17,6 +17,38 @@ function xFolded = fold(x,z0,dim,nFolded) ...@@ -17,6 +17,38 @@ function xFolded = fold(x,z0,dim,nFolded)
% created on 07.10.2019 by Simon Kerschbaum % created on 07.10.2019 by Simon Kerschbaum
found=0;
passed = [];
arglist = {'preserve','nFolded'};
% preserve: insert the folding point at both left and right part. Important for observer error!
% nFolded resolution of the folded states
if ~isempty(varargin)
if mod(length(varargin),2) % uneven number
error('When passing additional arguments, you must pass them pairwise!')
end
for index = 1:2:length(varargin) % loop through all passed arguments:
for arg = 1:length(arglist)
if strcmp(varargin{index},arglist{arg})
passed.(arglist{arg}) = varargin{index+1};
found=1;
break
end
end % for arg
% argument wasnt found in for-loop
if ~found
error([varargin{index} ' is not a valid property to pass to fold!']);
end
found=0; % reset found
end % for index
end
if ~isfield(passed,'preserve')
passed.preserve = 0;
end
if ~isfield(passed,'nFolded')
passed.nFolded = 0;
end
if isvector(x) if isvector(x)
x=x(:); x=x(:);
if dim == 2 if dim == 2
...@@ -26,8 +58,10 @@ function xFolded = fold(x,z0,dim,nFolded) ...@@ -26,8 +58,10 @@ function xFolded = fold(x,z0,dim,nFolded)
if nargin < 3 if nargin < 3
dim =1; dim =1;
end end
if nargin < 4 if passed.nFolded == 0
nFolded = size(x,dim); nFolded = size(x,dim);
else
nFolded = passed.nFolded;
end end
zDiscOld = linspace(0,1,size(x,dim)); zDiscOld = linspace(0,1,size(x,dim));
...@@ -35,9 +69,14 @@ function xFolded = fold(x,z0,dim,nFolded) ...@@ -35,9 +69,14 @@ function xFolded = fold(x,z0,dim,nFolded)
zLeft = zeros(1,nFolded); zLeft = zeros(1,nFolded);
zRight = zeros(1,nFolded); zRight = zeros(1,nFolded);
zRight(1,:) = linspace(z0,1,nFolded); if ~passed.preserve
dZ = diff(zRight(1,:)); zRight(1,:) = linspace(z0,1,nFolded);
zLeft(1,:) = linspace(0,z0-dZ(1),nFolded); dZ = diff(zRight(1,:));
zLeft(1,:) = linspace(0,z0-dZ(1),nFolded);
else
zRight(1,:) = linspace(z0,1,nFolded);
zLeft(1,:) = linspace(0,z0,nFolded);
end
sizX = size(x); sizX = size(x);
if dim~= 1 if dim~= 1
...@@ -61,14 +100,15 @@ function xFolded = fold(x,z0,dim,nFolded) ...@@ -61,14 +100,15 @@ function xFolded = fold(x,z0,dim,nFolded)
xFoldedResh(:,1) = xInterp.evaluate(fliplr(zLeft(1,:)),1:size(xResh,2)); xFoldedResh(:,1) = xInterp.evaluate(fliplr(zLeft(1,:)),1:size(xResh,2));
xFoldedResh(:,2) = xInterp.evaluate(zRight(1,:),1:size(xResh,2)); xFoldedResh(:,2) = xInterp.evaluate(zRight(1,:),1:size(xResh,2));
else else
xFoldedResh = zeros([size(xResh) 2]); xFoldedResh = zeros([nFolded size(xResh,2) 2]);
xFoldedResh(:,:,1) = xInterp.evaluate(fliplr(zLeft(1,:)),1:size(xResh,2)); xFoldedResh(:,:,1) = xInterp.evaluate(fliplr(zLeft(1,:)),1:size(xResh,2));
xFoldedResh(:,:,2) = xInterp.evaluate(zRight(1,:),1:size(xResh,2)); xFoldedResh(:,:,2) = xInterp.evaluate(zRight(1,:),1:size(xResh,2));
end end
if ~ismatrix(xPerm) if ~ismatrix(xPerm)
xFoldedPerm = reshape(xFoldedResh,[size(xPerm) 2]); sizXPerm = size(xPerm);
xFoldedPerm = reshape(xFoldedResh,[nFolded sizXPerm(2:end) 2]);
else else
xFoldedPerm = xFoldedResh; xFoldedPerm = xFoldedResh;
end end
......
...@@ -72,6 +72,24 @@ else ...@@ -72,6 +72,24 @@ else
LDisc = L; LDisc = L;
end end
end end
if isa(A0,'function_handle')
A0Disc = eval_pointwise(A0,zdisc);
else
if size(A0,1)~= ndisc % constant matrix
A0Disc(1:ndisc,:,:) = repmat(shiftdim(A0,-1),ndisc,1,1);
else
A0Disc = A0;
end
end
if isa(H,'function_handle')
HDisc = eval_pointwise(H,zdisc);
else
if size(H,1)~= ndisc % constant matrix
HDisc(1:ndisc,:,:) = repmat(shiftdim(H,-1),ndisc,1,1);
else
HDisc = H;
end
end
n = size(LDisc,2); n = size(LDisc,2);
T1 = [eye(n); zeros(n,n)]; % selection matrix of first solution state T1 = [eye(n); zeros(n,n)]; % selection matrix of first solution state
T2 = [zeros(n,n); eye(n,n)]; % selection matrix of second solution state T2 = [zeros(n,n); eye(n,n)]; % selection matrix of second solution state
...@@ -193,7 +211,7 @@ for jBlock = 1:length(sigMod.blockLength) ...@@ -193,7 +211,7 @@ for jBlock = 1:length(sigMod.blockLength)
for zeta_idx=1:z_idx for zeta_idx=1:z_idx
temp_prod(zeta_idx,:,:) = ... temp_prod(zeta_idx,:,:) = ...
reshape(Psi_mu(z_idx,zeta_idx,:,:),[2*n 2*n])... reshape(Psi_mu(z_idx,zeta_idx,:,:),[2*n 2*n])...
*T2*(reshape(LDisc(zeta_idx,:,:),[n n])\reshape(A0(zeta_idx,:,:),[n n]))*T1.'; *T2*(reshape(LDisc(zeta_idx,:,:),[n n])\reshape(A0Disc(zeta_idx,:,:),[n n]))*T1.';
end end
A0_til(z_idx,:,:) = trapz(zdisc(1:z_idx),temp_prod,1); A0_til(z_idx,:,:) = trapz(zdisc(1:z_idx),temp_prod,1);
Y(z_idx,:,:) = reshape(Psi_mu(z_idx,1,:,:),[2*n,2*n]) - reshape(A0_til(z_idx,:,:),[2*n 2*n]); Y(z_idx,:,:) = reshape(Psi_mu(z_idx,1,:,:),[2*n,2*n]) - reshape(A0_til(z_idx,:,:),[2*n 2*n]);
...@@ -226,7 +244,7 @@ for jBlock = 1:length(sigMod.blockLength) ...@@ -226,7 +244,7 @@ for jBlock = 1:length(sigMod.blockLength)
g2jk = G2*varphi_j_k; g2jk = G2*varphi_j_k;
for z_idx = 1:ndisc for z_idx = 1:ndisc
hjk(z_idx,:) = reshape(H(z_idx,:,:),[n nv])*varphi_j_k; hjk(z_idx,:) = reshape(HDisc(z_idx,:,:),[n nv])*varphi_j_k;
temp_prod = zeros(z_idx,2*n,1); % preallocate and reset temp_prod = zeros(z_idx,2*n,1); % preallocate and reset
for zeta_idx=1:z_idx for zeta_idx=1:z_idx
temp_prod(zeta_idx,:,1) = ... temp_prod(zeta_idx,:,1) = ...
...@@ -323,10 +341,10 @@ for i=1:n ...@@ -323,10 +341,10 @@ for i=1:n
end end
end end
for z_idx=1:ndisc for z_idx=1:ndisc
lhs(z_idx,:,:) = sd(L(z_idx,:,:))*sd(Xzz(z_idx,:,:))... lhs(z_idx,:,:) = sd(LDisc(z_idx,:,:))*sd(Xzz(z_idx,:,:))...
+ sd(M(z_idx,:,:))*sd(sol(z_idx,:,:))... + sd(muSys(z_idx,:,:))*sd(sol(z_idx,:,:))...
+ sd(sol(z_idx,:,:))*S ... + sd(sol(z_idx,:,:))*S ...
+ sd(A0(z_idx,:,:))*sd(sol(1,:,:)); + sd(A0Disc(z_idx,:,:))*sd(sol(1,:,:));
end end
% rhs = H; % rhs = H;
% err = lhs-rhs; % err = lhs-rhs;
......
...@@ -16,6 +16,9 @@ function xUnfolded = unfold(x,z0,dim,nNew) ...@@ -16,6 +16,9 @@ function xUnfolded = unfold(x,z0,dim,nNew)
% created on 07.10.2019 by Simon Kerschbaum % created on 07.10.2019 by Simon Kerschbaum
% warning('The folding point will appear twice! Think about fixing!')
% fixing is easy! Insert a spatial point very close to the folding point to achieve this
% double-point and then interpolate on a linear grid!!! SO easy!
if isvector(x) if isvector(x)
error('x needs at least a spatial coordinate and a folding coordinate with 2 entries.'); error('x needs at least a spatial coordinate and a folding coordinate with 2 entries.');
end end
...@@ -47,24 +50,32 @@ function xUnfolded = unfold(x,z0,dim,nNew) ...@@ -47,24 +50,32 @@ function xUnfolded = unfold(x,z0,dim,nNew)
xResh = reshape(xPerm,sizX(dim),[],2); xResh = reshape(xPerm,sizX(dim),[],2);
sizXResh = size(xResh); sizXResh = size(xResh);
xLeftInterp = numeric.interpolant(... xLeftInterp = numeric.interpolant(...
{linspace(zNew(z0Idx),0,nOld),... {linspace(z0,0,nOld),1:sizXResh(2)},...
1:sizXResh(2),1:2},xResh(:,:,1)); xResh(:,:,1));
xRightInterp = numeric.interpolant(... xRightInterp = numeric.interpolant(...
{linspace(zNew(z0Idx+1),1,nOld),... {linspace(z0+1e-13,1,nOld),1:sizXResh(2)},...
1:sizXResh(2),1:2},xResh(:,:,2)); xResh(:,:,2));
xUnfoldedResh = zeros(nNew,sizXResh(2)); xUnfoldedReshTemp = zeros(nNew,sizXResh(2));
else else
xResh = xPerm; xResh = xPerm;
sizXResh = size(xResh); sizXResh = size(xResh);
xLeftInterp = numeric.interpolant(... xLeftInterp = numeric.interpolant(...
{linspace(zNew(z0Idx),0,nOld),1:2},xResh(:,1)); {linspace(z0,0,nOld),1:2},xResh(:,1));
xRightInterp = numeric.interpolant(... xRightInterp = numeric.interpolant(...
{linspace(zNew(z0Idx+1),1,nOld),1:2},xResh(:,2)); {linspace(z0+1e-13,1,nOld),1:2},xResh(:,2));
xUnfoldedResh = zeros(nNew,1); xUnfoldedReshTemp = zeros(nNew,1);
end end
xUnfoldedResh(1:z0Idx,:) = xLeftInterp.eval({zNew(1:z0Idx),1:sizXResh(2)}); xUnfoldedReshTemp(1:z0Idx,:) = xLeftInterp.eval({linspace(0,z0,z0Idx),1:sizXResh(2)});
xUnfoldedResh(z0Idx+1:end,:) = xRightInterp.eval({zNew(z0Idx+1:end),1:sizXResh(2)}); xUnfoldedReshTemp(z0Idx+1:end,:) = xRightInterp.eval({linspace(z0+1e-13,1,nNew-z0Idx),1:sizXResh(2)});
% attention, xUnfoldedResh is defined on non-uniform grid:
xUnfoldedReshInterp = numeric.interpolant({...
[linspace(0,z0,z0Idx) linspace(z0+1e-13,1,nNew-z0Idx)],...
1:size(xUnfoldedReshTemp,2)
},...
xUnfoldedReshTemp);
xUnfoldedResh = xUnfoldedReshInterp.eval({zNew,1:size(xUnfoldedReshTemp,2)});
if ~ismatrix(xPerm) if ~ismatrix(xPerm)
xUnfoldedPerm = reshape(xUnfoldedResh,[nNew sizXPerm(2:end-1)]); xUnfoldedPerm = reshape(xUnfoldedResh,[nNew sizXPerm(2:end-1)]);
......
...@@ -27,6 +27,10 @@ properties (Dependent) ...@@ -27,6 +27,10 @@ properties (Dependent)
dt; dt;
end end
properties (Hidden)
derivatives;
end
methods methods
%% CONSTRUCTOR %% CONSTRUCTOR
function obj = BasicVariable(valueContinuous, varargin) function obj = BasicVariable(valueContinuous, varargin)
...@@ -41,24 +45,32 @@ function obj = BasicVariable(valueContinuous, varargin) ...@@ -41,24 +45,32 @@ function obj = BasicVariable(valueContinuous, varargin)
preParser.addParameter('N_t', 101); preParser.addParameter('N_t', 101);
preParser.addParameter('dt', []); preParser.addParameter('dt', []);
preParser.addParameter('N_diff', 1); preParser.addParameter('N_diff', 1);
preParser.addParameter('t', []);
preParser.parse(varargin{:}); preParser.parse(varargin{:});
if ~isempty(preParser.Results.dt) if ~isempty(preParser.Results.dt)
N_t = quantity.BasicVariable.setDt(preParser.Results.T, preParser.Results.dt); N_t = quantity.BasicVariable.setDt(preParser.Results.T, preParser.Results.dt);
preResults.T = preParser.Results.T; preResults.T = preParser.Results.T;
preResults.dt = preParser.Results.dt; preResults.dt = preParser.Results.dt;
preResults.N_diff = preParser.Results.N_diff; elseif ~isempty(preParser.Results.t)
else N_t = numel( preParser.Results.t );
preResults.T = preParser.Results.t(end);
preResults.dt = preParser.Results.t(2) - preParser.Results.t(1);
elseif ~isempty(preParser.Results.N_t)
N_t = preParser.Results.N_t; N_t = preParser.Results.N_t;
preResults.T = preParser.Results.T; preResults.T = preParser.Results.T;
preResults.N_t = preParser.Results.N_t; preResults.N_t = preParser.Results.N_t;
preResults.N_diff = preParser.Results.N_diff; else
error('No time domain specified!')
end end
preResults.N_diff = preParser.Results.N_diff;
prsr = misc.Parser(); prsr = misc.Parser();
prsr.addParameter('grid', {linspace(0, preParser.Results.T, N_t)'} ); prsr.addParameter('grid', {linspace(0, preParser.Results.T, N_t)'} );
prsr.addParameter('gridName', {'t'}); prsr.addParameter('gridName', {'t'});
prsr.parse(varargin{:}); preParserUnmateched = misc.struct2namevaluepair(preParser.Unmatched);
prsr.parse(preParserUnmateched{:});
uargin = misc.struct2namevaluepair(prsr.Unmatched); uargin = misc.struct2namevaluepair(prsr.Unmatched);
pargin = misc.struct2namevaluepair(prsr.Results); pargin = misc.struct2namevaluepair(prsr.Results);
...@@ -79,7 +91,7 @@ function obj = BasicVariable(valueContinuous, varargin) ...@@ -79,7 +91,7 @@ function obj = BasicVariable(valueContinuous, varargin)
obj.(parameter{1}) = p1.Results.(parameter{1}); obj.(parameter{1}) = p1.Results.(parameter{1});
end end
for parameter = fieldnames(preResults)' for parameter = fieldnames(preResults)'
obj.(parameter{1}) = preParser.Results.(parameter{1}); obj.(parameter{1}) = preResults.(parameter{1});
end end
obj.derivatives = obj; obj.derivatives = obj;
...@@ -87,30 +99,46 @@ function obj = BasicVariable(valueContinuous, varargin) ...@@ -87,30 +99,46 @@ function obj = BasicVariable(valueContinuous, varargin)
end end
function D = diff(obj, K) function D = diff(obj, K)
% DIFF compute the K-th derivative of obj
%
%
if nargin == 1 if nargin == 1
i = 1; K = 1;
end end
for j = 1:numel(obj) % for each requested derivative:
objj = obj(j);
for i = 1:numel(K) for i = 1:numel(K)
k = K(i); k = K(i);
if size(objj.derivatives, 1) < k+1 || isempty(objj.derivatives(k+1,:)) % for each component of this object:
% create a new object for the derivative of obj: for j = 1:numel(obj)
D = objj.copy(); obj_j = obj(j);
[D.name] = deal(['d/dt (' num2str(k) ') ']);
[D.valueDiscrete] = deal([]); % check if the derivative is not yet already computed:
for l = 1:numel(objj) if size(obj_j.derivatives, 1) < k+1 || isempty(obj_j.derivatives(k+1,:))
D(l).diffShift = D(l).diffShift + k; % create a new object for the derivative of obj:
Dij = obj_j.copy();
Dij.name = ['d/dt (' num2str(k) ') '];
Dij.valueDiscrete = [];
Dij.diffShift = Dij.diffShift + k;
% sort the created derivative into the list of computed
% derivatives
obj_j.derivatives(k+1,j) = Dij;
% call the function to evaluate the numerical values of
% this derivative
Dij.valueDiscrete;
Di(j) = Dij;
else % the derivative is already computed -> take this one:
Di(j) = obj_j.derivatives(k+1);
end end
objj.derivatives(k+1,:) = D;
D.valueDiscrete;
else
D(:,i) = objj.derivatives(k+1,:);
end end
Di = reshape(Di, [ 1 size(obj)]);
D(i,:) = Di;
end end
end
end end
%% PROPERTIES %% PROPERTIES
......
...@@ -5,10 +5,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi ...@@ -5,10 +5,6 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% Name of the domains that generate the grid. % Name of the domains that generate the grid.
gridName {mustBe.unique}; gridName {mustBe.unique};
% In this cell, already computed derivatives can be stored to avoid
% multiple computations of the same derivative.
derivatives;
end end
properties (Hidden, Access = protected, Dependent) properties (Hidden, Access = protected, Dependent)
...@@ -1060,7 +1056,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi ...@@ -1060,7 +1056,9 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
hold off; hold off;
end end
if obj.nargin() == 0 if isempty(obj(rowIdx, columnIdx, figureIdx))
warning('you are trying to plot an empty quantity');
elseif obj.nargin() == 0
plot(0, ... plot(0, ...
obj(rowIdx, columnIdx, figureIdx).valueDiscrete, ... obj(rowIdx, columnIdx, figureIdx).valueDiscrete, ...
additionalPlotOptions{:}); additionalPlotOptions{:});
...@@ -1214,7 +1212,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi ...@@ -1214,7 +1212,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
for it = 1 : numel(obj) for it = 1 : numel(obj)
newObj(it).valueDiscrete = obj(it).on(myGrid); newObj(it).valueDiscrete = obj(it).on(myGrid);
end end
[newObj.derivatives] = deal({}); %[newObj.derivatives] = deal({});
[newObj.grid] = deal(myGrid); [newObj.grid] = deal(myGrid);
end end
...@@ -1635,19 +1633,15 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi ...@@ -1635,19 +1633,15 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% ISEMPTY checks if the quantity object is empty % ISEMPTY checks if the quantity object is empty
% empty = isempty(obj) % empty = isempty(obj)
% check if there is any dimension which is zero % Check if there is any dimension which is zero
empty = any(size(obj) == 0); empty = any(size(obj) == 0);
% if the constructor is called without arguments, a % If the constructor is called without arguments, a
% quantity.Discrete is initialized without an initialization of % quantity.Discrete is created without an initialization of
% the grid. Thus, the grid is not a cell: % the grid. Thus, the quantity is not initialized and empty if
% the grid is not a cell.
empty = empty || ~iscell(obj(1).grid); empty = empty || ~iscell(obj(1).grid);
% in order to check this for arrays:
% scalar and array empty check has to be separated because the
% [obj.grid] does only work in the array case.
empty = empty || any(cellfun(@iscell, [obj.grid]));
end % isempty() end % isempty()
function P = ztzTimes(a, b) function P = ztzTimes(a, b)
......
...@@ -149,7 +149,6 @@ classdef Function < quantity.Discrete ...@@ -149,7 +149,6 @@ classdef Function < quantity.Discrete
end end
[mObj.name] = deal(['-' obj.name]); [mObj.name] = deal(['-' obj.name]);
[mObj.derivatives] = deal({}); % #FIXME copy derivatives
end end
function p = inner(A, B) function p = inner(A, B)
......
...@@ -7,6 +7,19 @@ end ...@@ -7,6 +7,19 @@ end
function testBasicVariableInit(testCase) function testBasicVariableInit(testCase)
%% %%
b = quantity.BasicVariable(@signals.gevrey.bump.g, 'order', 1.8, 'N_t', 101, 'T', 1); 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 end
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment