Commit 504b8caa authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Added fractional differentiation to the basic variable.

parent d508f367
......@@ -73,6 +73,72 @@ classdef BasicVariable < handle
end
end
function D = fDiff(obj, p)
% FDIFF compute the p-fractional derivative
% D = fDiff(obj, p) computes the fractional derivative of order p. To be
% concrete the Riemann-Liouville fractional derivative with initial point t0 = lower
% end of the considered time domain is computed.
% Because the basisc variable must be C^infinity function, we use the form
% of the fractional derivative described in eq. (2.80) in the book [Podlubny: Fractional
% differential equations]. Then, integration by parts is applied to eliminate the
% singularity:
%
n = ceil(p);
tDomain = obj(1).fun(1).domain;
t0 = tDomain.lower;
t = Discrete( tDomain );
tauDomain = tDomain.rename( tDomain.name + "_Hat" );
tau = Discrete( tauDomain );
if p > 0
D = quantity.Discrete.zeros(size(obj), tDomain);
for k = 0:n
f_dk = obj.diff(k);
if f_dk.at(t0) ~= 0
D = D + f_dk.at(t0) * (t - t0).^(k-p) / gamma( k-p + 1);
end
end
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 );
elseif p < 0 && p > -1
alpha = -p;
D = (t - t0).^(alpha) / gamma( alpha + 1) * obj.diff(0).at(t0) + ...
cumInt( (t-tau).^alpha * subs(obj.diff(1), tDomain.name, tauDomain.name), ...
tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha + 1 );
elseif p <= -1
alpha = -p;
D = cumInt( (t-tau).^(alpha-1) * subs( obj.diff(0), tDomain.name, tauDomain.name), ...
tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha );
elseif p == 0
D = [obj.fun()];
else
error("This should not happen")
end
end
function D = fLeibniz(obj, phi, p)
tau = phi.domain.name;
N = 10;
D = misc.binomial( p, 0 ) * phi * obj.fDiff(p);
for k = 1:N
D = D + misc.binomial(p,k) * phi.diff(tau, k) * obj.fDiff( p - k);
end
end
function d = domain(obj)
d = obj.fun(1).domain;
end
......
......@@ -127,4 +127,75 @@ testCase.verifyEqual(b(1).highestDerivative, g(1).highestDerivative)
end
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);
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.verifyEqual( b_fdiff.on(), f_fdiff.on(), 'AbsTol', 1e-2)
%% verify fractional integration
testCase.verifyEqual( b.fDiff(0).on(), F.on() )
for k = 1:5
alpha = 0.3 * k;
J = int( (t- sym("tau")).^(alpha -1) * subs(F.sym, "t", "tau"), "tau", t0, t) / gamma(alpha);
J_num = b.fDiff( -alpha );
J = quantity.Symbolic(J, timeDomain);
testCase.verifyEqual( J.on() , J_num.on, 'AbsTol', 2e-3 );
end
end
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);
end
p = signals.BasicVariable( poly(1), num2cell(poly(2:end)) );
p.fLeibniz(poly(1), 0.5)
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