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

Expanded symbolic computation of fractional derivatives

parent 504b8caa
function [DaY, t] = fDiff(y, p, t0, t, optArg)
% FDIFF(Y, A, t0, t) calculates the Caputo fractional derivative of
% order A of a symbolic function Y.
% ??? assuming that the input was 0 for t <= t0.
% The output DAY is symbolic.
arguments
y sym
p (1,1) double {mustBePositive};
t0 (1,1) double;
t (1,1) sym;
optArg.type = "RL"
end
if p < 0
error('a must be positive! Use fintegral instead?');
end
if misc.nearInt(p)
DaY = diff(y, round(p));
return
end
n = ceil(p);
tau = sym( string(t) + "_");
if optArg.type == "C"
DaY = 1/gamma(n-p) * int((t-tau)^(n-p-1) * subs(diff(y, t, n), t, tau), tau, t0, t);
%DaY = 1/gamma(1+n-a) * ((t - t0)^(n-a) * subs( diff(y, n), t, t0) ...
% + int((t-tau)^(n-a) * subs(diff(y, 1+n), t, tau), tau, t0, t)); % Integration by parts
elseif optArg.type == "RL"
DaY = 1 / gamma(n-p) * diff( int( (t-tau)^(n-p-1) * subs(y, t, tau), tau, t0, t), t, n);
else
error("The derivative of the type " + optArg.type + " is not implemented. Only C for Caputo and RL for Riemann-Liouville are available")
end
end
\ No newline at end of file
function [DaY, t] = fdiff(y, a, t0, t)
% FDIFF(Y, A, t0, t) calculates the Caputo fractional derivative of
% order A of a symbolic function Y.
% ??? assuming that the input was 0 for t <= t0.
% The output DAY is symbolic.
arguments
y sym
a (1,1) double {mustBePositive};
t0 (1,1) double;
t (1,1) sym;
end
if a < 0
error('a must be positive! Use fintegral instead?');
end
if misc.nearInt(a)
DaY = diff(y, round(a));
return
end
m = ceil(a);
syms tau;
%output = 1/gamma(m-a) * int((t-tau)^(m-a-1) * subs(diff(input, m), t, tau), tau, t0, t);
DaY = 1/gamma(1+m-a) * ((t - t0)^(m-a) * subs( diff(y, m), t, t0) ...
+ int((t-tau)^(m-a) * subs(diff(y, 1+m), t, tau), tau, t0, t)); % Integration by parts
end
\ No newline at end of file
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.
arguments
t (1,1) sym; % symbolic variable for the monomial
t0 (1,1) double; % expansion point of the monome
nu (1,1) double; % order of the monome
p (1,1) double; % order of the fractional derivative
optArg.type string = "RL"; % Type of the fractional derivative
end
m = (t-t0)^nu;
switch optArg.type
case "RL"
y = gamma( 1+nu ) / gamma( 1+nu-p ) * ( t - t0 )^(nu-p);
otherwise
error("Not jet implemented, only the RL - Rieman Liouville derivative is available")
end
end
\ No newline at end of file
......@@ -121,19 +121,18 @@ classdef BasicVariable < handle
else
error("This should not happen")
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)
%
tau = phi.domain.name;
N = 10;
D = misc.binomial( p, 0 ) * phi * obj.fDiff(p);
for k = 1:N
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);
end
end
......
......@@ -178,24 +178,23 @@ 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)
p = signals.BasicVariable( poly(2), num2cell(poly(3:end)) );
for k = 1:3
f = sym(p.diff(0)) * sym(poly(1));
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);
end
end
......@@ -6,14 +6,31 @@ end
function testDerivative(testCase)
t = sym("t");
t0 = 0.5;
for a = [0.1 0.4 0.7 0.9]
% compute the fractional derivative with the function
datt = fractional.fdiff(t^2, a, 0, t);
% define the analytic solution for the fractional derivative
datt_ = 2 / gamma(3 - a) * t^(2 - a);
% evaluate both and compare them
T = linspace(0, 1, 11);
testCase.verifyEqual( double(subs(datt, t, T)), double(subs(datt_, t, T)), 'AbsTol', 1e-15);
end
for p = [0.1 0.4 0.7 0.9]
% define the analytic solution for the fractional derivative
[y, m] = fractional.monomialFractionalDerivative( t, t0, 1, p);
% compute the fractional derivative with the function
f = fractional.fDiff(m, p, t0, t);
% evaluate both and compare them
dt = 0.1;
% skip the first value because there is a singularity that can not be computed by the
% default symbolic evaluation.
% TODO@ferdinand: allow the quantity.Symbolic to use the limit function for the evaluation
% at certain points.
tDomain = quantity.EquidistantDomain("t", t0+dt, t0+1, "stepSize", dt);
Y = quantity.Symbolic( y, tDomain);
F = quantity.Symbolic( f, tDomain);
testCase.verifyEqual( Y.on(), F.on(), 'AbsTol', 1e-15);
h = fractional.fDiff(m, p, t0, t, 'type', 'C');
H = quantity.Symbolic( h, tDomain);
testCase.verifyEqual( H.on(), Y.on(), 'AbsTol', 1e-15);
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