MATLAB求分数阶微分的数值解，G-L定义，R-L定义，Caputo定义

        分数阶微积分学是整数阶微积分学的直接拓展，将一阶导数、二阶导数、一重积分、二重积分等整数阶微积分拓展到0.75阶导数、阶导数等实数甚至是复数阶的导数或积分。这无疑拓展了微积分学的深度。

        对于整数阶微积分，一般可以具有简洁明确的物理意义，比如位移、速度和加速度可以很好地解释一个信号与其整数阶导数之间的关系。然而分数阶微积分却没有那么简洁易懂的物理解释。目前对于分数阶微积分的定义，比较应用广泛的是G-L定义，R-L定义和Caputo定义。

Grünwald-Letnikov定义：

用MATLAB语言编写出Grünwald-Letnikov分数阶微积分数值计算的函数：

function dy = glfdiff(y,t,gam)

if strcmp(class(y),'function_handel')
    y = y(t);
end

h = t(2)-t(1);
w = 1;
y = y(:);
t = t(:);

for j = 2:length(t)
    w(j) = w(j-1)*(1-(gam+1)/(j-1));
end
for i = 1:length(t)
    dy(i) = w(1:i)*[y(i:-1:1)]/h^gam;
end

        该函数的调用格式为：y1=glfdiff(y,t,r)。其中t为等间距的时间向量，r为阶次，y可以为函数f(t)在t各点的上的函数值，也可以为f(t)的函数句柄，若r为负值，则计算原函数的-r阶积分。

Riemann-Liouville定义：

f(t)的分数阶积分定义：

f(t)的分数阶微分定义应基于积分定义给出：

 下面给出一种R-L定义的分数阶微积分的数值解法的函数：

function [dy,t] = rlfdiff(y,t,r)

h = t(2)-t(1);
n = length(t);
dy1 = zeros(1,n);
y = y(:);
t = t(:);

if r>-1,
    m = ceil(r)+1;
    p = m-r;
    y3 = t.^(p-1);
elseif r==-1
    n = length(t);
    yy = 0.5*(y(1:n-1)+y(2:n)).*diff(t);
    dy(1) = 0;
    for i=2:n
        dy(i) = dy(i-1)+yy(i-1);
    end
    return
else
    m=-r;
    y3 = t.^(m-1);
end

for i = 1:n
    dy1(i) = y(i:-1:1).'*(y3(1:i));
end

if r>-1
    dy = diff(dy1,m)/(h^(m-1))/gamma(p);
    t = t(1:n-m);
else
    dy = dy1*h/gamma(m);
end

        该函数的调用格式为[y1,t]=rlfdiff(y,t,r)。其调用格式类似刚刚的glfdiff()函数，不同的是该函数返回了两个向量y1和t。因为该函数使用了MATLAB求差分的函数diff()，向量y1的实际长度可能小于y的长度。

Caputo定义：

        事实上Caputo分数阶微积分的定义与R-L定义是相同的，对于r阶微积分，记m=[r]，则二值之间的关系为：

Caputo分数阶微积分的数值解函数：

function dy = caputo(y,t,alfa,y0,L)

dy = glfdiff(y,t,alfa);
if nargin<=4
    L = 10;
end
if alfa>0
    q = ceil(alfa);
    if alfa<=1
        y0 = y(1);
    end
    for k = 0:q-1
        dy = dy-y0(k+1)*t.^(k-alfa)./gamma(k+1-alfa);
    end
    yy1 = interp1(t(L+1:end),dy(L+1:end),t(1:L),'spline');
    dy(1:L) = yy1;
end

        该函数的调用格式为y1=caputo(y,t,r,y0,L)。当r<=0时，则可以返回Caputo分数阶微积分；如果r<1，则y(t0)由向量y直接提取出来，就可以计算Caputo分数阶导数；如果r>1，则初值向量y0实际上的内容y0=[y(t0),y'(t0),y"(t0),...]直到[r]-1阶。

举两个例子观察一下不同定义下求得的分数阶导数：

        对于阶跃函数，其0.75阶导数的表达式为：t^(-0.75)/Gama(0.25)。比较G-L定义和R-L定义的微积分的数值解与精确解：

t0 = 0:0.01:3;
y0 = t0.^(-0.75)/gamma(0.25);
h = 0.01;
t1 = 0:h:3;
y = ones(size(t1));
y1 = glfdiff(y,t1,0.75);
[y2,t2] = rlfdiff(y,t1,0.75);
plot(t0,y0,t1,y1,'--',t2,y2,':');
ylim([0 6]);

运行结果：

根据运行结果，G-L定义的分数阶微分求得的数值解更接近于其精确解，这也可能是函数求解的方法的精度的差异，并不能依此说明G-L定义的分数阶积分更准确。

        对于函数y=sin(3t+1)，对比其G-L定义和Caputo定义下的0.3阶微分的数值解：

t = 0:0.01:pi;
y = sin(3*t+1);
d = t.^(-0.3)*sin(1)/gamma(0.7);

y1=glfdiff(y,t,0.3);
y2=caputo(y,t,0.3,0);

plot(t,y1,t,y2,'--',t,d,':')
axis([0 pi -3 3])   

运行结果：

由于函数在t0时刻的初值为sin1，故两个定义间的数值解求取出现了偏差。也就是说Caputo定义的分数阶积分在初值不为零时是与其他定义不等价的。

       事实上，分数阶微积分求数值解还有更多高精度的算法，本文只给出了比较基础的代码。但是无论什么方法，都很难避免会引入新的误差，分数阶微积分的深度还亟待更加的深入。