Commit 9003d24f authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Fractional derivative (caputo) for symbolic expressions

parent a9aa58f9
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 result = nearInt(number, tolerance)
if nargin == 1
result = numeric.near(number, round(number));
else
result = numeric.near(number, round(number), tolerance);
end
end
\ No newline at end of file
function [tests] = testFractional()
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
tests = functiontests(localfunctions);
end
function testDerivative(testCase)
t = sym("t");
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
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