From 3d80341bf20ea38d7d53390f021784aa9541e4bc Mon Sep 17 00:00:00 2001 From: Ferdinand Fischer <ferdinand.fischer@uni-ulm.de> Date: Thu, 30 Jul 2020 14:57:35 +0200 Subject: [PATCH] renamed Leibniz to productRule and merged fLeibniz and Leibniz --- +signals/BasicVariable.m | 49 +++++++++++++++++++------------ +unittests/+signals/signalsTest.m | 18 ++++++++++-- 2 files changed, 46 insertions(+), 21 deletions(-) diff --git a/+signals/BasicVariable.m b/+signals/BasicVariable.m index 290fe5d..e3f9348 100644 --- a/+signals/BasicVariable.m +++ b/+signals/BasicVariable.m @@ -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 diff --git a/+unittests/+signals/signalsTest.m b/+unittests/+signals/signalsTest.m index a70dbc3..568cee4 100644 --- a/+unittests/+signals/signalsTest.m +++ b/+unittests/+signals/signalsTest.m @@ -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 -- GitLab