From b6ab6f111e92dc1942326c79454afbca6899a0ce Mon Sep 17 00:00:00 2001
From: Ferdinand Fischer <ferdinand.fischer@fau.de>
Date: Fri, 6 Sep 2019 12:52:45 +0200
Subject: [PATCH] Moved helper functions for polynomial matrices to the new
package folder +polyMatrix
---
+misc/polynomial2coefficients.m | 1 +
+misc/powerSeries2Sum.m | 1 +
+polyMatrix/cTranspose.m | 13 +++++++
+polyMatrix/polynomial2coefficients.m | 52 ++++++++++++++++++++++++++
+polyMatrix/powerSeries2Sum.m | 40 ++++++++++++++++++++
+quantity/BasicVariable.m | 54 ++++++++++++++++++---------
+quantity/Discrete.m | 8 ++--
+quantity/Function.m | 4 +-
+quantity/Operator.m | 6 +--
+signals/GevreyFunction.m | 40 ++++++++++++--------
+signals/PolyBump.m | 41 ++++++++++++++++----
+signals/polyBump.m | 2 +-
+unittests/signalsTest.m | 7 +++-
13 files changed, 216 insertions(+), 53 deletions(-)
create mode 100644 +polyMatrix/cTranspose.m
create mode 100644 +polyMatrix/polynomial2coefficients.m
create mode 100644 +polyMatrix/powerSeries2Sum.m
diff --git a/+misc/polynomial2coefficients.m b/+misc/polynomial2coefficients.m
index 5aa20e5..5ecc914 100644
--- a/+misc/polynomial2coefficients.m
+++ b/+misc/polynomial2coefficients.m
@@ -1,4 +1,5 @@
function [ APoly ] = polynomial2coefficients(A, s_var, N)
+warning('DEPRICATED: Use polyMatrix.polynomial2coefficients instead')
% OPERATOR2POLYNOMIAL Returns the coefficient matrices of an operator
%
% Inputs:
diff --git a/+misc/powerSeries2Sum.m b/+misc/powerSeries2Sum.m
index 8db229a..6ece677 100644
--- a/+misc/powerSeries2Sum.m
+++ b/+misc/powerSeries2Sum.m
@@ -1,4 +1,5 @@
function [ S ] = powerSeries2Sum( M, s_var )
+warning('DEPRICATED: Use polyMatrix.powerSeries2Sum instead')
% POWERSERIES2SUM Sum of power series in s_var
%
% Description:
diff --git a/+polyMatrix/cTranspose.m b/+polyMatrix/cTranspose.m
new file mode 100644
index 0000000..947da96
--- /dev/null
+++ b/+polyMatrix/cTranspose.m
@@ -0,0 +1,13 @@
+function [M_T] = cTranspose(M)
+%CTRANSPOSE complex conjugate transpose of M
+% [M_T] = cTranspose(M) returns the complex conjugate transpose of a
+% polynomial matrix M = sum M_k s^k
+% M^T_k = sum (-1)^k M_k^T s^k
+% M has to be given as 3-dimensional array with the dimension 3 as
+% polynomial dimension.
+M_T = zeros(size(M,2), size(M, 1), size(M, 3));
+
+for k = 1:size(M, 3)
+ M_T(:,:,k) = (-1)^(k-1) * M(:,:,k)';
+end
+
diff --git a/+polyMatrix/polynomial2coefficients.m b/+polyMatrix/polynomial2coefficients.m
new file mode 100644
index 0000000..5aa20e5
--- /dev/null
+++ b/+polyMatrix/polynomial2coefficients.m
@@ -0,0 +1,52 @@
+function [ APoly ] = polynomial2coefficients(A, s_var, N)
+% OPERATOR2POLYNOMIAL Returns the coefficient matrices of an operator
+%
+% Inputs:
+% A Operator matrix with polynomials A \in C[s]^(n x n)
+% s_var symbolic variable
+% N maximal degree to be regarded for the coefficients. If the
+% polynomial contains higher degree entries than *N*,
+% the values are cut off.
+% Outputs:
+% APoly Coefficients
+%
+% Example:
+% -------------------------------------------------------------------------
+% import misc.*
+% syms s;
+% A = [10, s; 2*s^2, 3*s^3]
+% polynomial2coefficients(A, s)
+%
+% >> ans(:,:,1) = [10, 0]
+% [ 0, 0]
+% >> ans(:,:,2) = [ 0, 1]
+% [ 0, 0]
+% >> ans(:,:,3) = [ 0, 0]
+% [ 2, 0]
+% >> ans(:,:,4) = [ 0, 0]
+% [ 0, 3]
+% -------------------------------------------------------------------------
+import misc.*
+if ~exist('N', 'var')
+ % determine maximal degree of spatial variable s
+ sDegreeMax = -1;
+ for i=1:numel(A)
+ sDegreeMax = max([sDegreeMax, length(coeffsAll(A(i),s_var))]);
+ end
+ N=sDegreeMax;
+end
+
+% convert system-operator into polynomal expression
+APoly = sym(zeros([size(A), N]));
+for i=1:size(A,1)
+ for j=1:size(A,2)
+ tmp = coeffsAll(A(i,j), s_var);
+ for k=1:N
+ if k>length(tmp)
+ APoly(i, j, k) = sym(0);
+ else
+ APoly(i, j, k) = sym(tmp(k));
+ end
+ end
+ end
+end
\ No newline at end of file
diff --git a/+polyMatrix/powerSeries2Sum.m b/+polyMatrix/powerSeries2Sum.m
new file mode 100644
index 0000000..8db229a
--- /dev/null
+++ b/+polyMatrix/powerSeries2Sum.m
@@ -0,0 +1,40 @@
+function [ S ] = powerSeries2Sum( M, s_var )
+% POWERSERIES2SUM Sum of power series in s_var
+%
+% Description:
+% [ S ] = powerSeries2Sum( M, s_var )
+% Summation of all entries in *M* as power series with the symbolic
+% varivable *s_var*. *M* can be a cell or a three dimensional array. For
+% the 3-D array the dimension 1 and 2 are for the matrices that sould be
+% summed up and dimension 3 defines the exponent of the series. Its index
+% k corresponds to the coefficient for s_var^(k-1).
+%
+% Inputs:
+% M Matrix of interest
+% s_var Symbolic variable in which S should be developed
+% Outputs:
+% S Summation of power Series
+%
+% Example:
+% -------------------------------------------------------------------------
+% import misc.*
+% syms s
+% M = cat(3, [0 1], [2 0], [0 3]);
+% powerSeries2Sum(M, s)
+% >> ans = [2*s, 3*s^2 + 1]
+% -------------------------------------------------------------------------
+
+if iscell(M)
+
+ S=M{1};
+ for k=2:size(M,1)
+ S=S+M{k}*s_var^(k-1);
+ end
+
+else
+
+ S = M(:,:,1);
+ for k=2:size(M,3)
+ S = S + M(:,:,k)*s_var^(k-1);
+ end
+end
\ No newline at end of file
diff --git a/+quantity/BasicVariable.m b/+quantity/BasicVariable.m
index d8e3c16..78a573f 100644
--- a/+quantity/BasicVariable.m
+++ b/+quantity/BasicVariable.m
@@ -33,7 +33,7 @@ function obj = BasicVariable(valueContinuous, varargin)
parentVarargin = {};
- if nargin > 0
+ if nargin > 0 && ~isempty(varargin)
% make default grid:
preParser = misc.Parser();
@@ -55,14 +55,20 @@ function obj = BasicVariable(valueContinuous, varargin)
preResults.N_diff = preParser.Results.N_diff;
end
- grid = {linspace(0, preParser.Results.T, N_t)'};
- varargin = [varargin{:}, 'grid', {grid}];
- parentVarargin = {valueContinuous, varargin{:}, 'gridName', {'t'}};
+ prsr = misc.Parser();
+ prsr.addParameter('grid', {linspace(0, preParser.Results.T, N_t)'} );
+ prsr.addParameter('gridName', {'t'});
+ prsr.parse(varargin{:});
+
+ uargin = misc.struct2namevaluepair(prsr.Unmatched);
+ pargin = misc.struct2namevaluepair(prsr.Results);
+
+ parentVarargin = [{valueContinuous}, uargin(:)', pargin(:)'];
end
obj@quantity.Function(parentVarargin{:});
- if nargin > 0
+ if nargin > 0 && ~isempty(varargin)
% first parser
p1 = misc.Parser();
p1.addParameter('K', 1);
@@ -75,25 +81,31 @@ function obj = BasicVariable(valueContinuous, varargin)
for parameter = fieldnames(preResults)'
obj.(parameter{1}) = preParser.Results.(parameter{1});
end
+
+ obj.derivatives = obj;
end
end
-
- function n = nargin(obj)
- n = 1;
- end
-
+
function D = diff(obj, K)
if nargin == 1
i = 1;
end
-
+
for i = 1:numel(K)
- % create a new object for the derivative of obj:
- D(:,i) = obj.copy();
- [D.name] = deal(['d/dt (' num2str(K(i)) ') ']);
- [D.valueDiscrete] = deal([]);
- for l = 1:numel(obj)
- D(l,i).diffShift = D(l, 1).diffShift + K(i);
+ k = K(i);
+
+ if size(obj.derivatives, 1) < k+1 || isempty(obj.derivatives(k+1,:))
+ % create a new object for the derivative of obj:
+ D = obj.copy();
+ [D.name] = deal(['d/dt (' num2str(k) ') ']);
+ [D.valueDiscrete] = deal([]);
+ for l = 1:numel(obj)
+ D(l).diffShift = D(l).diffShift + k;
+ end
+ obj.derivatives(k+1,:) = D;
+ D.valueDiscrete();
+ else
+ D(:,i) = obj.derivatives(k+1,:);
end
end
end
@@ -126,6 +138,14 @@ function obj = BasicVariable(valueContinuous, varargin)
b = obj.copy();
b.valueDiscrete = K / obj.K * obj.valueDiscrete();
b.K = K;
+ b.derivatives(1) = b;
+ for i = 2:numel(b.derivatives)
+ bi = b.derivatives(i).copy;
+ bi.valueDiscrete = K / obj.K * bi.valueDiscrete();
+ bi.K = K;
+ b.derivatives(i) = bi;
+ end
+
end
end
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 4ca6aef..784dc04 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -5,10 +5,10 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.SymbolicII}) Discret
% Name of the domains that generate the grid.
gridName {mustBe.unique};
-
+
% In this cell, already computed derivatives can be stored to avoid
% multiple computations of the same derivative.
- derivatives cell = {};
+ derivatives;
end
properties (Hidden, Access = protected, Dependent)
@@ -325,12 +325,12 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.SymbolicII}) Discret
% find the finest grid of all input quantities by comparing
% gridSize for each iteratively.
- if isempty(a)
+ if isempty(a) || isempty(a(1).grid)
if nargin > 1
[referenceGrid, referenceGridName] = varargin{1}.getFinestGrid(varargin{2:end});
else
referenceGrid = {};
- referenceGridName = '';
+ referenceGridName = '';
end
return;
else
diff --git a/+quantity/Function.m b/+quantity/Function.m
index 0815054..5c7cc55 100644
--- a/+quantity/Function.m
+++ b/+quantity/Function.m
@@ -136,9 +136,7 @@ classdef Function < quantity.Discrete
mObj = obj.copy();
for k = 1:numel(obj)
- assert(isequal(obj(k).gridName, {'z'}), ['This only works for gridName == z, ', ...
- 'otherwise valueContinuous is wrong']);
- mObj(k).valueContinuous = @(z) - obj(k).valueContinuous(z);
+ mObj(k).valueContinuous = @(varargin) - obj(k).valueContinuous(varargin{:});
mObj(k).valueDiscrete = - obj(k).valueDiscrete;
end
diff --git a/+quantity/Operator.m b/+quantity/Operator.m
index 4bdd04a..813e2e9 100644
--- a/+quantity/Operator.m
+++ b/+quantity/Operator.m
@@ -136,7 +136,7 @@ classdef Operator < handle & matlab.mixin.Copyable
function C = adj(A)
assert(A(1).coefficient.isConstant())
- C = misc.polynomial2coefficients(misc.adj(A.powerSeries), A(1).s);
+ C = polyMatrix.polynomial2coefficients(misc.adj(A.powerSeries), A(1).s);
c = cell(1, size(C, 3));
for k = 1:size(C, 3)
c{k} = quantity.Discrete(double(C(:,:,k)), ...
@@ -149,7 +149,7 @@ classdef Operator < handle & matlab.mixin.Copyable
assert(A(1).coefficient.isConstant())
C = det(A.powerSeries);
- C = misc.polynomial2coefficients(C, A(1).s);
+ C = polyMatrix.polynomial2coefficients(C, A(1).s);
c = cell(1, size(C, 3));
for k = 1:size(C, 3)
c{k} = quantity.Discrete(double(C(:,:,k)), ...
@@ -196,7 +196,7 @@ classdef Operator < handle & matlab.mixin.Copyable
end
function S = powerSeries(obj)
- S = misc.powerSeries2Sum(squeeze(obj.M.on()), obj(1).s);
+ S =polyMatrix.powerSeries2Sum(squeeze(obj.M.on()), obj(1).s);
end
function A = ctranspose(obj)
diff --git a/+signals/GevreyFunction.m b/+signals/GevreyFunction.m
index 822a238..bc1b546 100644
--- a/+signals/GevreyFunction.m
+++ b/+signals/GevreyFunction.m
@@ -4,37 +4,45 @@ classdef GevreyFunction < quantity.BasicVariable
properties
% Gevrey-order
- order double = -1;
+ order double = -1;
end
-properties (Dependent)
- sigma;
-end
+ properties (Dependent)
+ sigma;
+ end
methods
function obj = GevreyFunction(varargin)
- obj@quantity.BasicVariable(@signals.gevrey.bump.g, varargin{:});
- prsr = misc.Parser();
- prsr.addParameter('order', 1.5 );
- prsr.parse(varargin{:});
- obj.order = prsr.Results.order;
+
+ if nargin == 0
+ g = [];
+ else
+ g = @signals.gevrey.bump.g;
+ 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;
+ end
end
- function s = get.sigma(obj)
- s = 1 ./ ( obj.order - 1);
- end
- function set.sigma(obj, s)
- obj.order = 1 + 1/s;
- 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 (Access = protected)
function v = evaluateFunction(obj, grid)
v = obj.K * obj.valueContinuous(grid, obj.diffShift, obj.T, obj.sigma);
- end
+ end
end
end
diff --git a/+signals/PolyBump.m b/+signals/PolyBump.m
index a306a3e..a385b85 100644
--- a/+signals/PolyBump.m
+++ b/+signals/PolyBump.m
@@ -4,6 +4,7 @@ classdef PolyBump < quantity.BasicVariable
a double;
b double;
c double;
+ sgn double;
end
properties ( Access = protected, Hidden)
@@ -12,31 +13,55 @@ classdef PolyBump < quantity.BasicVariable
end
methods
-
- function obj = PolyBump(varargin)
+ function obj = PolyBump(varargin)
+ % make default grid:
+ preParser = misc.Parser();
+ preParser.addParameter('T', 1);
+ preParser.addParameter('N_t', 101);
+ preParser.addParameter('dt', []);
+ preParser.addParameter('z', linspace(0, 1)');
+ preParser.parse(varargin{:});
+
+ if ~isempty(preParser.Results.dt)
+ Nt = quantity.BasicVariable.setDt(preParser.Results.T, preParser.Results.dt);
+ else
+ Nt = preParser.Results.N_t;
+ end
+
prsr = misc.Parser();
prsr.addParameter('a', 1);
prsr.addParameter('b', 1);
prsr.addParameter('c', 0);
+ prsr.addParameter('grid', {preParser.Results.z, linspace(0, preParser.Results.T, Nt)'});
+ prsr.addParameter('gridName', {'z', 't'});
+ prsr.addParameter('name', 'phi');
+ prsr.addParameter('sgn', 1);
+ prsr.addParameter('T', preParser.Results.T);
+ prsr.addParameter('N_t', Nt);
prsr.parse(varargin{:});
uargin = misc.struct2namevaluepair(prsr.Unmatched);
+ argin = [uargin(:)', ...
+ {'grid', prsr.Results.grid, ...
+ 'gridName', prsr.Results.gridName, ...
+ 'name', prsr.Results.name, ...
+ 'T', preParser.Results.T, ...
+ 'N_T', Nt}];
- obj@quantity.BasicVariable(@signals.polyBump, uargin{:});
+ obj@quantity.BasicVariable(@signals.polyBump, argin{:});
obj.a = prsr.Results.a;
obj.b = prsr.Results.b;
obj.c = prsr.Results.c;
+ obj.sgn = prsr.Results.sgn;
end
-
-
-
end
methods (Access = protected)
- function v = evaluateFunction(obj, grid)
- v = obj.K * obj.valueContinuous(grid, obj.diffShift, obj.T, obj.a, obj.b, obj.c);
+ function v = evaluateFunction(obj, z, t)
+ v = obj.K * obj.valueContinuous( t + obj.sgn * obj.c * z, obj.diffShift, obj.T, obj.a, obj.b, obj.c);
end
+
end
end
diff --git a/+signals/polyBump.m b/+signals/polyBump.m
index 5f8b03a..8a9a29f 100644
--- a/+signals/polyBump.m
+++ b/+signals/polyBump.m
@@ -14,5 +14,5 @@ else
end
% restict the output to the considered domain
-out = o .* ( t >= c & t <= (T-c) );
+out = o .* ( t >= c & t <= (T - c) );
end
diff --git a/+unittests/signalsTest.m b/+unittests/signalsTest.m
index 451365d..0a4a178 100644
--- a/+unittests/signalsTest.m
+++ b/+unittests/signalsTest.m
@@ -35,10 +35,15 @@ function testGevreyBump(testCase)
% Test script to check if the class for the bumb function returns the right
% signals.
-g = quantity.BasicVariable(@signals.gevrey.bump.g, 'diffShift', 1);
+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);
%% Test scaling of Gevrey function
verifyEqual(testCase, b, g.on(), 'AbsTol', 1e-13);
+% test the derivatives of the bump
+g = signals.GevreyFunction('diffShift', 0);
+g.diff(0:5);
+
+
end
\ No newline at end of file
--
GitLab