derivative.m 1.27 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
function [h, t] = derivative(y, p, t0, t, optArg)
% DERIVATIVE calculates a fractional derivative
%	[h, t] = derivative(y, p, t0, t, optArg) comput the p-th fractional derivative of a symbolic
%	expression y. The lower terminal is t0 and the symbolic independent variable is t. By the
%	name-value-pair argument, the "type" of the fractional derivative can be specified. Use "RL" for
%	the Riemann-Liouville derivative and "C" for the Caputo derivative.
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)
	h = diff(y, round(p));
	return
end

n = ceil(p);
tau = sym( string(t) + "_");

if optArg.type == "C"	
	h = 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"
	h = 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