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

Test for the fractional leibniz rule for basic variables

parent 6140352b
function [y, m] = monomialFractionalDerivative(t, t0, nu, p, optArg)
%MONOMIALFRACTIONALDERIVATIVE compute fractional derivative of a monomial
% [y, m] = monomialFractionalDerivative(t, a, nu, p, optArg) compute the p-fractional derivative of
% the monomial m = (t-a)^nu.
% the monomial m = (t-a)^nu. At this, t must be a symbolic variable, t0 is the lower terminal as
% double value, nu is the exponent as double, p is the fractional order as double the type of the
% derivative can be specified by the name-value-pair "type". So far online the "RL" :
% Riemann-Liouville derivative is implemented.
arguments
t (1,1) sym; % symbolic variable for the monomial
t0 (1,1) double; % expansion point of the monome
......
classdef BasicVariable < handle
%BasicVariable Class to handle basic variables for the trajectory plannning
%BasicVariable Class to handle quantities with a priori known derivatives
% This class is used to simplify the usage of quantities for which a lot of computations must be
% done on its derivatives. So a BasicVariable consists of a basic function and its derivatives,
% which are precomputed and safed. So computations can be performed without recalculation of the
% derivatives.
properties ( SetAccess = protected )
derivatives cell;
......@@ -101,7 +105,8 @@ classdef BasicVariable < handle
end
end
assert( ~ any( isinf( obj.diff(n+1).on() ) ), "conI:BasicVariable:fDiff", ...
"The function has a singularity! Hence, the fractional derivative can not be computed")
D = D + cumInt( ( t - tau ).^(n-p) * subs(obj.diff(n+1), tDomain.name, tauDomain.name), ...
tauDomain, tauDomain.lower, tDomain.name) / gamma( n+1-p );
......@@ -123,17 +128,24 @@ classdef BasicVariable < handle
end
end
function D = fLeibniz(obj, phi, p)
% FLEIBNIZ compute the fractional derivative of multiplication of two functions:
% D^p ( phi(t) * obj.fun(t) )
% = sum_i=0^infty (p, k) * dt^k phi(t) * D^(p-k) obj.fun(t)
%
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;
D = misc.binomial( p, 0 ) * phi * obj.fDiff(p);
for k = 1:N % TODO@ferdinand: make while loop with lower threshold.
D = D + misc.binomial(p,k) * phi.diff(tau, k) * obj.fDiff( p - k);
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
end
......
......@@ -132,27 +132,33 @@ function test_BasicVariable_fDiff(testCase)
t = sym('t');
t0 = 0;
timeDomain = quantity.Domain('t', linspace(t0, t0+1, 1e2));
beta = 2;
F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
d1F = F.diff("t", 1);
d2F = F.diff("t", 2);
%% verify singularity case
beta = 1.5;
F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
b(1,:) = signals.BasicVariable(F, {F.diff("t", 1), F.diff("t", 2), F.diff("t", 3)});
p = 0.3;
nu = 2;
for k = 1:5
% compute the fractional derivative with the formula for the polynomial (see Podlubny: 2.3.4)
f_fdiff(k,1) = quantity.Symbolic( gamma(1+nu) / gamma( 1 + nu - k*p ) * ( t - t0 ).^(nu - k*p), timeDomain);
% compute the fractional derivative numerically
b_fdiff(k,1) = b.fDiff(k*p);
end
% plot(f_fdiff - b_fdiff)
testCase.verifyError( @() b.fDiff(p), "conI:BasicVariable:fDiff")
testCase.verifyEqual( b_fdiff.on(), f_fdiff.on(), 'AbsTol', 1e-2)
%% verify regular case
for beta = [1, 2, 3]
F = quantity.Symbolic(( t - t0 )^beta, timeDomain);
b(1,:) = signals.BasicVariable(F, {F.diff("t", 1), F.diff("t", 2), F.diff("t", 3)});
alpha = 0.3;
for k = 1:5
p = k*alpha;
% compute the fractional derivative with the formula for the polynomial (see Podlubny: 2.3.4)
f_fdiff(k,1) = quantity.Symbolic( fractional.monomialFractionalDerivative(t, t0, beta, p), timeDomain);
% compute the fractional derivative numerically
b_fdiff(k,1) = b.fDiff(p);
end
% plot(f_fdiff - b_fdiff)
testCase.verifyEqual( b_fdiff.on(), f_fdiff.on(), 'AbsTol', 2e-2)
end
%% verify fractional integration
testCase.verifyEqual( b.fDiff(0).on(), F.on() )
......@@ -164,7 +170,7 @@ function test_BasicVariable_fDiff(testCase)
J = quantity.Symbolic(J, timeDomain);
testCase.verifyEqual( J.on() , J_num.on, 'AbsTol', 2e-3 );
testCase.verifyEqual( J.on() , J_num.on, 'AbsTol', 3e-3 );
end
end
......@@ -174,27 +180,55 @@ function test_BasicVariable_fLeibniz(testCase)
t = quantity.Domain("t", linspace(0,1));
% choose a polynomial basis
N = 5;
for k = 1:N
poly(k) = quantity.Symbolic(sym("t").^(k-1), t);
beta = 1;
h(1,1) = quantity.Symbolic(sym("t").^beta, t);
poly(1,1) = quantity.Symbolic( sym("t"), t);
for k = 2:N
h(k,1) = h(k-1).diff();
poly(k,1) = quantity.Symbolic( sym("t")^(k-1) / factorial(k-1), t);
end
p = signals.BasicVariable( poly(2), num2cell(poly(3:end)) );
b = signals.BasicVariable( h(1), num2cell(h(2:end)) );
for k = 1:3
f = sym(p.diff(0)) * sym(poly(1));
f = sym(b.diff(0)) * sym(poly(k));
alpha = 0.3 * k;
n = ceil(alpha);
frac_f = diff( ...
int( (sym("t") - sym("tau"))^(n-alpha-1) * subs(f, "t", "tau"), "tau", 0, "t"), ...
"t", n) / gamma(n - alpha);
frac_f = quantity.Symbolic(frac_f, t);
frac_p = p.fLeibniz(poly(1), alpha);
frac_f = quantity.Symbolic(fractional.fDiff(f, alpha, t.lower, t.name), t);
frac_p = b.fLeibniz(poly(k), alpha);
testCase.verifyEqual(frac_f.on(), frac_p.on(), 'AbsTol', 4e-3)
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