第3章
分数阶微积分：定义与计算
正如以前指出的那样，分数阶微积分学可以追溯到 Newton和 Leibniz刚创建传统微积分学的年代。由于没有统一并被广泛接受的定义，分数阶微积分学在其早期发展过程中并没有很好的进展。直到 19世纪中叶，一些著名学者包括法国数学家 Liouville（1809－1882）于 1834年、德国数学家 Riemann（1826－1866）于 1847年、捷克数学家 Grünwald（1838－1920）于 1867年和俄国数学家 Letnikov（1837-1888）于 1868年分别提出了有意义的分数阶微积分的定义，分数阶微积分学的领域才真正被建立起来。
后来，Liouville和 Riemann的定义统一成 Riemann–Liouville定义，而 Grün-wald和 Letnikov的定义被整合成 Grünwald–Letnikov定义。这两个定义一直在分数阶微积分学领域广泛使用，后来被证明，这两个定义在工程应用中是等效的。这两个定义都适合描述初始值为零函数的分数阶微积分问题。
对于非零初值的问题，1967年，意大利数学家 Caputo提出的定义是非常有用的，特别是在分数阶微分方程的研究中，Riemann–Liouville与 Grünwald–Letnikov两个定义需要已知信号分数阶导数在初始时刻的值，而 Caputo定义需要已知信号及其整数阶导数在初始时刻的值，更接近实际应用。
除了这几个广泛应用的定义之外，文献中还有其他分数阶导数的定义，如 Erdélyi–Kober定义、Hadamard定义、Marchaud定义、Riesz定义、Riesz–Miller定义、Miller–Ross定义、Weyl定义等。本书不涉及这些分数阶微积分的定义，但可以将本书介绍的数值算法拓展到其他定义中。
本章研究的主要问题是：若原函数 f(t)或其样本点已知，如何得出该函数的分数阶导数或积分。后面章节还将介绍如果原函数未知，如何构造信号的分数阶导数与积分信号。
这里将引入一个统一的分数阶微分、积分的算子 t0 Dtα ，并广泛应用于本书的叙述，其中，α限于实数，t为自变量，t0为该变量的下边界。
·48·分数阶微积分学——数值算法与实现
说明 3-1统一的微积分算子
（1）若 α . 0且 t0 =0，则可以省略记号 t0，如果自变量为 t且没有其他变量，则 t也可以省略。
（2）若 α> 0，t0 Dtα 算子表示函数对自变量 t的 α阶导数，α =0表示原信号，若 α< 0，则表示 .α阶积分。
（3）如果 α为复数，则它的实部决定微分还是积分，这超出本书的范围。
本章 3.1节首先对整数阶 Cauchy积分公式进行拓展，给出分数阶的 Cauchy积分公式，但后续内容中不再考虑该定义。在 3.2节和 3.3节中，分别对传统微积分定义进行直接扩展，介绍分数阶微积分的 Grünwald–Letnikov定义和 Riemann– Liouville定义，并给出这些定义的数值算法，编写 MATLAB的通用求解函数。3.4节对非零初值函数的分数阶微积分概念进行拓展，给出 Caputo分数阶微积分定义。
3.5节将给出各种微积分定义之间的关系，并介绍基于 Caputo定义的分数阶微积分数值计算方法。3.6节将给出分数阶微积分的性质，并试图给出分数阶积分的物理解释方法。
必须指出的是：
（1）本章给出的分数阶微积分计算方法只是常规的算法，用于分数阶微积分的概略计算，一般精度较低。第 4章将全面介绍各种定义下函数分数阶微积分的高精度数值方法与 MATLAB实现。
（2）这里提到的数值解方法是在原函数或其样本点已知时提出的，如果原函数未知，其分数阶微积分信号的构造可以参见第 5章。
3.1分数阶Cauchy积分公式
在复分析领域有一个特别重要的 Cauchy积分公式，可以将已知复函数 f(z)的高阶导数问题转换为封闭曲线的积分问题。本节将探讨将 Cauchy积分公式拓展到分数阶导数的分数阶 Cauchy积分公式，并给出相应分数阶微积分的性质。
3.1.1 Cauchy积分公式
首先观察 Cauchy积分公式的整数阶导数表示
ddtnn f(t)= 2nπ! jC (τ . f(tτ)) n+1 dτ（3-1-1）
式中，C为光滑的封闭曲线，在该区域内函数 f(t)是单值解析函数。
如果 n被非整数 γ取代，在 t = τ点处将出现孤立奇点。在计算封闭路径积分时应该移除该奇点，并将阶乘替换成 Gamma函数，然后即可直接拓展式（ 3-1-1）处理分数阶导数的问题。
3.1.2常用函数的分数阶微分与积分公式
在传统微积分中，众所周知，三角函数的 n阶导数可以如下计算： ddtnn sin at = a n sin at + n2π （3-1-3） dn
dtn cos at = a n cos at + n2π （3-1-4）如果 n由任意实数 α取代，则上面的定义可以直接扩展，得出三角函数的分数阶微积分公式 Dtα sin at = a α sin at + α2π （3-1-5） απ
Dα t cos at = a α cos at + 2 （3-1-6）
此外，指数函数的分数阶导数可以由下式求出：
Dα t eλt = λ.αE1,1.α(λt) （3-1-7）
幂函数 tm的 n阶导数可以表示成
dn dtn tm = m(m . 1) · · · (m . n + 1) tm.n = m! (m . n)! tm.n （3-1-8）
用 Gamma函数取代阶乘，则该式可以直接扩展到分数阶导数：
Γ(m + 1)
Dtαtm = Γ(m . α + 1) tm.α , m> .1（3-1-9）
·50·分数阶微积分学——数值算法与实现
这里，m并不局限于整数，它也可以是非整数。
文献 [1]还给出了更多常用函数的分数阶导数公式，感兴趣的读者可以参考该著作。不过，该定义本身的条件也是比较苛刻的，因为该定义要求 f(t)在封闭区域内为单值解析函数。本书将不再进一步探讨与使用基于 Cauchy积分公式的分数阶微积分定义。
3.2 Grünwald–Letnikov分数阶微积分定义与计算
Grünwald–Letnikov的分数阶微积分定义是使用最广泛的分数阶微积分定义之一。本节先给出整数阶高阶导数的一般表示形式，并对其拓展，给出 Grünwald– Letnikov的分数阶导数定义。本节还给出一些基于该定义的数值求解算法，并通过例子比较各种算法的优劣。
3.2.1高阶整数阶导数的推导
在给出 Grünwald–Letnikov分数阶导数定义之前，先观察一下整数阶导数的问题。先考虑函数 f(t)的一阶导数公式
d1 f(t)= lim 1 f(t) . f(t . h)（3-2-1）
dt1 h→0 h
由上述的结果可以很容易地推导出二阶导数公式
d2 f(t)= lim 1 f(t) . 2f(t . h)+ f(t . 2h)（3-2-2）
dt2 h→0 h2
循环使用上述方法，最终可以得出函数的 n阶导数为
n
ddtnn f(t)= lim h1 n . (.1)j nj f(t . jh)（3-2-3）
h→0
j=1
式中，二项式展开式可以写成
nnn
. ..
(1 . z)n =(.1)j nj zj = j(! (.n 1). jnj!)! zj = wjzj（3-2-4）j=0 j=0 j=0
且二项式系数可以由下式直接计算：
wj =(.1)j nj =(.1)j j!(nn. ! j)! ,j =0, 1, 2, ··· ,n（3-2-5）
3.2.2 Grünwald–Letnikov分数阶微分的定义
可以容易地将上面的 n阶导数公式直接拓展到非整数 α的情形。和整数阶不同的是，二项式表达式不再是有限项的和，而变成了无穷级数的形式，即
∞∞
..
(1 . z)α =(.1)j αj zj = wjzj（3-2-6）
j=0 j=0
这样，拓展的二项式系数变成
α (.1)jΓ(α + 1)
wj =(.1)j j = Γ(j + 1)Γ(α . j + 1) ,j =0, 1, 2, ···（3-2-7）
比较式（ 3-2-5）与式（ 3-2-7）中的二项式系数可见，后者只是将前者的阶乘替换成了更一般的 Gamma函数，从而扩大了其适用范围。
假设 t . t0时函数 f(t)的值为零，则无限项的和可以转换成有限项的和，这样可以引入 Grünwald–Letnikov分数阶导数公式。
说明 3-2 Grünwald–Letnikov定义
（1）算子左上角的 GL记号表示 Grünwald–Letnikov定义，没有冲突时可略。
（2）可以看出整数阶微分只使用当前和前几个有限步长内的函数值，而分数阶微分涉及从 t0开始的所有函数值，因而可以认为分数阶导数是有记忆的。
（3）若 α为整数，式（ 3-2-8）与式（ 3-2-3）完全一致。
（4）该定义同样适用于 α> 0的微分和 α< 0的积分。另外，若 α =0，由定义
GLD0
可知 t0 t f(t)= f(t)。
（5）Grünwald–Letnikov定义满足定义 3-1。
后面将介绍各种计算二项式系数的方法，并将介绍 Grünwald–Letnikov分数阶微积分的数值计算方法。
3.2.3 Grünwald–Letnikov分数阶微分与积分的数值计算
如果选择的计算步长 h足够小，则式（ 3-2-8）中的求极限操作可以忽略，这样， Grünwald–Letnikov定义下的分数阶导数与积分可以由下面的式子直接计算：
[(t.t0)/h]
.
GL
t0 Dtαf(t) ≈ h1 α wjf(t . jh)（3-2-9）
j=0
其中，二项式系数可以由式（ 3-2-7）直接计算。这样，似乎可以很自然地给出如下的计算方法，并写出 MATLAB求解函数。
·52·分数阶微积分学——数值算法与实现
根据算法 3-1可以编写出如下的 MATLAB函数 glfdiff0()：
function df=glfdiff0(f,t,gam)
arguments, f(:,1), t(:,1) double, gam(1,1) double, end
[f,h,n]=fdiffcom(f,t); J=0:(n-1); a0=f(1);
if a0~=0 && gam>0, dy(1)=sign(a0)*Inf; end
w=gamma(gam+1)*(-1).^J./gamma(J+1)./gamma(gam-J+1);
for i=2:n, dy(i,1)=w(1:i)*f(i:-1:1)/h^gam; end
end
该函数的调用格式为 f1 =glfdiff0(f,t,α)，其中， t为等间距的时间列向量， α为阶次，这里的 f既可以为函数 f(t)在 t各点上函数值构成的列向量，也可以为 f(t)的函数句柄，换句话说，这时原始函数可以用匿名函数表示。分数阶导数的数值解在变元 f1列向量中返回。若 t向量不是列向量，函数将自动转换 t为列向量。
在该函数中使用了下一级公用函数 fdiffcom()，其内容为
function [y,h,n]=fdiffcom(f,t)
y=f; if strcmp(class(f),'function_handle'), y=f(t); end
h=t(2)-t(1); n=length(t);
end
遗憾的是，这样的算法有明显的缺陷：其中直接涉及大数的 Gamma函数值，而 Γ(172)和以后各项在 MATLAB双精度数据结构下已经成为无穷大项 Inf，所以后面各项的影响被完全忽略，导致最终计算结果出现大误差，所以分数阶导数的计算需要更可靠的方法。
上面算法的最大问题是在计算二项式系数时使用了 Gamma函数，因此导致不可避免的误差。如果能避开 Gamma函数，使用更有效的方法计算二项式系数，则可以得出更可靠的结果。下面将给出二项式系数的递推算法。
证明由式（ 3-2-7）可见， w0 =1。对二项式系数当前项比前一项，并利用性质
Γ(x+1) = xΓ(x)，就可以如下推导公式：
wj (.1)jΓ(α +1) (.1)j.1Γ(α + 1)
=
Γ(j + 1)Γ(α . j + 1) jΓ(j)Γ(α . j + 1) α . j +1 α +1
wj.1 Γ(j + 1)Γ(α . j + 1)Γ(j)Γ(α . j + 2) Γ(j)Γ(α . j + 2)
Γ(j)(α . j + 1)Γ(+ 1) j.α = . = .
= . =1 .
jj等式两端同时乘以 wj.1，就可以证明式（ 3-2-10）中的递推公式。
这样，利用式（ 3-2-10）递推求出二项式系数，就可以由式（ 3-2-9）直接计算给定函数分数阶导数。由于递推算法成功地避免了 Gamma函数的直接计算，所以解决了算法 3-1中的问题。可以证明，下面算法的精度达到 o(h) [2]。
基于上述算法，可以用 MATLAB语言编写出 Grünwald–Letnikov分数阶微积分数值计算的函数：
function dy=glfdiff(y,t,gam)
arguments, y(:,1), t(:,1) double, gam(1,1) double, end
[y,h,n]=fdiffcom(y,t); w=[1,zeros(1,n-1)]; a0=y(1); dy(1)=0;
if a0~=0 && gam>0, dy(1)=sign(a0)*Inf; end
for j=2:n, w(j)=w(j-1)*(1-(gam+1)/(j-1)); end
for i=1:length(t), dy(i,1)=w(1:i)*y(i:-1:1)/h^gam; end
end
该函数的调用格式为 y1 =glfdiff(y,t,γ)，与前面介绍的 glfdiff0()函数完全一致。和 glfdiff0()函数相比，这里只修改了 w向量的计算语句，别的没有变化。如果 γ为负值，则该函数可以直接计算原函数的 .γ阶积分。
例 3-1试求出常数的 0.75阶 Cauchy导数和 Grünwald–Letnikov导数，另外，求出这两个定义下常数的 0.75阶积分。 解在整数阶微积分中，常数的各阶导数均为零，其一阶积分是一条斜线，很自然地人们想知道，常数的分数阶微积分到底是多少。
其实常数 1可以表示成幂函数 f(t)=1= t0，这样由式（3-1-9）可见，其分数阶导数与积分分别满足
.0.75 .0.75 t0.75
D0.75 Γ(1) ttD.0.75 Γ(1) t0.75
y(t)= = ,y(t)= =
tt
Γ(1 . 0.75) Γ(0.25) Γ(1+0.75) Γ(1.75)
·54· 分数阶微积分学——数值算法与实现
由下面的语句可以得出常数的 Cauchy定义与 Grünwald–Letnikov定义的 0.75阶导数与积分，如图 3-1所示。可以看出，得出的两组曲线吻合度是很好的。
0.75阶微分
图 3-1常数的 0.75阶导数与积分
>> t=0:0.001:1; y=ones(size(t));
y1=glfdiff(y,t,0.75); y1a=t.^(-0.75)/gamma(0.25);
y2=glfdiff(y,t,-0.75); y2a=t.^(0.75)/gamma(1.75);
y2a=gamma(1)*t.^(-0.75)/gamma(1-0.75);
plot(t,y1,t,y1a,'--',t,y2,t,y2a,'--'), ylim([0 5])
选择不同的阶次 γ =0, 1/2, 3/4, 1和 5/4，函数的各阶导数如图 3-2所示。可见，当 γ等于整数 0和 1时，结果与原函数及其一阶导数完全一致。
>> t=0:0.001:0.5; y=ones(size(t)); gam=0:0.25:1.25;
for g=gam, y1=glfdiff(y,t,g); plot(t,y1), hold on; end
ylim([-3 5]), hold off
γ = 3/4
γ = 1/2
γ =1/4
γ =0
γ =1
γ = 5/4
图 3-2常数的各阶分数阶导数
函数 glfdiff()还可以用于直接计算分数阶积分，如果积分阶次分别设置为 0, 1/4, 1/2, 3/4, 1和 5/4，即设置阶次向量为 γ = [0, .1/4, .1/2, .3/4, .1, .5/4]，则可以由下面的语句绘制出分数阶积分，如图 3-3所示。可见，当 γ =0和 γ = .1时，得出的结果与原函数和整数阶积分完全一致。
>> t=0:0.01:1.5; y=ones(size(t)); gam=-1.25:0.25:0;
for g=gam, y1=glfdiff(y,t,g); plot(t,y1), hold on; end
图 3-3不同阶次的分数阶积分曲线
例 3-2选择不同的计算步长求取阶跃函数的 0.5阶导数并评价解的精度。 解选择计算步长为 h =0.0001，0.001，0.01和 0.1，就可以用下面的语句求取原函数的 0.5阶导数，并和理论值对比得出其误差，结果在表 3-1中列出。可见，理论上的 o(h)精度是正确的。
>> t0=0.2:0.2:1; y0=gamma(1)*t0(:).^(-0.5)/gamma(1-0.5);
t=0:0.0001:1; y=ones(size(t)); y1=glfdiff(y,t,0.5);
t=0:0.001:1; y=ones(size(t)); y2=glfdiff(y,t,0.5);
t=0:0.01:1; y=ones(size(t)); y3=glfdiff(y,t,0.5);
t=0:0.1:1; y=ones(size(t)); y4=glfdiff(y,t,0.5);
y1=y1(2001:2000:10001); y2=y2(201:200:1001);
y3=y3(21:20:101); y4=y4(3:2:11);
T=[[t0, t0]; [y1 y2 y3 y4]', [y0-y1 y0-y2 y0-y3 y0-y4]']
表 3-1不同计算步长下的分数阶导数的误差
计算步长 tk时刻的导数 0.2 tk时刻的误差 0.4 0.6 0.8
h 0.2 0.4 0.6 0.8
y0 1.2616 0.8921 0.7284 0.6308
0.0001 1.2615 0.8920 0.7284 0.6308 7.9×10.5 2.8×10.5 1.5×10.5 7.1×10.6
0.001 1.2608 0.8918 0.7282 0.6307 0.000786 0.000279 0.000152 7.1×10.5
0.01 1.2537 0.8893 0.7269 0.6298 0.0079 0.00278 0.00152 0.00070
0.1 1.1859 0.8647 0.7134 0.6210 0.07571 0.02738 0.015 0.00701
·56· 分数阶微积分学——数值算法与实现
在前面的例子中，因为分数阶导数的解析表达式已知，所以可以和理论值相比评价计算误差。而在一般应用中理论值是未知的，这时可以选择某计算步长得出结果，再选择一个更小的步长重新计算，看看二者是否一致。如果一致则说明结果是可靠的；若不一致则应该再减小计算步长重新计算，直到得出一致的结果。 例 3-3现在考虑给定的函数 f(t)= e.t sin(3t + 1)，t ∈ (0, π)，可以看出，该函数的初始值不是零，试计算该信号的分数阶导数。 解分别选择计算步长为 h =0.01，h =0.001与 h =0.0005，就可以由下面的语句得出已知信号的 0.5阶 Grünwald–Letnikov导数，如图 3-4所示。可以看出，除了在 t =0时刻附近的一个极小的区域外，几条曲线吻合度很高，说明 h =0.01计算步长下的结果比较精确，h =0.001与 h =0.0005的结果几乎完全重合。
>> f=@(t)exp(-t).*sin(3*t+1); %原函数还可以由匿名函数表示 t=0:0.001:pi; dy=glfdiff(f,t,0.5); %在不同步长下求导 t1=0:0.01:pi; dy1=glfdiff(f,t1,0.5); t2=0:0.0005:pi; dy2=glfdiff(f,t2,0.5); plot(t,dy,t1,dy1,'--',t2,dy2,':'), axis([0,pi,-1.5 6])
图 3-4不同计算步长下结果的比较
选择不同的微分阶次 γ =0，1/4，1/2，3/4，1，5/4，就可以由下面语句计算出该函数的分数阶微分，如图 3-5所示。
>> t=0:0.01:pi; y=exp(-t).*sin(3*t+1); gam=0:0.25:1.25;
for g=gam, y1=glfdiff(y,t,g); plot(t,y1); hold on; end
axis([0 pi, -3 3]), hold off
下面的语句还可以直接得出不同阶次下的分数阶积分，如图 3-6所示。可以看出，该函数的分数阶微分积分比整数阶微积分含有更丰富的信息。
>> for gam=-1:0.25:0
y1=glfdiff(y,t,gam); plot(t,y1); hold on
end
xlim([0 pi]), hold off
图 3-5不同阶次的分数阶导数曲线
γ =0
γ = .1
图 3-6不同阶次的积分曲线
例 3-4试用 Cauchy定义和 Grünwald–Letnikov定义分别求 f(t)= sin(3t + 1)函数的 0.75阶导数并比较其结果。
解由式（3-1-5）给出的 Cauchy积分公式可知，函数的 0.75阶导数为 0 Dt 0 .75 f(t)=
30 .75 sin (3t+1+0.75π/2)，而 Grünwald–Letnikov定义的数值解可以由 glfdiff()函数求出，二者的曲线如图 3-7所示。
>> t=0:0.01:pi;
y=sin(3*t+1); y1=3^0.75*sin(3*t+1+0.75*pi/2);
y2=glfdiff(y,t,0.75); plot(t,y1,t,y2,'--')
axis([0 pi -4 8])
比较这两条导数曲线可见，在 Cauchy定义下， t =0时并没有突然的跳变，这是因为函数 y(t)在 t . 0时也被认为满足 f(t)= sin(3t + 1)，而在 Grünwald–Letnikov定义中，假设 t< 0时的初始值为零。
·58·分数阶微积分学——数值算法与实现
图 3-7不同定义下的分数阶导数曲线
3.2.4 Podlubny的矩阵算法
Podlubny教授提出的矩阵算法也可以用来求给定函数的 Grünwald–Letnikov分数阶导数与积分[3]，其具体算法如下：
注意，式（ 3-2-11）中的矩阵其实是一个旋转的 Hankel矩阵，这样，在 MATLAB下可以实现矩阵算法，求出给定函数的分数阶导数。
function dy=glfdiff_mat(y,t,gam)
arguments, y(:,1), t(:,1) double, gam(1,1) double, end
[y,h,n]=fdiffcom(y,t); w=[1,zeros(1,n-1)];
for j=2:n, w(j)=w(j-1)*(1-(gam+1)/(j-1)); end
dy=rot90(hankel(w(end:-1:1)))*y/h^gam;
end
可以看出，这里给出的方法在数学上更简洁，也更适合于信号的序贯求导，如求解 DαDβDγy(t)的问题可以转换成矩阵乘积进行计算。该算法的局限性是，如果总的计算点数 N过大，则可能需要连 MATLAB都不能处理的大型矩阵运算。
3.2.5短时记忆效应及其探讨
一般情况下，如果计算步长 h选择得过小，或 [(t.t0)/h]的值过大，则式（ 3-2-9）参与求和的点数会特别庞大，可能导致计算量显著增加，最终得不出计算结果，这样应该考虑减少计算点数。在实际应用中计算分数阶导数不一定非得使用以前 [(t . t0)/h]所有的信息，用最近时间区间 [t . L, t]内的信息就能减少计算量[2]：
t0 Dtαf(t) ≈ (t.L)Dtαf(t)（3-2-12）
这种方法又称为短时记忆效应（short-memory effect）[2，4]，采用该方法，Grünwald– Letnikov分数阶导数可以近似为
N.(t)
y(t) ≈ h1 α wjf(t . jh)（3-2-13）
j=0
式中， N(t)= min t . t0 ,L （3-2-14）
hh
其中， L称为记忆时长。这样，如何选择这个记忆时长 L将成为最关键的问题。假设感兴趣时间区间 (t0,T )内函数 f(t)的值满足 |f(t)| . M，则可见近似误差为
.(t)= Dαf(t) . (t.L)Dαf(t)
t0 tt
=1 L f ′ (τ) dτ . ML.α （3-2-15）Γ(1 . α)(t . τ)α |Γ(1 . α)|
t0
如果期望误差小于预先指定的正数 ε，即 .(t) <ε，则记忆时长应该选为
1/α
L . ε|Γ(1 M . α)| （3-2-16）
不过，该公式在选择记忆时长时可能过于苛刻。比如，若选择误差容限 ε = 10.6，记忆时长将是个非常大的数值，除非降低误差要求。例如，如果选择一个相对较大的误差容限 ε = 10.3，阶次 α =0.5，函数的上界为 M =1，则需要的记忆时长为 L = 318309，若选择 ε = 10.4，则记忆时长 L = 31830988，如果阶次更小，如 α =0.1，则记忆时长将是天文数字，所以在精度要求较高时，式（ 3-2-16）可能过于保守，有时短时记忆效应未必有实际意义。
基于短时记忆的思想，作者改写了 glfdiff()函数，增加了短时记忆功能，编写了新的 glfdiff_mem()函数，其调用格式为 y1 =glfdiff_mem(y,t,γ,L0)，其中，y，t，γ，y1与 glfdiff()完全一致，增加的 L0为计算点数，即 L0 = L/h。如果不给出 L0，则得出的是 Grünwald–Letnikov的分数阶微积分。
function dy=glfdiff_mem(y,t,gam,L0)
·60· 分数阶微积分学——数值算法与实现
arguments, y(:,1), t(:,1), gam(1,1), L0=length(t); end
[y,h,n]=fdiffcom(y,t); w=[1,zeros(1,n-1)]; dy(1)=0;
for j=2:L0+1, w(j)=w(j-1)*(1-(gam+1)/(j-1)); end
for i=1:n
L=min([i,L0]); dy(i,1)=w(1:L)*y(i:-1:i-L+1)/h^gam;
end, end
也可以考虑将矩阵算法 3-3中的 w向量后面的值设置为零，但这对提高矩阵运算的速度没有帮助。后面将给出例子探讨短时记忆效应的实际应用与效果。
例 3-5在分数阶导数计算公式中，二项式系数 wj表明导数值对以往函数值的依赖程度，试绘制 α =0.1与 α =0.8时的 wj曲线。 解利用式（3-2-10）中的递推公式可以直接计算出不同 α值下的 wj系数，从而绘制出曲线，如图 3-8所示，可以看出，分数阶导数对近处几个点的依赖程度很高，对远处的点依赖程度虽然减小，但也有些依赖。例如，当 α =0.1时，w100 = .5.9076×10.4。所以不能完全忽略以前点的信息。
>> t=0:0.1:3; a1=0.1; a2=0.8;
w1=1; for j=2:length(t), w1(j)=w1(j-1)*(1-(a1+1)/(j-1)); end
w2=1; for j=2:length(t), w2(j)=w2(j-1)*(1-(a2+1)/(j-1)); end
plot(t,w1,t,w2,'--')
图 3-8不同阶次 α下的二项式系数 wj曲线
例 3-6再重新考虑例 3-3中的函数 f(t)= e.t sin(3t + 1)。试利用短时记忆效应求
取函数的 0.5阶微分，并探讨短时记忆效应的近似效果。 解选择计算步长 h =0.001，并选择容许误差为 ε = 10.3。因为函数的最大值为 M =1，由式（3-2-16）计算出来的最小记忆时长为 L = 1/(10.3Γ(1 . 0.5)) 1/0.5 =
3.2×105。当然，这样的要求有些苛刻。
选择短一些的记忆时长，探讨近似的效果，对比结果在表 3-2中给出。可以看出，如果记忆时长选择为 L0 = 1200，基本上可以得出满意的近似结果。因为全时长总共 3142
个点，选择记忆时长 L0 = 1200将减少超过 2/3的计算时间，得出的曲线如图 3-9所示。
表 3-2当计算步长为 h =0.001时记忆时长与误差的关系
记忆时长 L0 误差范数 记忆时长 L0 误差范数 记忆时长 L0 误差范数
100 0.71736 1100 0.063559 2100 0.028305
200 0.38972 1200 0.057211 2200 0.026644
300 0.26417 1300 0.051879 2300 0.025141
400 0.19737 1400 0.047346 2400 0.023777
500 0.15597 1500 0.043449 2500 0.022532
600 0.12789 1600 0.04007 2600 0.021099
700 0.10767 1700 0.037115 2700 0.018628
800 0.092463 1800 0.034513 2800 0.015194
900 0.080647 1900 0.032206 2900 0.011026
1000 0.071227 2000 0.030149 3000 0.006443
>> t=0:0.001:pi; y=exp(-t).*sin(3*t+1);
y0=glfdiff(y,t,0.5); L0=[100:100:3000]; T=[];
for L=L0
dy=glfdiff_mem(y,t,0.5,L); T=[T; L, norm(y0-dy,inf)];
end
dy=glfdiff_mem(y,t,0.5,1200);
plot(t,y0,t,dy,'--'); axis([0 pi -1 5])
图 3-9 L0 = 1200时的短时记忆效果
在同一个例子中，如果终止仿真时间变成 50s，则用普通 Grünwald–Letnikov算法将需要 2.74s才能得出分数阶导数，而使用记忆时长为 L0 = 1200的短时记忆算法则只需 0.21s。在实际应用中，在特别耗时的应用中可以考虑采用短时记忆算法，但得出的结果需要检验。
>> t=0:0.001:50; y=exp(-t).*sin(3*t+1);
tic, y0=glfdiff(y,t,0.5); toc
·62· 分数阶微积分学——数值算法与实现
tic, dy=glfdiff_mem(y,t,0.5,1200); toc
下面给出一个短时记忆效应的反例。
例 3-7试尝试用短时记忆方法求出阶跃函数的分数阶积分。
解阶跃信号的 0.75阶积分可以由 Cauchy积分公式直接得出 t0.75/Γ(1 + 0.75)。选择计算步长 h =0.001，则总共需要 5s，产生 5000个数据点。现在尝试一个较大的记忆时长，如 L0 = 4000，则可以采用下面的语句直接计算分数阶积分，如图 3-10所示。可以看出，即使采用这么大的记忆时长，还会产生很大的计算误差，说明在这个例子中短时记忆算法是不适合使用的，必须使用全部信息。
>> t=0:0.001:10; y=ones(size(t));
y0=gamma(1)*t.^(0.75)/gamma(1+0.75); % Cauchy积分公式
dy=glfdiff_mem(y,t,-0.75,4000); %短时效应
dy1=glfdiff(y,t,-0.75); plot(t,dy,t,dy1,'--');
图 3-10失效的短时记忆效应
从例 3-7可以看出，尽管很多人认为短时记忆效应是一种高效的近似方法，事实上并不然。因为使用该方法还是有风险的，有时甚至会得出错误或误导性的结果。由于式（ 3-2-16）过于苛刻，要保证不太苛刻的精度要求，也经常出现计算出的 L远远大于实际计算点数的情况。所以除非特别必要，请勿采用该算法，而应该考虑更好的计算方法。
3.3 Riemann–Liouville分数阶微积分定义与计算
前面给出的 Grünwald–Letnikov定义虽然利于数值计算，但在理论推导或公式证明方面并不容易处理，所以也应该考虑其他数学定义。例如本节要介绍的由整数阶积分拓展而来的 Riemann–Liouville定义等。
3.3.1高阶整数阶积分公式
先考虑整数阶积分。很显然，给定函数 f(t)的一阶积分可以表示为
t
ddt..11 f(t)= f(τ)dτ（3-3-1）
0
从该结果再求一次积分，就可以得出原函数的二阶积分：
tτ t
ddt..22 f(t)= f(τ)dτdt = f(τ)(t . τ)dτ（3-3-2）
00 0
类似地，可以推导出 n阶积分公式为
tτ τ
d.n
f(t)= ··· f(τ) dτ1dτ2 ··· dτ
dt.n
00 0
n重 n重（3-3-3） tt
11 f(τ)
= f(τ)(t . τ)n.1dτ = dτ
(n . 1)! (n . 1)! (t . τ)1.n
00
3.3.2 Riemann–Liouville分数阶微积分定义
如果整数阶次 n由实数 α取代，则可以直接给出如下的 Riemann–Liouville分数阶积分公式。
Riemann–Liouville定义是分数阶微积分领域应用较广的定义之一，算子 D两端的下标为积分的上下限[5]。
分数阶微分不能将积分中的 α替换成 .α直接定义，而是应该基于积分定义给出，现在考虑 β阶分数阶微分问题。如果 n . 1 <β . n，记 n = .β.，则
RL β RL .(n.β)
t0 Dt f(t)= ddtnn t0 Dt f(t)（3-3-5）基于积分定义，可以给出正式的 Riemann–Liouville分数阶导数的定义。
·64·分数阶微积分学——数值算法与实现
3.3.3常用函数的 Riemann–Liouville微积分公式
首先考虑下限 t0 =0时幂函数和指数函数的 Riemann–Liouville分数阶导数的问题。
.
.
证明指数函数 eλt的 Taylor级数展开为
∞∞
eλt =(λt)k = λktk （3-3-9）
k! Γ(k + 1)
k=0 k=0
利用式（ 3-3-7）的性质，对式（ 3-3-9）右侧的项逐项计算 Riemann–Liouville分数阶微积分，即可以证明本定理：
.
∞
RL eλt RL
0 Dα = 0 Dα
tt
k=0
.
∞
Γ(k + 1) Γ(
Γ(k +1 . α)
λktk
=
k=0
.
∞
Γ(k +1 . α)
k=0
(λt)k = t.αE1,1.α(λt)（3-3-10）
.α
= t
定理 3-4 .常用函数的 Riemann–Liouville分数阶导数
一些常用函数的 Riemann–Liouville分数阶导数公式[2]为
(t . a).α RLDtαH(t . a) = Γ(1 . α) H(t . a)（3-3-11）
RL 0 DαH(t . a)f(t)= H(t . a) RLDαf(t)（3-3-12）
t at
t.α.1
RL 0 Dtαδ(t)= Γ(.α)（3-3-13）
t.α.n.1
RL 0 Dtαδ(n)(t)= Γ(.α . n)（3-3-14）
√
RL 0 Dtαcosh λt = t.αE2,1.α λt2（3-3-15）
在上面的定理中，H(t . a)为 Heaviside函数，即阶跃函数，其含义为
H(t . a)=1,t . a （3-3-21）0,其他
函数 δ(t)为冲激（impulsive）函数，又称 Dirac函数，其定义为
∞
t .δ(t)dt =1
=0时 δ(t)=0.，且 （3-3-22）
.∞
Dirac函数是 Heaviside函数的一阶导数。
此外，引入中间变量 r = λ2 + μ2，tan φ = μ/λ，还可以得出
RL
.∞Dtαeλt sin μt = r αeλt sin (μt + αφ)（3-3-23）
RL
.∞Dtαeλt cos μt = r αeλt cos (μt + αφ)（3-3-24）
必须指出的是，目前只有很少几类函数的 Riemann–Liouville分数阶导数的解析解存在。
从 Riemann–Liouville定义来看，因为变量 t在积分边界与被积函数中同时出现，所以一般函数 f(t)的分数阶导数的解析解很难求出，只能使用数值方法求解 Riemann–Liouville分数阶导数，当然，数值解的效果也不佳。
3.3.4初始时刻平移的性质
在前面分数阶导数的叙述中，一些公式涉及的初始时刻 t0 =0。如果想把这些公式的初始时刻都一般化为 t0，则应该对感兴趣区间做出必要的变换。本节将引入一个定理，允许读者做相应的变换。
重新回顾一下 Riemann–Liouville分数阶导数的定义：
t
RL 1 dn f(τ)
t0 Dtαf(t)= Γ(n . α) dtnt0 (t . τ)1+α.n dτ（3-3-25）
·66· 分数阶微积分学——数值算法与实现
证明令 τ = θ + t0，则 dτ = dθ，θ ∈ [0,t . t0]，代入式（ 3-3-25）的右侧，则
t t.t0
1 dn f(τ)1 dn f(θ + t0)
dτ = dθ
1+α.n 1+α.n
Γ(n . α) dtnt0 (t . τ)Γ(n . α) dtn 0 (t . t0 . θ)（3-3-27）
再令 s = τ + t0，则 dt = ds，故 t . t0 = s，代入式（ 3-3-27）的右侧，则
t.t0 s
1 dn f(θ + t0)1 dn f(θ + t0)
dθ = dθ
1+α.n 1+α.n
Γ(n . α) dtn 0 (t . t0 . θ)Γ(n . α) dsn 0 (s . θ)（3-3-28）
令 F (θ)= f(θ + t0)，对式（ 3-3-28）进行变换，可以直接证明定理。
s
1 dn F (θ) RL
Γ(n . α) dsn 0 (s . θ)1+α.n dθ = 0 DsαF (s)= RL 0 Dsαf(s + t0)（3-3-29）
例 3-8已知指数函数的分数阶导数为 RL 0 Dtαeλt = t.αE1,1.α(λt)，试推导下限为 t0的 Riemann–Liouville分数阶导数。 解根据定理 3-5，记 s = t . t0，经过变量替换可以得出
RLDαeλt RL 0 Dαeλ(s+t0) = eλt0 RL 0 Dαeλs = eλt0 st0 t = ss .αE1,1.α λs （3-3-30）
= eλt0 (t . t0).αE1,1.α
λ(t . t0)
假设 λ = .1，t0 =0.5，α =0.7，可以由下面的语句直接计算 Grünwald–Letnikov导数，并利用式（3-3-30）计算精确解，如图 3-11所示。可以看出，这样绘制的两条曲线完全重合，说明由平移公式导出的解析式是正确的。
>> lam=-1; t0=0.5; alfa=0.7;
t=t0:0.01:5; u=exp(lam*t); y1=glfdiff(u,t,alfa);
y2=exp(lam*t0)*(t-t0).^(-alfa).*ml_func([1,1-alfa],lam*(t-t0));
plot(t,y1,t,y2,'--')
3.3.5 Riemann–Liouville定义的数值计算
这里给出一种 Riemann–Liouville分数阶微积分的数值计算方法。如果计算步长为 h，则可以由下面的步骤计算函数的 α阶微积分。
图 3-11函数的 Grünwald–Letnikov导数与平移公式比较
上述算法的 MATLAB实现由下面给出。事实上在算法实现时作者加入了其他内容，使得分数阶与整数阶积分也可以直接计算出来。
function [dy,t]=rlfdiff(y,t,r)
arguments, y(:,1), t(:,1) double, r(1,1) double, end
[y,h,n]=fdiffcom(y,t); dy=zeros(n,1); dy1=dy;
if r>-1, m=ceil(r)+1; p=m-r; y3=t.^(p-1);
elseif r==-1, yy=0.5*(y(1:n-1)+y(2:n)).*diff(t);
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
end
该函数的调用格式为 [y1,t]=rlfdiff(y,t,γ)，其格式类似于 glfdiff()函数，不同的是该函数返回了两个向量 y1和 t。因为该函数内部使用了 MATLAB求差分的函数 diff()，向量 y1的实际长度可能小于向量 y的长度。
例 3-9已知阶跃函数 0.75阶导数的解析解为 Dt 0.75y(t)= t.0.75/Γ(0.25)，试比较
·68· 分数阶微积分学——数值算法与实现
Grünwald–Letnikov与 Riemann–Liouville分数阶导数数值解的精度。 解选择计算步长 h =0.01，下面的语句可以计算出阶跃函数 0.75阶导数的理论值和 Grünwald–Letnikov与 Riemann–Liouville分数阶导数的数值解，如图 3-12所示。可以看出，直接计算的 Riemann–Liouville分数阶导数有较大误差。在实际应用中不建议
采用该函数计算 Riemann–Liouville分数阶微积分。
图 3-12阶跃函数的 0.75阶 Riemann–Liouville导数
>> 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])
3.3.6 Riemann–Liouville微积分的符号计算
利用 MATLAB的符号运算功能，可以直接实现定义 3-4和定义 3-5的 Riemann– Liouville分数阶积分与微分计算。可以编写如下的解析计算函数：
function dy=riemannsym(f,alpha,t0,t)
arguments, f(1,1), alpha(1,1), t0(1,1)=0, t(1,1)=symvar(f), end
syms tau; alpha=sym(alpha);
switch class(f)
case 'sym', F=subs(f,t,tau);
case 'symfun', F=f(tau);
otherwise, error('f must be a symbolic expression')
end
if alpha<0, alpha=-alpha;
dy=1/gamma(alpha)*int(F/(t-tau)^(1-alpha),tau,t0,t);
elseif alpha==0, dy=f;
else, n=ceil(alpha);
f1=int(F/(t-tau)^(1+alpha-n),tau,t0,t);
dy=diff(f1,t,n)/gamma(n-alpha);
end, dy=simplify(dy);
end
其中， f为符号表达式或符号函数， α为微分的阶次，若其值为负则求 Riemann– Liouville积分。 t0为积分下限，而 t为符号变量。由此可以计算某些函数 Riemann– Liouville微积分的解析解。
3.4 Caputo分数阶微积分定义
在前面介绍的 Grünwald–Letnikov和 Riemann–Liouville分数阶微积分定义中，如果初始条件非零，有时微分方程描述可能出现困难，所以在非零初值系统的研究与仿真中，需要引入另一种分数阶微积分的定义—— Caputo定义[6]。
3.4.1 Caputo微积分定义
可以看出，Caputo分数阶积分的定义与 Riemann–Liouville定义完全一致。
3.4.2常用的 Caputo导数公式
作为例子，先考虑一下指数函数 f(t)= eλt的 Caputo分数阶导数问题。
证明因为 α> 0，m = .α.，γ = m . α，且指数函数的 Riemann–Liouville分
·70·分数阶微积分学——数值算法与实现
数阶积分可以由定理 3-3得出，该定理重新写成 RL 0 Dt .γeλt = tγE1,1+γ(λt)（3-4-4）由于 Caputo与 Riemann–Liouville分数阶积分的定义是完全一致的，经过下面的推导可以发现，指数函数的 Caputo分数阶导数满足定理的结论。 dm
C eλt RL .γ eλt = λm RL .γeλt
0 Dα = 0 D0 D= λmtγE1,1+γ(λt)
tt t
dtm
必须指出的是，上面第一步使用的性质是后面将介绍的定理 3-18。
推论 3-1三角函数函数的分数阶导数公式还可以写成
C0Dtα sin λt =(jλ2)j m tγ[E1,γ+1(jλt) . (.1)mE1,γ+1(.jλt)]（3-4-9）
C0Dtα cos λt =(jλ2)m tγ[E1,γ+1(jλt)+(.1)mE1,γ+1(.jλt)]（3-4-10） 证明由著名的 Euler公式就可以直接证明该推论。 sin λt = 21 jejλt . e.jλt , cos λt = 12 ejλt + e.jλt（3-4-11）
令 m =0，γ = .α，可以直接用上述公式计算常用函数的 Caputo分数阶积分。后面将介绍 Caputo分数阶微积分的数值计算。
MATLAB可以编写如下的符号运算函数： symstau;alpha=sym(alpha);m=ceil(alpha); Leffler函数等特殊函数，所以有时得出的结果与本书的不完全一致。这种情况下，可以通过绘制曲线对比，验证结果的正确性。C3-100504 例 试利用符号运算计算和。(2+1)(2+1)DD..tt00tt caputosym() 解当前版本的函数在多项式函数求导方面只能处理整数阶多项式，positive,f(t)=2*t+1; %f=2*t+1>>t 或由命令定义函数syms 06得出的结果为。R=4 /R=10/(3Γ(06))π.tt，12 .C 3-643 解对这里给出的问题，由定理 可知， 故的解.5γ 07λ 3f()D.tm===，，，.0t ECaputo.caputosym()07析解为。现在尝试用函数求其导数：..243(3).tt117,.positive,f=exp(-3*t);F(t)=caputosym(f,sym(4.3),0,t)>>tsyms
3.4.3 Caputo定义的符号运算
C
事实上，有些函数的 Caputo导数是可以通过计算机解析计算的。由定义 3-6，
function dy=caputosym(f,alpha,t0,t)
arguments, f(1,1), alpha(1,1), t0(1,1)=0, t(1,1)=symvar(f), end
if alpha<=0, dy=riemannsym(f,alpha,t0,t); return; end
switch class(f)
case 'sym', F=diff(subs(f,t,tau),m);
case 'symfun', F=diff(f(tau),m);
otherwise, error('f must be a symbolic expression')
end
C
dy=1/gamma(m-alpha)*int(F/(t-tau)^(1+alpha-m),tau,t0,t);
dy=simplify(dy);
end
在函数内部先识别给定函数是符号变量还是符号函数，依据情况分别处理。
必须指出的是，由于这个函数底层使用了 MuPAD提供的符号运算，很多涉及
分数阶的特定问题无法通过解析积分得出；另外，由于 MuPAD内核没有 Mittag-
不能处理伪多项式。
R1=caputosym(f,sym(0.5),0,t), R2=caputosym(f,sym(0.4),0,t)
例 3-11已知指数函数 f(t)= e3t，试求 0Dt 4.3f(t)。
得出的结果极其烦琐，可以通过手工整理，得
√
81 10 √√
e.3t
F (t) = 27 Γ(0.7) . Γ(0.7, .3t) (1+ 5)i +10 . 25
4 Γ(0.7)
可见，这里的结果比前面推出的解析解复杂得多。不过，通过绘制 t ∈ (0, 5)区间的
曲线，可见二者绘制的曲线完全重合，由此表明，符号运算得出的结果也是正确的。
·72· 分数阶微积分学——数值算法与实现
>> t0=[0:0.01:5]'; y0=F(t0); y0=double(y0);
y1=-243*t0.^(0.7).*ml_func([1,1.7],-3*t0);
plot(t0,real(y0),t0,y1,'--') %略去微小的虚部元素
3.5各种不同分数阶微积分定义之间的关系
前面介绍了几种不同的分数阶微积分的定义，除了 Cauchy积分公式之外，其他几个定义之间是有密切关系的。在某些条件下，有些定义是相互等效的，有些定义是可以在某些条件下相互转换的。本节将介绍不同定义之间的相互关系，并介绍 Caputo定义的数值计算方法。
3.5.1 Grünwald–Letnikov与Riemann–Liouville定义的关系
文献 [2]指出，对大量的函数而言， Grünwald–Letnikov定义与 Riemann–Liou-ville定义是等效的，这里具体给出二者之间的关系。
证明定理的证明可以参见文献 [2]。对式（ 3-5-1）右端求 m阶导数，可以得出
m.1
dm . t f(m)(τ)dτ
f(j)(t0)(t . t0)m+j.α 1
GLDαf(t)= +
t0 t
(t . τ)α.2m+1
dtm Γ(1+m+j.α) Γ(2m . α)
t0
j=0
（3-5-2）若 m . 1，则上式的第一项最高为 tm.1，故求 m阶导数时该项将完全消失。使用分部积分法对上式求 m次积分，则
t
GL 1 f(m)(τ)dτ RL
t0 Dtαf(t)= Γ(m . α) t0 (t . t0)1+α.m = t0 Dtαf(t)（3-5-3）
可见，该结果与 Riemann–Liouville定义完全一致。
3.5.2 Caputo与 Riemann–Liouville定义的关系
当 0 <α< 1，可知 m =1，这样，式（ 3-4-1）可以改写成
t
C
t0 Dtα y(t) = Γ(11 . α)(ty . ′ (ττ))α dτ（3-5-4）
t0
如果函数 y(t)的初值非零，比较 Caputo与 Riemann–Liouville分数阶导数定
义可得 C Dα y(t)= RLDα y(t) . y(t0)（3-5-5）
t0 tt0 t
式中，常数 y(t0)的分数阶导数可以写成 y(t0)
RL
t0 Dtα y(t0) = Γ(1 . α)(t . t0).α（3-5-6）可以看出，Caputo与 Riemann–Liouville分数阶微积分的区别为 C Dtα y(t)= RLDtα y(t) . Γ(1 y(. t0) α)(t . t0).α（3-5-7）
t0 t0
更一般地，如果阶次 α满足 α> 0，则可以通过下面的定理描述 Caputo与 Riemann–Liouville定义之间的关系。
说明 3-3 Caputo分数阶导数公式
（1）可以看出， 0 . α . 1只是上面关系的一个特例。
（2）正如前面介绍的，如果 α< 0，则 Riemann–Liouville分数阶积分与 Caputo分数阶积分是完全等效的，在实际应用中不区分二者。
3.5.3 Caputo分数阶微分的数值计算
从 3.5.1节的叙述可知， Caputo与 Riemann–Liouville或 Grünwald–Letnikov定义的区别在于对初始值的处理，二者之差可以由式（ 3-5-8）直接计算。这样，以函数 glfdiff()为基础就可以写出 Caputo分数阶导数 C t0 Dtαy(t)计算的代码。
function dy=caputo(y,t,gam,y0,L)
arguments, y(:,1), t(:,1), gam(1,1), y0, L(1,1)=10, end
dy=glfdiff(y,t,gam);
if gam>0, q=ceil(gam); if gam<=1, y0=y(1); end
for k=0:q-1, dy=dy-y0(k+1)*t.^(k-gam)./gamma(k+1-gam); end
·74· 分数阶微积分学——数值算法与实现
yy1=interp1(t(L+1:end),dy(L+1:end),t(1:L),'spline'); dy(1:L)=yy1; end, end该函数的调用格式为 y1 =caputo(y,t,α,y0,L)，当 α . 0，则可以返回 Grün-wald–Letnikov分数阶积分，即 Caputo分数阶积分；如果 α< 1，则 y(t0)由向量 y直接提取出来，就可以计算 Caputo分数阶导数；如果 α> 1，则初值向量 y0实际上的内容为 y0 = y(t0),y ′ (t0), ··· ,y(q.1)(t0)，q = .α.。
说明 3-4 Caputo分数阶微积分
（1）在实际计算中，返回的 Caputo导数前几个初始点会有较大的误差，这样采用样条插值的方法重新计算前 L个点的值， L的默认值为 10，且如果增大分数阶的阶次，则 L的值也应该增加。
（2）因为在前几个点中使用样条插值计算，所以有时得出的结果误差很大，这个函数并不是很实用，后面将给出更好的计算方法及其实现，而不建议采用这里给出的算法与函数。
例 3-12考虑例 3-4中的正弦函数的分数阶微分问题，f(t)= sin(3t + 1)。试求出其
0.3，1.3与 2.3阶 Caputo分数阶导数。 解可以看出，函数在 t =0时刻的初值为 sin 1，这样，两个定义之间的区别为 d(t)= t.0.3 sin 1/Γ(0.7)。通过下面的语句可以绘制出 Grünwald–Letnikov与 Caputo定义的曲线，并绘制出 d(t)函数的曲线，如图 3-13所示。
>> 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])
图 3-13 0.3阶导数
因为需要 C 0Dt 2.3y(t)函数，所以还需要的初始值是 y ′ (0)与 y ′′ (0)，这些值可以由原
函数通过符号运算的方式得出，结果可以变换回双精度的向量，就可以得出函数的 1.3阶与 2.3阶分数阶导数曲线。为了得到更好的近似效果，分数阶导数前几个点的值用样条插值方法重建，得出的曲线如图 3-14所示。
>> syms t; y(t)=sin(3*t+1); y00=sin(1);
y1=diff(y); y10=double(y1(0)); y2=diff(y1); y20=double(y2(0));
t=0:0.01:pi; y=@(t)sin(3*t+1); y1=caputo(y,t,1.3,[y00 y10],10);
y2=caputo(y,t,2.3,[y00,y10,y20],30);
axis left, plot(t,y1), axis right, plot(t,y2,'--')
C 1.3 0 t
D y(t)
C 2.3 0 t
D y(t)
图 3-14 1.3和 2.3阶导数双纵轴曲线
事实上，这里给出的求解函数有两个隐患：
（1）这里使用的插值结果引入了其他误差。
（2）实际应用中只有 y(t)的采样点数据，不可能有函数各阶导数的初值。所以，这里给出的函数不适合 Caputo导数计算，也不建议使用。第 4章将回到这个问题重新介绍，并给出高精度的数值解。
3.6分数阶微积分的性质与几何解释
和整数阶微积分类似，分数阶微积分也有很多自己的性质。本节首先介绍一些分数阶微积分的常用性质 ,然后，试图给出简单分数阶积分的物理解释与几何解释。不过，迄今为止，分数阶微积分尚没有被广泛接受的物理和几何解释。
3.6.1分数阶微积分的性质
·76·分数阶微积分学——数值算法与实现
证明这里先考虑证明 Grünwald–Letnikov分数阶微积分是线性的，记 τ = t . jh，可以通过下面的公式直接证明线性性质：
[(t.t0)/h]
.
GLDα af(t)+bg(t)= lim 1(.1)j α af(τ)+bg(τ)
t0 th→0 hα j
j=0 [(t.t0)/h] [(t.t0)/h]
..
1 α 1 α
= a lim (.1)j f(τ)+ b lim (.1)j g(τ)
h→0 hα j h→0 hα j
j=0 j=0 GLDα Dα
= af(t)+ b GL t g(t)
t0 tt0
对 Riemann–Liouville分数阶微积分定义而言，可以通过下面的推导证明它也是线性的：
t
RL 1 dk af(τ)+ bg(τ)
Dα af(t)+ bg(t)= dτ
t0 t (t . τ)1+α.k
Γ(k . α) dtk
t0 tt
a dk f(τ) b dk g(τ)
= dτ + dτ
(t.τ)1+α.k (t . τ)1+α.kΓ(k.α) dtk Γ(k . α) dtk
t0 t0 RLDα Dα
= af(t)+ b RL g(t)
t0 tt0 t
同样也可以证明 Caputo分数阶微积分是线性的。除了线性性质，分数阶微积分还有其他性质，这里将不加证明地列出一些常用的性质。关于这些性质的证明可以参见其他参考文献。
推论 3-2 (序贯导数)如果 n为整数，则下面的序贯分数阶微分满足：
RL RL
ddtnn t0 Dtαf(t)= t0 Dtn+αf(t)（3-6-5） n.1
dn . (t . t0)k.α.n
RLDα RLDn+α
t0 t dtn f(t)= t0 t f(t) . Γ(1 + k . n . α) f(k)(t0)（3-6-6）
k=0
证明由定理 3-10与推论 3-2可以直接证明本定理。这是一个非常重要的定理，将是后面将介绍的 Caputo微分方程框图求解方法的理论基础之一。
3.6.2分数阶积分的几何解释
整数阶微积分有很好的物理解释和几何解释。如果一个函数表示质点的位移，则其一阶导数表示该质点的瞬时速度；如果函数表示空间位置，则其一阶导数表示
·78·分数阶微积分学——数值算法与实现
切线方向。
相比之下，尽管分数阶微积分学也是一个有着 300多年悠久历史的学科，但分数阶微分与积分的定义至今尚没有被广泛接受的几何解释与物理解释。 Podlubny教授曾经给出了一个 Riemann–Liouville分数阶积分的合理几何与物理解释 [7，8]，这里将介绍该解释。
先回顾一下 Riemann–Liouville分数阶积分的定义
t
RL .γ
0 Dt f(t) = Γ(1 γ)(t . f(ττ)) 1.γ dτ（3-6-9）
0
可以引入一个辅助函数
g(τ)= Γ(γ 1+ 1) tγ . (t . τ)γ（3-6-10）
这样，分数阶积分可以直接写成 Stieltjes积分的形式：
t
RL .γ
0 Dt f(t)= f(τ)dg(τ)（3-6-11）
0
现在选择 3个坐标轴 τ，gt(τ)和 f(τ)建立起三维坐标系，则函数 f(t)可以分别在 τ-f(τ)坐标系、 g(τ)-f(τ)坐标系和整个三维坐标系下表示出来。这样，在 g(τ)-f(τ)坐标系下围成的面积可以认为是 f(t)函数分数阶积分，而 τ-f(τ)坐标系下的面积则为该函数的一阶积分。作者编写了一个 MATLAB函数 fence_shadow()绘制三维坐标系下的图形：
function fence_shadow(t,f,gam,key)
arguments
t(1,:), f, gam(1,1), key(1,1){mustBeMember(key,[1,2,3])};
end
x=t; z=f(x); tn=x(end);
switch key
case 1, y=(tn^gam-(tn-x).^gam)/gamma(1+gam);
axis([minmax(x),minmax(y),minmax(z)]), hold on
for i=1:length(x)-1, ii=[i,i,i+1,i+1]; x0=x(ii);
z0=[0 z(i:i+1) 0]; y0=zeros(1,4); patch(x0,y0,z0,'c')
y0=y(ii); patch(x0,y0,z0,'g')
x0=zeros(1,4); patch(x0,y0,z0,'r')
end, view(-37.5000,30)
case {2,3}, axis([minmax(x),minmax(z)]), hold on
for i=1:length(x)-1, x1=0:0.01:x(i+1);
tn=x1(end); y=(tn^gam-(tn-x1).^gam)/gamma(1+gam);
if key==2, plot(x1,y);
else, z=f(x1); plot(y,z,[y(end),y(end)],[z(end),0]); end, end, end, hold off end该函数的调用格式为 fence_shadow(t,f,γ,key)，其中， t为时间点向量， f为函数句柄，γ为积分阶次，key为标示，其取值对应文献 [7]中的 3种图形。这里的函数 g(τ)可以看作观测者自己的时间坐标系，该坐标系是物理时间 t扭曲的时间坐标系，如果在观测者自己的时间坐标系 g(t)下观测到物体的运动速度为 f(t)，则该观测者观测到的该物体位移就是 f(t)函数的 α阶积分。
例 3-13给定函数 f(t)= t +0.5 sin t，且 t ∈ (0, 10)，试绘制其 0.75阶 Riemann– Liouville积分的几何解释图。 解根据上面叙述，可以生成向量 t和原函数 f(t)，这样，可以绘制出 Riemann– Liouville分数阶积分的几何解释图，如图 3-15所示。可以看出，事实上图形是由填充的颜色块构造出来的。
图 3-15分数阶积分的几何解释示意图
>> t=0:0.5:10; f=@(t)t+sin(t)/2;
clf, fence_shadow(t,f,0.75,1), grid
set(gca,'xdir','reverse'), set(gca,'ydir','reverse')
可以看出，在左侧墙上的阴影部分面积是函数 f(t)的一阶积分，而右侧墙上的阴影部分则为所需的分数阶积分。如果选择 key为 2和 3，则可以得出如图 3-16和图 3-17所示的其他几何解释图。
>> clf, fence_shadow(t,f,0.75,2)
set(gca,'xdir','reverse','ydir','reverse'), ylim([0 2*pi])
figure, fence_shadow(t,f,0.75,3), xlim([0 2*pi])
如果图 3-15中的几何解释形象地将分数阶积分描述成栅栏上的影子，则在 τ增加时栅栏影子的长度也随之变化，图 3-16描述的是栅栏影子长度的变化过程，即在底面 τ-gt(τ)上的投影，而图 3-17给出的则是在墙面 gt(τ)-f(τ)上投影的变化过程。
·80·分数阶微积分学——数值算法与实现
f(τ)g t (τ)
时间 τ
图 3-16栅栏影子长度的变化示意图
gt(τ)
图 3-17栅栏在墙上影子变化的示意图
尽管 Podlubny教授的几何解释较好地描述了分数阶积分问题，但要想得出更易于理解的各种定义下分数阶微积分的几何与物理解释尚需更多努力。
本章习题
（1）设函数 u(t)= C，C为任意常数。根据定义 3-4和定义 3-5，用 MATLAB的符号
RLD0.5 RLD.0.5
运算功能计算 t u(t)和 0 t u(t)的表达式。
（2）试应用数值方证公式 法验0
Γ(v + 1)
GLDα (t . a)v =(t . a)v.α
at
Γ(v +1 . α)
其中，a为任意常数，α< 0时，v> .1；α> 0时，v> .α.。
（3）列举一个公式（3-3-2）的实例，并应用 MATLAB的符号运算证明该实例。例如，列举如下的实例： tτ sin τdτdt = t (t . τ ) sin τdτ
00 0
（4）试随机生成一组 λ、α的值，并用数值方法验证下面的等式：
RL eλt
0 Dα = t.αE1,1.α(λt)
t
（5）设函数 u(t)在定义域上连续，则其 Riemann–Liouville分数阶积分满足下面的关系式，试给定一个符合条件的函数 u(t)，用数值方法验证下面的等式：
RL RL .β RL .β RL RL .α.β
0 D.α 0 Du(t)= 0 D0 D.α u(t)= 0 Du(t)
tt tt t
其中，α> 0，β>0。
（6）根据 Riemann–Liouville分数阶导数的定义式，函数 u(t)的 Riemann–Liouville分数阶导数可以写成下面的形式：
dn
RL RL
0 Dα u(t)= 0 Dα.n u(t)
tt
dtn 其中，α> 0，n = .α.。试列举一个函数 u(t)，并用数值方法验证上面的等式。
（7）如果函数 f(t)具有任意阶的 Riemann–Liouville分数阶导数和积分，则有下面的等式[9]： RL RL 0 Dα.1
0 Dα [tf(t)] = t 0 Dαf(t)+ α RL f(t)
t tt
其中，α> 0。试列举一个满足上述条件的函数 f(t)，并用数值方法验证该等式。
（8）试利用 MATLAB的符号运算功能证明下面的等式[9]：
Γ(v + 1)
RL v
④ L 0 D.αt= , (α> 0,v > .1)
t
sα+v+1
RL eat 1
叫 L 0 Dt .α = sα(s . a) , (α> 0)
RL
③ L 0 D.α cos at =1 , (α> 0)
t
sα(s2 + a2)
（9）根据 Caputo分数阶导数的定义式，函数 u(t)的 Caputo分数阶导数可以写成下面的形式： C 0Dα u(t)= RL u
0 Dα.n (n)(t)
tt
其中，α> 0，n = .α.。试列举一个函数 u(t)，并用数值方法验证上面的等式。
（10）根据定理 3-6可以得到下面的等式：
C 0 D1.3e.0.5t =(.0.5)2 t0.7 E1,1.7 (.0.5t)
t
试用函数 caputo()和第 2章提供的函数 ml_func()检验上面的等式。
（11）由例 3-1可知常数的 Grünwald–Letnikov导数非零，其 Caputo定义的 0.5阶导数是什么？
（12）根据定理 3-9可知，如果函数 f(t)满足该定理的条件，则函数 f(t)的 Grünwald– Letnikov分数阶积分和导数与 Riemann–Liouville分数阶积分和导数等效。试列举满足条件的函数，并应用函数 glfdiff()和 rlfdiff()验证这一结论。
（13）根据定理 3-10可以得到下面的等式：
C RL
0 D1.5 sin t = 0 D1.5 sin t +1
tt
Γ(0.5) t0.5
应用数值算法检验上面的等式。
·82· 分数阶微积分学——数值算法与实现
（14）应用函数 rlfdiff()和 caputo()计算 RLDtαt2和 CDtαt2，分别取 α =0.5, 1, 1.5,
2, 2.5，比较相同阶数的Riemann–Liouvill0e分数阶导0数和 Caputo分数阶导数，得出二者之间的关系。
（15）根据定理 3-2可以得到下面的等式：
√
π
RL 0.5
0 D.0.5
t t=2 t
试用函数 rlfdiff()验证这一等式，选择不同的步长计算上面的等式，找出步长与计算误差的关系。
（16）证明下面的等式，并用函数 caputo()和 rlfdiff()验证该等式。
2
d2.k
. 2! dk
C 2e.t RL 2 e.t
0D1.50 D.0.5
t= tdt2.k
tt
k!(2 . k)! dtk
k=0
（17）证明下面的等式，并用函数 caputo()和 rlfdiff()验证这一等式。
π
C RL α..α.
0Dα sin t = 0 Dsin t + .α.
tt 2
其中，α> 0。
参考文献
[1] Oldham K B，Spanier J. The fractional calculus：Theory and applications of differ-entiation and integration to arbitary order[M]. San Diego：Academic Press，1974.
[2] Podlubny I. Fractional differential equations[M]. San Diego：Academic Press，1999.
[3] Podlubny I. Matrix approach to discrete fractional calculus[J]. Fractional Calculus and Applied Analysis，2000，3（4）：359–386.
[4] Dor.ák L. Numerical models for simulation the fractional-order control systems[R]. UEF-04-94. The Academy of Sciences Institute of Experimental Physics，Kosice, Slovak Republic，1994. [2022-10-2]. http://arxiv.org/pdf/math/0204108.pdf.
[5] Hilfer R. Applications of fractional calculus in physics[M]. Singapore：World Scientific， 2000.
[6] Caputo M. Linear models of dissipation whose Q is almost frequency independent II[J]. Geophysical Journal International，1967，13（5）：529–539.
[7] Podlubny I. Geometric and physical interpretations of fractional integration and differentiation[J]. Fractional Calculus and Applied Analysis，2001，5（4）：230–237.
[8] Bai L，Xue D Y，Meng L. Geometric interpretation for Riemann–Liouville fractional-order integral[C]. Proceedings of the Chinese Control and Decision Conference. Hefei， China，2020，3225–3230.
[9] 郭柏灵，蒲学科，黄凤辉.分数阶偏微分方程及其数值解 [M].北京：科学出版社，2011.