Commit 3d80341b authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

renamed Leibniz to productRule and merged fLeibniz and Leibniz

parent 145edeb2
......@@ -128,25 +128,36 @@ classdef BasicVariable < handle
end
end
function h = fLeibniz(obj, phi, p)
% FLEIBNIZ compute the fractional derivative of multiplication of two functions
% f = fLeibniz(obj, phi, p)D = fLeibniz(obj, phi, p) computes the fractional leibniz
% rule of the multiplication of the function f = obj.fun and phi of fractional order
% p.
% p \~~ oo / p \ (k) p-k
% D (phi(t) f(t)) = > | | | phi (t) D f(t)
% a t /__ k=0 \ k / a t
tau = phi.domain.name;
N = 10;
h = misc.binomial( p, 0 ) * phi * obj.fDiff(p);
h0 = h;
k = 1;
while max( abs( h0 ) ) / max( abs( h ) ) > 1e-9
h0 = misc.binomial(p,k) * phi.diff(tau, k) * obj.fDiff( p - k);
h = h + h0;
k = k+1;
end
function h = productRule(obj, phi, n)
% LEIBNIZ compute the k-derivative of the product of the basic variable and phi
% h = productRule(obj, phi, n) computes the k-th derivative of obj.fun*phi, i.e.,
% h = d_t^n ( obj.fun(t) * phi )
% using the product rule. This is usefull if the derivatives of the basic variable obj
% are known exactly and the derivatives of phi, but not the multiplication
% obj.diff(0)*phi.
t = phi.domain;
h = quantity.Discrete.zeros( size(obj), t );
if numeric.near(round(n), n)
% for integer order derivatives
for k = 0:n
h = h + misc.binomial(n, k) * obj.diff(n-k) * phi.diff(t, k);
end
else
% for fractional order derivatives: [see Podlubny]
% p \~~ oo / p \ (k) p-k
% D (phi(t) f(t)) = > | | | phi (t) D f(t)
% a t /__ k=0 \ k / a t
h = misc.binomial( n, 0 ) * phi * obj.fDiff(n);
h0 = h;
k = 1;
% do the iteration on the loop until the difference is very small.
while max( abs( h0 ) ) / max( abs( h ) ) > 1e-9
h0 = misc.binomial(n,k) * phi.diff(t, k) * obj.fDiff( n - k);
h = h + h0;
k = k+1;
end
end
end
......
......@@ -175,7 +175,7 @@ function test_BasicVariable_fDiff(testCase)
end
function test_BasicVariable_fLeibniz(testCase)
function test_BasicVariable_fractional_ProductRule(testCase)
t = quantity.Domain("t", linspace(0,1));
% choose a polynomial basis
......@@ -200,7 +200,7 @@ for k = 1:3
n = ceil(alpha);
frac_f = quantity.Symbolic(fractional.fDiff(f, alpha, t.lower, t.name), t);
frac_p = b.fLeibniz(poly(k), alpha);
frac_p = b.productRule(poly(k), alpha);
testCase.verifyEqual(frac_f.on(), frac_p.on(), 'AbsTol', 4e-3)
end
......@@ -208,7 +208,21 @@ end
end
function test_BasicVariable_ProductRule(testCase)
t = quantity.Domain("t", linspace(0,1,11));
phi = quantity.Symbolic( sin( sym("t") * pi), t);
bfun = quantity.Symbolic( cos(sym("t") * pi), t);
b = signals.BasicVariable( bfun, { bfun.diff(t, 1), bfun.diff(t, 2), bfun.diff(t, 3) });
for k = 0:3
h1 = b.productRule(phi, k);
h2 = diff( phi * bfun, t, k);
testCase.verifyEqual( h1.on, h2.on, 'AbsTol', 1e-13 );
end
end
......
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