第1页 / 共6页
第2页 / 共6页
第3页 / 共6页
第4页 / 共6页
第5页 / 共6页
第6页 / 共6页
Caputo分数阶导数的一种二阶逼近方法.pdf
Caputo fractional derivatives
Huang Fenghui
GuangZhou 510641)
(Department of Mathematics, School of Sciences, South China University of Technology,
Abstract: Fractional derivatives and fractional integrals provide more accurate models of the systems
than ordinary derivatives and integrals do. However, effective numerical method for fractional
derivatives and integrals is still very few, and most of them offer worse than second-order accuracy. In
this paper we propose a second-order accurate algorithm for Caputo fractional derivative base on the
algorithm for numerical fractional Riemann-Liouville integration. The method can be considered the
generalization of the three points formula for the first-order derivative. Two numerical examples are
given，and illustrative the effectiveness of the method.
Keywords:Numerical solution; Caputo fractional derivative; fractional integral
中国科技论文在线
http://www.paper.edu.cn
Caputo 分数阶导数的一种二阶逼近方法#
黄凤辉*
（华南理工大学理学院数学系，广州 510641）
摘要：分数阶微积分描述的模型往往比经典的整数阶微积分模型能更精确反映复杂的实际问
题。然而，当前行之有效的分数阶微积分的数值算法并不多，而且能达到二阶精度的就更少
了。 本文根据 Riemann-Liouville 分数阶积分的一种二阶精度的逼近方法，导出了求解
Caputo 分数阶导数的一种二阶方法。该方法可以看作是一阶导数二阶精度的三点公式的推
广。最后给出两个数值例子，验证了方法的有效性及其收敛精度。
关键词：数值解；Caputo 分数阶导数；分数阶积分
中图分类号：O175.2
A second-order accurate numerical approximation for
0 引言
分数阶微积分是传统的整数阶微积分概念的推广，是专门研究任意阶积分和微分的数学
性质及其应用领域。最早提出这一思想的数学家是 Leibniz（1695 年），在他之后也有许多
数学家先后在该领域做出了杰出的贡献。近三百年来，分数阶微积分这一重要的纯数学分支
已渐成体系，但 是 发 展 极 为 缓 慢 ， 而 且 分 数 阶 微 积 分 理 论 还 不 完 善 。 其 主 要 原 因
是 分 数 阶 微 积 分 的 多 种 定 义 之 间 的 不 等 价 性 及 其 非 局 部 性 ， 使 得 其 缺 乏 明 显 的
几 何 意 义 ， 应 用 一 时 受 到 限 制 。 直到 Mandelbrot 提出分形学说，将 Riemann-Liouville
分数阶微积分用以分析和研究分形媒介中的布朗运动，分数阶微积分才在许多学科的现代工
程计算中得以广泛关注和应用。另外，Podlubny 和 Machado 分别给出了分数阶导数的几何
及应用解释[1,2] 。 由 于 应 用 的 推 广 ， 分 数 阶 微 积 分 在 近 几 年 得 到 了 比 较 大 的 一 个
发 展 ，已 被 广泛应用于物理、化学、地下水模拟、系统生物、金融、水文、生态环境等诸
多领域[3-6]，引起了越来越多的专家和学者的广泛兴趣和关注。然而，分数阶微积分在数 值
算 法 方 面 还 是 比 较 有 限 。 现 有 的 算 法 包 括 基 于 Grunwald-Letnikov 分数阶导数的逼近
公式（称之为 G 算法，1 阶精度）[7]；Lubich 于 1986 年提出的所谓 R 算法（计算分数阶积
分，1-2 阶精度）和 L 算法（计算α阶分数阶导数，0
< 时为 1 阶精度）[8-9]；以及 Lubich
分数阶线性多步法[10]；还有 Diethelm 于 1997 年提出的所谓 D 算法（计算α阶分数阶导数，
基金项目：高等学校博士学科点专项科研基金（新教师基金，编号：20070561040）；华南理工大学中央高
校基本科研业务费专项资金(编号：2009ZM0050)
作者简介：黄凤辉，（1974-），女，副教授，主要研究方向：微分方程数值解. E-mail: huangfh@scut.edu.cn
1α<
- 1 -
中国科技论文在线
http://www.paper.edu.cn
1α<
0
< 时为 2 α− 阶精度）[11]。这些方法最终都是表示成一些给定的离散点上的函数值的
加权和，可以看作是整数阶导数的差商逼近的推广，统称为分数阶导数的差商逼近。本文正
是将一种求解 Riemann-Liouville 分数阶积分的二阶精度的 R 算法推广应用到求解 Caputo 分
数阶导数上，得到一种新的二阶精度方法。
1 分数阶微积分介绍
1.1 定义
一直以来，人们不断探讨分数阶微积分的定义及性质研究。在不同科学技术领域的应用
问题中，分数阶微积分采用的定义也有所不同。主要有Riemann-Liouville, Weyl, Reize, Caputo,
compos 等定义。下面给出本文要用到的分数阶微积分的定义。
f x x > ，若存在实数
p
Rµ∈ ，对于实函数 ( ),
∞ ，使
∈
0
)
那么称函数 ( )
f x 属于空间Cµ ；且 ( )
,
µ>
f x C
( )
(0,
1
f x 属于空间 ,mC m N
µ
∈ ，当且仅
1
得
定义 1：设
x f x
f x
( ),
( )
p
=
x ∈ Cµ 时。
mf
)( )
当 (
f x ∈
定义 2：若 ( )
1
∫
(
)
Γ
α
aJ f x
f x
( )
( ).
0
a = 时，简记 0
特别地，当 0
J f x
( )
α
a
)
1
α
−
=
−
=
x
y
(
a
x
Cµ µ≥ − ，则
1
,
0α≥ 阶 Riemann-Liouville 分数阶积分定义为：
f y dy
( )
,
α
>
0,
x
>
a
,
J f x
( )
α
α=
J f x
( )
。Riemann-Liouville 分数阶积分算子具有
如下一些性质，
f x ∈
引理 1：若 ( )
J J f x
( )
α β
a
a
（1）
=
（3）
aJ
(
α
x a
−
)
γ
=
J
（5）
α
cos(
ax
)
=
(
1),
α β+
a
mCµ µ
J
≥ −
f x
( )
1)
(
Γ +
γ
(
1)
+ +
Γ
α γ
ax
( 1) (
k
−
∞
∑
k
(
2
Γ
+
α
x
α
(
=
0
k
αβ
,
>
0, 0
γ
1,
≠ ≥ − 则有
J J f x
( )
α β
a
a
=
， （2）
x a
−
+
)
α γ
， (4)
ax
J e
α
=
x
α
∞
∑
k
=
0
)
+
k
2
1)
, (6)
α
J
sin(
ax
)
=
，
J J f x
( )
β α
a
a
ax
(
)
k
k
(
1)
α
+ +
ax
( 1) (
k
−
∞
∑
k
(
2
Γ +
α
Γ
=
0
k
x
α
,
2
)
+
k
1
+
2)
.
有关 Riemann-Liouville 分数阶积分算子的性质请参见相关文献[7] 。
∈ 定义α阶 Caputo 分数阶导数为：
定义 3：函数
f x C m N
( )
,
,
m
−∈
1
D f x
( )
α
*
a
=
J
)
(
x
a
)
m
f
[
(
y
x
α
−
−
m
a
∫
)
α
x
( )]
1
⎧
⎪Γ
m
(
−
⎪
= ⎨
d
m
⎪
⎪⎩
dx
a D f x D f x
( )
a = 时，简记 *
α
,
( )
1
f x ∈
α=
*
f x
( ),
m
特别地，当 0
定义 4：若函数 ( )
m
1
α
− −
f
(
m
)
y dy m
( )
,
1
− <
<
α
m
,
(1)
α
=
m
.
.
Cµ µ≥ − ，定义其α阶 Riemann-Liouville 分数阶导数为：
α
(
D f x
)
=
a
d
m
dx
m
(
( )
J f x
α
a
)
=
d
m
dx
m
⎡
⎢
Γ⎣
1
(
)
α
x
∫
a
(
x
−
y
1
α
−
)
(
(
f y dy m
)
⎤
⎥
⎦
1
− <
≤
α
m
)
.
- 2 -
中国科技论文在线
http://www.paper.edu.cn
引理 2: 若 1m
mα− < ≤ , 函数
f x C −∈
( )
m
1
，并设
( )
a D f xα
存在，那么
D f
α
a
( )
x
=
D f
α
*
a
( )
x
+
1.2 分数阶积分的 R2 算法
m
1
−
∑
k
=
0
x
Γ
k
α
−
(
k
f
−
k
( ) (
)
a
)
1
+
α
. (2)
如何数值计算分数阶积分、微分方程，特别是线性和非线性的 Abel 积分方程，以及含
分数阶积分算子的积分微分方程，其关键是寻找好的数值方法求解 Riemann-Liouville 积分
算子。Lubich 提供了一类方法（R 算法），其中一种具有二阶精度，即所谓的 R2 算法，由
如下定理给出。另外，为了书写简洁，我们取 a=0, 但不是必要的。
kh h T N
,
定理 1：（R2 算法[8-9]） 设
f x C
( )
[0,
则
∈
=
=
/
,
2
T x
k
],
(
J f x
α
( ))(
x
k
)
=
其中
k
∑
j
=
0
a
j k
,
f x
(
j
)
+
E f h
,
,
)
α
J
(
，
− −
j
1)
1
α
+
,
j
=
1
≤ ≤ −
j
=
0
j
k
k
1
, cα是与α有关的常数。
⋅
(
(
j k
,
k
a
=
+
Γ
2)
h
α
(
α
+
1)
1
α
+
j
− +
k
k
(
,
1)
α
α
− −
−
k
j
2(
1)
)
1
1
α
+
α
+
−
−
k
(
⎧ −
⎪
k
(
⎨
⎪
1,
⎩
E f h
c
,
,
αα
J
注：当 1α= 时，R2 算法退化成复化梯形公式：
b
∫
所以 R2 算法也可以看作是复化梯形公式的推广。
h O h
(
)
2
α
f b
( ) 2
f x dx
( )
f a
( )
O h
(
∑
f x
(
h
2
x
k
x
k
)]
≤
=
),
+
=
+
+
)
(
)
f
1
−
k
1
=
[
M
∞
2
a
2
k
''
= +
a kh h
,
=
(
b a M k
,
−
) /
=
0,1,
M
,
2 Caputo 分数阶导数的数值逼近
2.1 数值逼近公式
由 Caputo 分数阶导数的定义 3，并应用 R2 算法可以得到如下求解 Caputo 分数阶导数
的方法（在定理 1 中，α用 m α− 替代， f 用 (
mα
,
定理 2：设
f x C
( )
)mf 取代）。
=
T m
1
− <
并令
[0,
],
≤
2
kx
kh h T N
=
/
,
,
则
(
D f x
α
*
( ))(
x
k
)
=
Γ
(
{
[(
k
m
1
α
− +
−
1)
−
(
k m
− + −
α
1)
k
m
]
α
−
f
(
m
)
(
x
0
)
+
f
(
m
)
(
x
k
)
（3）
2)
− +
j
1)
m
1
α
− +
−
2(
k
−
j
)
m
1
α
− +
+
(
k
− −
j
1)
m
1
α
− +
(
m
)
]
f
(
x
j
}
)
+
E
D
(
f h
,
)
,
α
m
−
α
m
+∈
h
m
α
− +
k
1
−
∑
[(
+
k
j
1
=
其中
E f h
,
)
,
α
D
(
≤
c
m
−
α
f
(
m
+
2)
(
x
k
∞
m
)
h O h
(
2
−
α
=
2
)
, mc α− 是与α有关的常数。
在具体进行数值计算时，（3）中整数阶导数可以用一些高阶方法（大于等于 2 阶精度）
≤ 为例讨论，此时 1m = ，（3）中的各个节点
1α<
离散，最后仍保证二阶精度。下面以0
上的一阶导数分别采用如下二阶方法离散：
f x
(
'
0
)
=
f x
3 (
−
0
f x
) 4 (
+
1
h
2
)
−
f x
(
2
)
+
O h
(
2
)
， （4）
- 3 -
中国科技论文在线
http://www.paper.edu.cn
f x
(
k
−
2
f x
(
'
k
)
=
f x
(
j
1
+
f x
(
'
)
=
j
f x
) 4 (
−
k
h
2
f x
(
j
1
−
)
−
h
2
f x
) 3 (
+
k
)
1
−
+
O h
(
2
),
k
≥ （5）
2
)
+
O h
(
2
), 1
≤ ≤ − ， （6）
1
k
j
特别地，当 1α= 时，公式退化成一阶导数的三点逼近公式(5)。
另外，对于 Riemann-Liouville 分数阶导数，根据它与 Caputo 分数阶导数的关系式(2), 可
以推出其对应的数值逼近格式，这里就不多讨论。
2.2 精度及收敛性测试
为了验证方法的有效性及其收敛性，我们选择了存在解析导数的函数来进行数值模拟。
≤ 的情况，且对(3)中的一阶导数采用二阶逼近格式(4)-(6)离散。
这里，我们只讨论 0
另外，下面所有的数值结果都是计算 1x = 处的导数值。
1α<
数值算例 1 考虑函数
f x
( )
E
µ=
(
−
xµ
)
,µ为给定的常数。由分数阶微积分的性质（引理
1）及 Caputo 分数阶导数的定义，可以推出
D E
α
µ
(
µ
−
x
)
=
x E
α
−
µ α
,1
−
(
−
x
µ
)
−
其中
,E
zαβ 是 Mittag-Leffler 函数，它由如下级数形式定义：
( )
α
−
x
(1
Γ −
)
α
z
k
k
∞
∑
E
αβ
( )
z
=
,
=
(
)
α β
Γ
+
k
0
特别地， ( )
( )
E z
,1E
z
表 1 给出了
α=
α
,
α β
>
0,
>
0
.
.
0.5α= 时不同步长 h 和µ下的绝对误差值。由表格可以看出，对于不同的
µ，绝对误差随着步长 h 的减小而减小，即方法收敛。为了进一步观察方法的精度，我们描
出了误差与步长 h 的关系，如图 1 所示，其中两坐标均为对数坐标， 2.5µ= 。虚线分别表
x= 的直线。由图可以看出，不
示不同α取值下步长与绝对误差的关系，实线则是表示 2
x= 的直线平行，表明方
同α取值下，绝对误差的对数与步长对称均为直线关系，且与 2
法是二阶收敛精度，即
y
y
2
,
α =
)
O h
(
)
。
DE f h
,
(
h
0.1
0.05
0.01
0.005
0.001
0.0005
DE f h
,
(
(
6.3)µ=
,0.5)
0.00001935610003
0.00001195487295
0.00000096195163
0.00000027206927
0.00000001260561
0.00000000325442
表 1 不同步长 h 下的绝对误差
Tab. 1 Absolute errors varying with step h
DE f h
,
,0.5)
DE f h
,
(
,0.5)
(
3.5)µ=
(
(
2)µ=
0.00064730487006
0.00023492814824
0.00001342483722
0.000003598329595
0.000000156905360
0.000000039996080
0.00108202906989
0.00037182040156
0.00002022289064
0.00000537095047
0.000000231643458
0.000000058905831
- 4 -
中国科技论文在线
http://www.paper.edu.cn
)
r
o
r
r
e
(
g
o
l
0
1
α=0.1
α=0.5
α=0.8
α=1
y=2x
-4
-6
-8
-10
-12
-14
-16
-18
µ=2.5
-20
-8
-7
-6
-5
log10(h)
-4
-3
-2
图 1 误差与步长的关系
Fig. 1 Errors as a function of the time step
)
f x
数值算例 2 考虑函数 ( )
=
e
ax
+sin(
bx
, 其中 ,a b 为给定的常数。由分数阶微积分的
性质（引理 1）及 Caputo 分数阶导数的定义，可以推出其解析导函数为
x
−
α
∞
∑
k
1
=
(
k
k
ax
)
α
− +
+
)
1
Γ
(
bx
1
−
α
∞
∑
k
=
0
2
( 1) (
k
−
(
2
Γ
bx
)
α
− +
k
.
k
2
)
表 2 给出了不同步长 h 和α下的绝对误差值。由表格可以看出，对于不同的α，绝对误
差随着步长 h 的减小而减小，即方法收敛。图 2 描出了误差与步长 h 的关系，其中两坐标均
为对数坐标， 0.5α= 。虚线分别表示不同 ,a b 取值下步长与绝对误差的关系，实线则是表示
f x ，绝对误差对数值与步长对数值
y
均为直线关系，且与 2
。
)
x= 的直线。由图可以看出，对于不同的具体函数 ( )
x= 的直线平行，表明方法是二阶收敛精度，即
,
α =
O h
(
2
y
)
2
DE f h
,
(
以上两个算例均表明了方法的有效性及其二阶收敛精度。
h
0.10
0.05
0.025
0.0125
0.00625
0.001
0.0005
表 2 不同步长 h 下各分数阶导数的绝对误差 (
DE f h α
)
,
,
(
(
(
Tab. 2 Absolute errors varying with step h for various
DE f h
,
,0.1)
0.00821495383730
0.00235706967277
0.00062925200200
0.00016249195729
0.00004128751458
0.00000107171059
0.00000026829336
DE f h
,
,0.5)
0.00595904064952
0.00180792962313
0.00049937230388
0.00013210543216
0.00003417121998
0.00000091209006
0.00000022966936
DE f h
DE f h
,
,
0.00620936876143
0.00286500073621
0.00144931994993
0.00102775098895
0.00035004321333
0.00030900489624
0.00008600952365
0.00008649712280
0.00000787679963
0.00002131681139
0.000000300586901 0.00000054174604
0.000000082659080 0.00000013534202
,0.7)
,1)
(
- 5 -
中国科技论文在线
http://www.paper.edu.cn
fx)=e-2x
f(x)=sin x
f(x)=e-2x+sin x
f(x)=ex+sin(2x)
y=2x
-2
-4
-6
-8
-10
-12
-14
-16
)
r
o
r
r
e
(
0
1
g
o
l
α=0.5
-18
-8
-7
-6
-5
log10(h)
-4
-3
-2
图 2 误差与步长的关系
Fig. 2 Errors as a function of the time step
3 结论
1α<
本文给出了一种数值求解 Caputo 分数阶导数的二阶精度方法。数值测试结果表明了
≤ 时方法的有效性。特别地，当 1α= 时，方法退化为一阶导数的三点差商逼近格式
0
(5)。另外，该方法也可以推广应用到数值计算 Riemann-Liouville 分数阶导数。
[参考文献] (References)
[1] Podlubny I. Geometric and physical interpretation of fractional integration and fractional differentiation[J].
Fractional Calculus and Applied Analysis, 2002, 5(4): 367～386 .
[2] Tenreiro Machado J A. Fractional Derivatives: Probability Interpretation and Frequency Response of Rational
Approximations [J]. Communications in Nonlinear Science and Numerical Simulations, 2009, 14(9–10): 3492～
3497.
[3] Diethelm K, Ford A D. On the solution of nonlinear fractional-order differential equations used in the modeling
of viscoplasticity[A]. F. Keil, W. Mackens, H. Vo , J. Werther (Eds.). Scientisc Computing in Chemical
Engineering II. Computational Fluid Dynamics, Reaction Engineering, and Molecular Properties[C], Heidelberg,
Springer-Verlag, 1999. 217～224.
[4] Benson D A, Wheatcraft S W and Meerschaert M M. Application of a fractional advection-dispersion
equation[J], Water Resources Research, 2000, 36(6): 1403～1412.
[5] Mainardi F, Raberto M, Gorenflo R, Scalas E. Fractional calculus and continuous-time finance II: the
waiting-time distribution[J], Physica A, 2000, 287: 468～481.
[6] Su N, Liu F and Anh V. Simulating seawater intrusion in aquifers using modified Fokker-Planck equation and
Boussinesq equation subject to phase-modulated tidal waves[A]. In: C. Gulati, Y.-X. Lin, S. Mishra & J. Rayner
(eds.): Advances in Advances in Statistics, Combinatorics & Related Areas[C]. Singapore, World Scientific, 2002,
320～331.
[7] Podlubny I. Fractional differential equations[M]. San Diago, Academic Press, 1999.
[8] Lubich C H. Discretized fractional calculus[J]. SIAM J. Math. Anal.,1986, 17:704～719.
[9]Diethelm K, Ford N and Freed A D. Detailed error analysis for a fractional Adams method[J]. Numerical
Algorithms, Kluwer Academic Publishers，2004, 36(l):31～52.
[10] Lubich C. Fractional linear multistep methods for Abel-Volterra integral equations of the second kind[J].
Math. Comp., 1985, 45:462～469.
[11] Diethelm K. An algorithm for the numerical solution of differential equations of fractional order[J].
Electron. Trans. Numer. Anal., 1997, 5:1～6.
- 6 -
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
