Commit b6ab6f11 authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Moved helper functions for polynomial matrices to the new package folder +polyMatrix

parent 8ebc5ee0
function [ APoly ] = polynomial2coefficients(A, s_var, N)
warning('DEPRICATED: Use polyMatrix.polynomial2coefficients instead')
% OPERATOR2POLYNOMIAL Returns the coefficient matrices of an operator
%
% Inputs:
......
function [ S ] = powerSeries2Sum( M, s_var )
warning('DEPRICATED: Use polyMatrix.powerSeries2Sum instead')
% POWERSERIES2SUM Sum of power series in s_var
%
% Description:
......
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
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
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
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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)
......
......@@ -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
......@@ -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
......
......@@ -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
......@@ -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
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