分数阶微积分学是整数阶微积分学的直接拓展,将一阶导数、二阶导数、一重积分、二重积分等整数阶微积分拓展到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定义的分数阶积分在初值不为零时是与其他定义不等价的。
事实上,分数阶微积分求数值解还有更多高精度的算法,本文只给出了比较基础的代码。但是无论什么方法,都很难避免会引入新的误差,分数阶微积分的深度还亟待更加的深入。