derivative.m 1.27 KB
 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``````