Commit 1c1ac195 authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Multiplication and product rule for basic variables

parent e76d20c0
...@@ -132,32 +132,41 @@ classdef BasicVariable < handle ...@@ -132,32 +132,41 @@ classdef BasicVariable < handle
function h = productRule(obj, phi, n) function h = productRule(a, b, n)
% LEIBNIZ compute the k-derivative of the product of the basic variable and phi % PRODUCTRULE compute the k-derivative of the product of the basic variable a and b
% h = productRule(obj, phi, n) computes the k-th derivative of obj.fun*phi, i.e., % h = productRule(a, b, n) computes the k-th derivative
% h = d_t^n ( obj.fun(t) * phi ) % h = d_t^n ( a * b )
% using the product rule. This is usefull if the derivatives of the basic variable obj % using the product rule. This is usefull if the derivatives of the basic variable a
% are known exactly and the derivatives of phi, but not the multiplication % are known exactly and the derivatives of b, but not the multiplication a*b
% obj.diff(0)*phi. t = b.domain;
t = phi.domain; assert( size(a,2) == size(b,1) || all( size(b) == 1) || all( size(a) == 1) )
h = quantity.Discrete.zeros( size(obj), t );
if size(a,2) == size(b,1)
h = quantity.Discrete.zeros( [size(a, 1), size(b,2)], t );
elseif all( size(a) == 1)
h = quantity.Discrete.zeros( size(b), t );
elseif all( size(b) == 1 )
h = quantity.Discrete.zeros( size(a), t );
else
error("The product rule for quantities of the size " + size(a) + " and " + size(b) + " is not possible.")
end
if numeric.near(round(n), n) if numeric.near(round(n), n)
% for integer order derivatives % for integer order derivatives
for k = 0:n for k = 0:n
h = h + misc.binomial(n, k) * obj.diff(t, n-k) * phi.diff(t, k); h = h + misc.binomial(n, k) * a.diff(t, n-k) * b.diff(t, k);
end end
else else
% for fractional order derivatives: [see Podlubny] % for fractional order derivatives: [see Podlubny]
% p \~~ oo / p \ (k) p-k % p \~~ oo / p \ (k) p-k
% D (phi(t) f(t)) = > | | | phi (t) D f(t) % D (phi(t) f(t)) = > | | | phi (t) D f(t)
% a t /__ k=0 \ k / a t % a t /__ k=0 \ k / a t
h = misc.binomial( n, 0 ) * phi * obj.fDiff(n); h = misc.binomial( n, 0 ) * b * a.fDiff(n);
h0 = h; h0 = h;
k = 1; k = 1;
% do the iteration on the loop until the difference is very small. % do the iteration on the loop until the difference is very small.
while max( abs( h0 ) ) / max( abs( h ) ) > 1e-9 while max( abs( h0 ) ) / max( abs( h ) ) > 1e-9
h0 = misc.binomial(n,k) * phi.diff(t, k) * obj.fDiff( n - k); h0 = misc.binomial(n,k) * b.diff(t, k) * a.fDiff( n - k);
h = h + h0; h = h + h0;
k = k+1; k = k+1;
end end
...@@ -208,31 +217,60 @@ classdef BasicVariable < handle ...@@ -208,31 +217,60 @@ classdef BasicVariable < handle
% todo: verify that the pre computed values in b are used. % todo: verify that the pre computed values in b are used.
arguments arguments
a double; a
b signals.BasicVariable b
end end
assert( size(a,2) == size(b,1), "dimensions of the terms for multiplication do not match") if isa( a, "signals.BasicVariable")
assert( numel(size(a)) <= 2 ) I = a.domain;
assert( numel(size(b)) <= 2 ) elseif isa( b, "signals.BasicVariable")
I = b.domain;
degree = max([b.highestDerivative]);
F = quantity.Discrete(b);
% do the multiplication for each derivative:
for i = 1:degree+1
tmp(i,:) = a * shiftdim(F(i,:));
end end
if size(a,2) == size(b,1)
O = quantity.Discrete.zeros( [size(a, 1), size(b,2)], I );
elseif all( size(a) == 1)
O = quantity.Discrete.zeros( size(b), I );
elseif all( size(b) == 1 )
O = quantity.Discrete.zeros( size(a), I );
else
error("The product rule for quantities of the size " + size(a) + " and " + size(b) + " is not possible.")
end
% restore the result as signals.BasicVariable if isa( a, "double") && isa( b, "signals.BasicVariable" )
for k = 1: (size(a,1) * size(b,2)) % assert( size(a,2) == size(b,1), "dimensions of the terms for multiplication do not match")
for l = 2:size(tmp,1) assert( numel(size(a)) <= 2 )
derivs_{l-1} = tmp(l,k); assert( numel(size(b)) <= 2 )
degree = max([b.highestDerivative]);
% do the multiplication for each derivative:
for i = 1:degree+1
diffs{i} = a * b.diff(I, i-1);
end end
c(k) = signals.BasicVariable( tmp(1,k), derivs_);
elseif isa( a, "signals.BasicVariable") && isa( b, "double" )
c = ( b.' * a.' ).';
return;
elseif isa( a, "signals.BasicVariable" ) && isa( b, "signals.BasicVariable" )
assert( a(1).domain.isequal(b(1).domain) , "Both BasicVariables must have the same domain");
order = max( [a.highestDerivative], [b.highestDerivative] );
for n = 0 : order
diffs{n+1} = O;
for k = 0 : n
diffs{n+1} = diffs{n+1} + nchoosek(n, k) * a.diff(I, k) * b.diff(I, n-k);
end
end
else
error("The multiplication of " + class(a) + " and " + class(b) + " is not yet implemented");
end end
c = reshape(c, [size(a,1), size(b,2)]);
c = signals.BasicVariable(diffs{1}, diffs(2:end));
c = reshape(c, size(O) );
end end
end % methods (Access = public) end % methods (Access = public)
......
...@@ -57,6 +57,28 @@ for k = 0:3 ...@@ -57,6 +57,28 @@ for k = 0:3
testCase.verifyEqual( h1.on, h2.on, 'AbsTol', 1e-13 ); testCase.verifyEqual( h1.on, h2.on, 'AbsTol', 1e-13 );
end end
%% test product rule for a scalar and a vector:
a = [phi; bfun];
a = signals.BasicVariable( a, { a.diff(t, 1), a.diff(t, 2), a.diff(t, 3) });
for k = 0:3
c = a.productRule( bfun, k );
testCase.verifyEqual( c.on(), diff( a.diff(t,0) * bfun, t, k).on(), 'AbsTol', 1e-13);
c = b.productRule( a, k);
testCase.verifyEqual( c.on(), diff( bfun * a.diff(t,0), t, k).on(), 'AbsTol', 1e-13);
end
%% test product rule for two vectors:
for k = 0:3
c = a.productRule(a',k);
testCase.verifyEqual( c.on(), diff( a.diff(t,0) * a.diff(t,0).', t, k ).on(), 'AbsTol', 1e-13);
c = productRule(a', a, k);
testCase.verifyEqual( c.on(), diff( a.diff(t,0).' * a.diff(t,0), t, k ).on(), 'AbsTol', 1e-13);
end
end end
function diffShiftTest(testCase) function diffShiftTest(testCase)
...@@ -235,4 +257,25 @@ for k = 1:3 ...@@ -235,4 +257,25 @@ for k = 1:3
testCase.verifyEqual(frac_f.on(), frac_p.on(), 'AbsTol', 4e-3) testCase.verifyEqual(frac_f.on(), frac_p.on(), 'AbsTol', 4e-3)
end end
end
function testProductBasicVariables(testCase)
t = quantity.Domain("t", linspace(0,1,11));
b = quantity.Symbolic( [cos(sym("t") * pi); sin(sym("t"))], t);
b = signals.BasicVariable( b, { b.diff(t, 1), b.diff(t, 2), b.diff(t, 3) });
c = b.' * b;
for k = 0:3
h1 = productRule(b.', b, k);
testCase.verifyEqual( h1.on(), c.diff(c(1).domain, k).on(), 'AbsTol', 1e-13 );
end
c = b * b.';
for k = 0:3
h1 = productRule(b, b.', k);
testCase.verifyEqual( h1.on(), c.diff(c(1).domain, k).on(), 'AbsTol', 1e-13 );
end
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