Commit 960b0569 authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Fixed quantity.*/diff to be callable with a quantity.Domain. Secondly, the...

Fixed quantity.*/diff to be callable with a quantity.Domain. Secondly,  the quantity.Discrete/diff and quantity.Symbolic/diff were merged so that the common part is in quantity.Discrete/diff and the particular part for each class is in quantity.*/diff_inner
parent 1d947711
......@@ -1742,6 +1742,16 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
diffGridName = obj(1).gridName;
end
if isa(diffGridName, 'quantity.Domain')
% a quantity.Domain is used instead of a grid name
% -> get the grid name from the domain object
%
% #todo@domain: make the default case to call with a
% quantity.Domain instead of a grid name. Then, the
% section about the grid selection and so can be simplified
%
diffGridName = {diffGridName.name};
end
% diff for each element of diffGridName (this is rather
% inefficient, but an easy implementation of the specification)
......@@ -1755,37 +1765,7 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end
end
else
gridSelector = strcmp(obj(1).gridName, diffGridName);
gridSelectionIndex = find(gridSelector);
spacing = gradient(obj(1).grid{gridSelectionIndex}, 1);
assert(numeric.near(spacing, spacing(1)), ...
'diff is currently only implemented for equidistant grid');
permutationVector = 1 : (numel(obj(1).grid)+ndims(obj));
objDiscrete = permute(obj.on(), ...
[permutationVector(gridSelectionIndex), ...
permutationVector(permutationVector ~= gridSelectionIndex)]);
if size(objDiscrete, 2) == 1
derivativeDiscrete = gradient(objDiscrete, spacing(1));
else
[~, derivativeDiscrete] = gradient(objDiscrete, spacing(1));
end
rePermutationVector = [2:(gridSelectionIndex), ...
1, (gridSelectionIndex+1):ndims(derivativeDiscrete)];
result = quantity.Discrete(...
permute(derivativeDiscrete, rePermutationVector), ...
'size', size(obj), 'grid', obj(1).grid, ...
'gridName', obj(1).gridName, ...
'name', ['(d_{', diffGridName, '}', obj(1).name, ')']);
if k > 1
% % if a higher order derivative is requested, call the function
% % recursivly until the first-order derivative is reached
result = result.diff(k-1, diffGridName);
end
result = obj.diff_inner(k, diffGridName);
end
end
......@@ -2241,6 +2221,40 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
end %% (Static)
methods(Access = protected)
function result = diff_inner(obj, k, diffGridName)
gridSelector = strcmp(obj(1).gridName, diffGridName);
gridSelectionIndex = find(gridSelector);
spacing = gradient(obj(1).grid{gridSelectionIndex}, 1);
assert(numeric.near(spacing, spacing(1)), ...
'diff is currently only implemented for equidistant grid');
permutationVector = 1 : (numel(obj(1).grid)+ndims(obj));
objDiscrete = permute(obj.on(), ...
[permutationVector(gridSelectionIndex), ...
permutationVector(permutationVector ~= gridSelectionIndex)]);
if size(objDiscrete, 2) == 1
derivativeDiscrete = gradient(objDiscrete, spacing(1));
else
[~, derivativeDiscrete] = gradient(objDiscrete, spacing(1));
end
rePermutationVector = [2:(gridSelectionIndex), ...
1, (gridSelectionIndex+1):ndims(derivativeDiscrete)];
result = quantity.Discrete(...
permute(derivativeDiscrete, rePermutationVector), ...
'size', size(obj), 'grid', obj(1).grid, ...
'gridName', obj(1).gridName, ...
'name', ['(d_{', diffGridName, '}', obj(1).name, ')']);
if k > 1
% % if a higher order derivative is requested, call the function
% % recursivly until the first-order derivative is reached
result = result.diff(k-1, diffGridName);
end
end
function [idx, permuteGrid] = computePermutationVectors(a, b)
% Computes the required permutation vectors to use
% misc.multArray for multiplication of matrices
......
......@@ -19,7 +19,25 @@ classdef Operator < handle & matlab.mixin.Copyable
end
methods
function obj = Operator(A, varargin)
% OPERATOR initialization of an operator object
% obj = Operator(A, varargin) initializes an operator 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) + Ak dz^k x(z,t)
% The input parameter A can be:
% cell-array of doubles
% this is good to initialize operators with constant
% coefficients
% A[x](t) = a_0 x(t) + a_1 dt x(t) + ...
%
% cell-array of quantities
% this helps to initialize operators with variable
% coefficients:
% A[x(t)](z) = A0(z) x(z,t) + A1(z) dz x(z,t) + ...
% A2(z) dz^2 x(z,t) + Ak dz^k x(z,t)
%
%
%
if nargin > 0
prsr = misc.Parser();
......@@ -71,6 +89,29 @@ classdef Operator < handle & matlab.mixin.Copyable
I = quantity.Operator(II);
end
function value = applyOn(obj, y, varargin)
% APPLYON evaluation of the differential parameterization
% VALUE = applyOn(OBJ, Y, varargin) applies the operator OBJ
% on the function Y, i.e.,
% value = A0(t) y(t) + A1(t) dt y(t) + A2(t) dz^2 y(t)...
% is evaluated.
myParser = misc.Parser();
myParser.addParameter('domain', y.domain);
myParser.addParameter('n', length(obj));
myParser.parse(varargin{:});
myDomain = myParser.Results.domain;
n = myParser.Results.n;
value = 0;
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
......@@ -106,6 +147,7 @@ classdef Operator < handle & matlab.mixin.Copyable
end
C = quantity.Operator(C);
end
function C = plus(A, B)
if ~isa(A, 'quantity.Operator')
......@@ -131,9 +173,11 @@ classdef Operator < handle & matlab.mixin.Copyable
end
C = quantity.Operator(C);
end
function C = minus(A, B)
C = A + (-1)*B;
end
function C = adj(A)
assert(A(1).coefficient.isConstant())
......@@ -146,6 +190,7 @@ classdef Operator < handle & matlab.mixin.Copyable
C = quantity.Operator(c);
end
function C = det(A)
assert(A(1).coefficient.isConstant())
......@@ -159,13 +204,16 @@ classdef Operator < handle & matlab.mixin.Copyable
C = quantity.Operator(c);
end
function C = uminus(A)
C = (-1)*A;
end
function s = size(obj, varargin)
s = [length(obj), size(obj(1).coefficient)];
s = s(varargin{:});
end
function l = length(obj)
l = numel(obj);
end
......@@ -416,6 +464,7 @@ classdef Operator < handle & matlab.mixin.Copyable
function newOperator = changeGrid(obj, gridNew, gridNameNew)
newOperator = quantity.Operator(obj.M.changeGrid(gridNew, gridNameNew));
end
function newOperator = subs(obj, varargin)
newOperator = obj.M.subs(varargin{:});
......
......@@ -531,34 +531,6 @@ classdef Symbolic < quantity.Function
'name', ['sqrtm(', x(1).name, ')']);
end % sqrtm()
function result = diff(obj, k, diffGridName)
% diff applies the kth-derivative for the variable specified with
% the input gridName to the obj. If no gridName is specified, then diff
% applies the derivative w.r.t. to all gridNames / variables.
if nargin == 1 || isempty(k)
k = 1; % by default, only one derivatve per diffGridName is applied
end
if nargin <= 2 % if no diffGridName is specified, then the derivative
% w.r.t. to all gridNames is applied
diffGridName = obj(1).gridName;
end
result = obj.copy();
if iscell(diffGridName)
for it = 1 : numel(diffGridName)
result = result.diff(k, diffGridName{it});
end
else
diffVariable = obj.gridName2variable(diffGridName);
[result.name] = deal(['(d_{' char(diffVariable) '} ' result(1).name, ')']);
[result.valueDiscrete] = deal([]);
for l = 1:numel(obj)
result(l).valueSymbolic = diff(obj(l).valueSymbolic, diffVariable, k);
result(l).valueContinuous = obj.setValueContinuous(result(l).valueSymbolic, obj(1).variable);
end
end
end % diff()
function b = flipGrid(a, myGridName)
if ~iscell(myGridName)
myGridName = {myGridName};
......@@ -892,6 +864,16 @@ classdef Symbolic < quantity.Function
end
end
function result = diff_inner(obj, k, diffGridName)
result = obj.copy();
[result.name] = deal(['(d_{' diffGridName '} ' result(1).name, ')']);
[result.valueDiscrete] = deal([]);
for l = 1:numel(obj)
result(l).valueSymbolic = diff(obj(l).valueSymbolic, diffGridName, k);
result(l).valueContinuous = obj.setValueContinuous(result(l).valueSymbolic, obj(1).variable);
end
end
end
methods (Static)
......
......@@ -35,6 +35,30 @@ A = quantity.Operator(a, 's', s);
testCase.TestData.A = A;
testCase.TestData.a = a;
end
function testApplyOn(testCase)
A = [0 1 0 0; 0 0 1 0; 0 0 0 1; 1:4];
B = [0 0 0 1]';
C = [1 0 0 0];
sys = ss(A, B, C, []);
G = tf(sys);
denom = flip( G.Denominator{:} );
det = quantity.Operator(num2cell(denom));
t = sym('t', 'real')
y = quantity.Symbolic( t.^3 .* (1 - t).^3, ...
'grid', linspace(0,1), 'gridName', 't');
u = det.applyOn(y)
end
function testCTranspose(testCase)
......
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