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

Moved the quantity.Operator to signals.Operator

parent 1d3f5e37
......@@ -208,12 +208,12 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
o = reshape(o, size(obj));
end
function o = quantity.Operator(obj)
function o = signals.PolynomialOperator(obj)
A = cell(size(obj, 3), 1);
for k = 1:size(obj, 3)
A{k} = obj(:,:,k);
end
o = quantity.Operator(A);
o = signals.PolynomialOperator(A);
end
function o = quantity.Symbolic(obj)
......@@ -1836,7 +1836,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
% result = DIFF(obj, k, diffGridName) applies the
% 'k'th-derivative for the variable specified with the input
% 'diffGridName' to the obj. If no 'diffGridName' is specified,
% then diff applies the derivative w.r.t. to all gridNames.
% then diff applies the derivative w.r.t. all gridNames.
if nargin == 1 || isempty(k)
k = 1; % by default, only one derivatve per diffGridName is applied
else
......@@ -1852,7 +1852,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
if nargin < 3 % if no diffGridName is specified, then the derivative
% w.r.t. to all gridNames is applied
% w.r.t. all gridNames is applied
diffGridName = obj(1).gridName;
end
......
......@@ -3,7 +3,6 @@ classdef TimeDelayOperator < handle & matlab.mixin.Copyable
% D[h](z,t) = h( t + coefficient(z) )
% This comes from the inverse Laplace transfromed value of
% exp(-s coefficient(z) ) * H(s) <-> h(t + coefficient(z))
%
properties
coefficient {mustBe.quantity};
......
classdef Operator < handle & matlab.mixin.Copyable
classdef PolynomialOperator < handle & matlab.mixin.Copyable
%OPERATOR define operator matrices
% This class describes operators of the form
% A[x(t)](z) = A0(z) x(z,t) + A1(z) dz x(z,t) + A2(z) dz^2 x(z,t)
......@@ -18,7 +18,7 @@ classdef Operator < handle & matlab.mixin.Copyable
gridName;
end
methods
function obj = Operator(A, varargin)
function obj = PolynomialOperator(A, varargin)
% OPERATOR initialization of an operator object
% obj = Operator(A, varargin) initializes an operator of the
% form
......@@ -102,7 +102,7 @@ classdef Operator < handle & matlab.mixin.Copyable
for k = 1:size(obj, 1)
II{k} = int(obj(k).coefficient, varargin{:});
end
I = quantity.Operator(II);
I = signals.PolynomialOperator(II);
end
function value = applyTo(obj, y, varargin)
......@@ -125,20 +125,19 @@ classdef Operator < handle & matlab.mixin.Copyable
for k = 1:n
value = value + obj(k).coefficient * y.diff(k - 1, myDomain );
end
end
function C = mtimes(A, B)
% if isnumeric(A) && numel(A) == 1
% A = quantity.Operator(eye(size(B(1).coefficient, 1)) * A);
% A = signals.PolynomialOperator(eye(size(B(1).coefficient, 1)) * A);
% end
%
if ~isa(A, 'quantity.Operator')
A = quantity.Operator(A);
if ~isa(A, 'signals.PolynomialOperator')
A = signals.PolynomialOperator(A);
end
if ~isa(B, 'quantity.Operator')
B = quantity.Operator(B);
if ~isa(B, 'signals.PolynomialOperator')
B = signals.PolynomialOperator(B);
end
[n, m] = size(A(1).coefficient * B(1).coefficient);
......@@ -161,16 +160,16 @@ classdef Operator < handle & matlab.mixin.Copyable
end
end
end
C = quantity.Operator(C);
C = signals.PolynomialOperator(C);
end
function C = plus(A, B)
if ~isa(A, 'quantity.Operator')
A = quantity.Operator(A);
if ~isa(A, 'signals.PolynomialOperator')
A = signals.PolynomialOperator(A);
end
if ~isa(B, 'quantity.Operator')
B = quantity.Operator(B);
if ~isa(B, 'signals.PolynomialOperator')
B = signals.PolynomialOperator(B);
end
for k = 1 : max(length(A), length(B))
......@@ -187,7 +186,7 @@ classdef Operator < handle & matlab.mixin.Copyable
end
end
end
C = quantity.Operator(C);
C = signals.PolynomialOperator(C);
end
function C = minus(A, B)
......@@ -197,20 +196,20 @@ classdef Operator < handle & matlab.mixin.Copyable
function C = adj(A)
assert(A(1).coefficient.isConstant())
C = polyMatrix.polynomial2coefficients(misc.adj(A.powerSeries), A(1).s);
C = polyMatrix.polynomial2coefficients(misc.adj(sym( A )), A(1).s);
c = cell(1, size(C, 3));
for k = 1:size(C, 3)
c{k} = quantity.Discrete(double(C(:,:,k)), ...
'gridName', {}, 'grid', {}, 'name', ['adj(' A(1).coefficient.name ')']);
end
C = quantity.Operator(c);
C = signals.PolynomialOperator(c);
end
function C = det(A)
assert(A(1).coefficient.isConstant())
C = det(A.powerSeries);
C = det( sym(A) );
C = polyMatrix.polynomial2coefficients(C, A(1).s);
c = cell(1, size(C, 3));
for k = 1:size(C, 3)
......@@ -218,7 +217,7 @@ classdef Operator < handle & matlab.mixin.Copyable
'gridName', {}, 'grid', {}, 'name', ['det(' A(1).coefficient.name ')']);
end
C = quantity.Operator(c);
C = signals.PolynomialOperator(c);
end
function C = uminus(A)
......@@ -239,9 +238,10 @@ classdef Operator < handle & matlab.mixin.Copyable
function d = double(obj)
d = double(obj.M);
end
function [i, m, s] = near(obj, B, varargin)
if isa(B, "quantity.Operator")
if isa(B, "signals.PolynomialOperator")
b = B.M;
else
b = B;
......@@ -252,6 +252,7 @@ classdef Operator < handle & matlab.mixin.Copyable
i = s;
end
end
function i = isempty(obj)
i = numel(obj)==0 || isempty(obj(1).coefficient);
end
......@@ -261,14 +262,19 @@ classdef Operator < handle & matlab.mixin.Copyable
end
function S = powerSeries(obj)
S =polyMatrix.powerSeries2Sum(squeeze(obj.M.on()), obj(1).s);
warning('DEPRICATED! use PolynomialOperator.sym() instead!')
S = polyMatrix.powerSeries2Sum(squeeze(obj.M.on()), obj(1).s);
end
function S = sym(obj)
S = polyMatrix.powerSeries2Sum( obj.M.on() , obj(1).s);
end
function A = ctranspose(obj)
for k = 1:length(obj)
M{k} = (-1)^(k-1) * obj(k).coefficient'; %#ok<AGROW>
end
A = quantity.Operator(M);
A = signals.PolynomialOperator(M);
end
function [Phi, Psi] = stateTransitionMatrixByOdeSystem(A, varargin)
......@@ -280,7 +286,7 @@ classdef Operator < handle & matlab.mixin.Copyable
% 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 quantity.Operator object by
% if B(z,s) is defined as signals.PolynomialOperator object by
% name-value-pair.
% To solve this problem the recursion formula from [1]
% dz w_k(z) = sum_j=0 ^d A_j(z) * w_{k-j}(z) + B_k nu
......@@ -318,7 +324,7 @@ classdef Operator < handle & matlab.mixin.Copyable
ip = misc.Parser();
ip.addOptional('N', 1);
ip.addOptional('B', quantity.Operator.empty(1, 0));
ip.addOptional('B', signals.PolynomialOperator.empty(1, 0));
ip.parse(varargin{:});
B = ip.Results.B;
N = ip.Results.N;
......@@ -365,8 +371,8 @@ classdef Operator < handle & matlab.mixin.Copyable
Psi(:,:,k) = p(idxk,:);
end
Phi = quantity.Operator(Phi);
Psi = quantity.Operator(Psi);
Phi = signals.PolynomialOperator(Phi);
Psi = signals.PolynomialOperator(Psi);
end
function [Phi, Psi, F0] = stateTransitionMatrix(A, varargin)
......@@ -379,14 +385,14 @@ classdef Operator < handle & matlab.mixin.Copyable
% 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 quantity.Operator object by
% if B(z,s) is defined as signals.PolynomialOperator object by
% name-value-pair.
n = size(A(1).coefficient, 2);
ip = misc.Parser();
ip.addOptional('N', 1);
ip.addOptional('B', quantity.Operator.empty(1, 0));
ip.addOptional('B', signals.PolynomialOperator.empty(1, 0));
ip.addParameter('F0', []);
ip.parse(varargin{:});
B = ip.Results.B;
......@@ -473,12 +479,12 @@ classdef Operator < handle & matlab.mixin.Copyable
% TODO: warning schreiben, dass überprüft ab welcher Ordnung
% die Matrizen nur noch numerisches Rauschen enthalten!
Phi = quantity.Operator(Phi);
Psi = quantity.Operator(Psi);
Phi = signals.PolynomialOperator(Phi);
Psi = signals.PolynomialOperator(Psi);
end
function newOperator = changeGrid(obj, gridNew, gridNameNew)
newOperator = quantity.Operator(obj.M.changeGrid(gridNew, gridNameNew));
newOperator = signals.PolynomialOperator(obj.M.changeGrid(gridNew, gridNameNew));
end
function newOperator = subs(obj, varargin)
......@@ -486,7 +492,7 @@ classdef Operator < handle & matlab.mixin.Copyable
if isnumeric(newOperator)
newOperator = squeeze( num2cell(newOperator, [1 2]) );
newOperator = quantity.Operator(newOperator);
newOperator = signals.PolynomialOperator(newOperator);
end
end
......@@ -501,7 +507,7 @@ classdef Operator < handle & matlab.mixin.Copyable
methods (Static)
function e = empty(varargin)
o = quantity.Discrete.empty(varargin{:});
e = quantity.Operator(o);
e = signals.PolynomialOperator(o);
end
end
......
function [tests] = testOperator()
function [tests] = testPolynomialOperator()
%TESTGEVREY Summary of this function goes here
% Detailed explanation goes here
tests = functiontests(localfunctions);
......@@ -30,7 +30,7 @@ for k = 1:3
a{k} = quantity.Symbolic(A(:,:,k), 'grid', Z, 'gridName', 'z');
end
A = quantity.Operator(a, 's', s);
A = signals.PolynomialOperator(a, 's', s);
testCase.TestData.A = A;
testCase.TestData.a = a;
......@@ -44,8 +44,8 @@ function testApplyTo(testCase)
C = [1 0 0];
sys = ss(A, B, C, []);
detsIA = quantity.Operator(num2cell(flip(charpoly(A))));
adjB = quantity.Operator( adjoint(sym('s')*eye(size(A)) - A) * B );
detsIA = signals.PolynomialOperator(num2cell(flip(charpoly(A))));
adjB = signals.PolynomialOperator( adjoint(sym('s')*eye(size(A)) - A) * B );
t = sym('t', 'real');
sol.y = quantity.Symbolic( t.^3 .* (1 - t).^3, ...
......@@ -57,8 +57,8 @@ function testApplyTo(testCase)
[sim.y, sim.t, sim.x] = lsim(sys, sol.u.on(), sol.y.domain.grid, 'foh');
% deviation of final values:
dev.y = sim.y(end)
dev.x = sim.x(end, :)
dev.y = sim.y(end);
dev.x = sim.x(end, :);
end
......@@ -66,14 +66,13 @@ function testCTranspose(testCase)
At = testCase.TestData.A';
testCase.verifyTrue(At(3).coefficient == testCase.TestData.A(3).coefficient');
end
%
function testFundamentalMatrixSpaceDependent(testCase)
%
a = 20;
z = sym('z', 'real');
Z = linspace(0, 1, 5)';
A = quantity.Operator(...
A = signals.PolynomialOperator(...
{quantity.Symbolic([1, z; 0, a], 'gridName', 'z', 'grid', Z)});
F = A.stateTransitionMatrix();
......@@ -96,8 +95,8 @@ 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');
A = quantity.Operator({A0, A1});
B = quantity.Operator(quantity.Symbolic([-1 -1; 0 0], 'grid', Z, 'gridName', 'z'));
A = signals.PolynomialOperator({A0, A1});
B = signals.PolynomialOperator(quantity.Symbolic([-1 -1; 0 0], 'grid', Z, 'gridName', 'z'));
[Phi1, Psi1] = A.stateTransitionMatrix(N, B);
[Phi2, Psi2] = A.stateTransitionMatrixByOdeSystem(N, B);
......@@ -119,7 +118,6 @@ function testChangeGrid(testCase)
end
end
function testMTimes(testCase)
z = linspace(0, pi, 3)';
......@@ -130,8 +128,8 @@ function testMTimes(testCase)
A2 = quantity.Discrete(reshape((repmat([9 10; 11 12], length(z), 1)), length(z), 2, 2), ...
'gridName', 'z', 'grid', z, 'name', 'A');
A = quantity.Operator({A0, A1, A2});
B = quantity.Operator({A0});
A = signals.PolynomialOperator({A0, A1, A2});
B = signals.PolynomialOperator({A0});
C = A*B;
......@@ -166,8 +164,8 @@ z = linspace(0, pi, 5)';
A2 = quantity.Discrete(reshape((repmat([9 10; 11 12], length(z), 1)), length(z), 2, 2), ...
'gridName', 'z', 'grid', z, 'name', 'A');
A = quantity.Operator({A0, A1, A2});
B = quantity.Operator({A0});
A = signals.PolynomialOperator({A0, A1, A2});
B = signals.PolynomialOperator({A0});
C = A + B;
......@@ -183,7 +181,48 @@ z = linspace(0, pi, 5)';
end
function testDet(testCase)
% test the det of OP via the characteristic polynomial of a matrix A
% det(sI - A) = det(OP)
A = [0 -1; 1 0];
OP = signals.PolynomialOperator( {-A, [1 0; 0 1]} );
d1 = det(OP);
testCase.verifyEqual( charpoly([0 -1; 1 0]), double( [d1.coefficient] ) );
end
function testAdj(testCase)
% adj(sI - A) / det(sI - A) = (sI - A)^-1
A = [0 -1; 1 0];
OP = signals.PolynomialOperator( {-A, eye(2) } );
invA = sym( adj(OP) ) / sym( det(OP) );
invSymA = inv( OP(1).s * eye(2) - A );
testCase.verifyEqual( invA, invSymA )
end
function testInitialization(testCase)
A = signals.PolynomialOperator( {[0 -1; 1 0], [0 1; 0 0]} );
verifyClass(testCase, A, 'signals.PolynomialOperator');
end
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);
o = signals.PolynomialOperator({A0, A1});
mo = -o;
%%
verifyEqual(testCase, o(1).coefficient.on , -mo(1).coefficient.on);
end
......
function [tests] = testQuantityOpertor()
%TESTQUANTITYOPERTOR Summary of this function goes here
% Detailed explanation goes here
tests = functiontests(localfunctions);
end
function testTimes(testCase)
%%
syms z s
A0 = [z, z.^2; z.^3, z.^4];
A1 = [sin(z), cos(z); tan(z), sinh(z)];
o = quantity.Operator(A0 + A1*s, 'grid', {linspace(0,1)'});
end
function testDet(testCase)
%%
A = quantity.Operator( cat(3, [0 -1; 1 0], [0 1; 0 0]), 'grid', {0} );
d1 = det(A);
d2 = det(A.valueContinuous);
%%
verifyEqual(testCase, d1.valueContinuous, d2);
end
function testAdj(testCase)
A = quantity.Operator( cat(3, [0 -1; 1 0], [0 1; 0 0]), 'grid', {0} );
a1 = adj(A);
a2 = misc.adj(A.valueContinuous);
%%
verifyEqual(testCase, a1.valueContinuous, a2);
end
function testInitialization(testCase)
A = quantity.Operator( cat(3, [0 -1; 1 0], [0 1; 0 0]), 'grid', {0} );
verifyClass(testCase, A, 'quantity.Operator');
end
function testVariable(testCase)
syms z s
o = quantity.Operator([z * s; 1 + s], 'grid', {linspace(0,1)'});
verifyEqual(testCase, o.variable, z);
verifyEqual(testCase, o.operatorVar, s);
end
function testEvaluate(testCase)
%%
syms z s
A0 = [z, z.^2; z.^3, z.^4];
A1 = [sin(z), cos(z); tan(z), sinh(z)];
o = quantity.Operator(A0 + A1*s, 'grid', {linspace(0,1)'});
value = o.valueDiscrete();
Z = o.grid{1};
%%
for k = 1:2
for l = 1:2
verifyLessThan(testCase, max(value(:,k,l,1) - double(subs(A0(k,l), z, Z))), 2 * eps);
verifyLessThan(testCase, max(value(:,k,l,2) - double(subs(A1(k,l), z, Z))), 2 * eps);
end
end
end
function testUminus(testCase)
%%
syms z s
A0 = [z, z.^2; z.^3, z.^4];
A1 = [sin(z), cos(z); tan(z), sinh(z)];
o = quantity.Operator(A0 + A1*s, 'grid', {linspace(0,1)'});
mo = -o;
value = o.valueDiscrete();
mValue = mo.valueDiscrete();
Z = o.grid{1};
%%
verifyEqual(testCase, value, -mValue);
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