白鹭^1,2, 薛定宇¹
1. 东北大学 信息科学与工程学院, 辽宁 沈阳 110819;
2. 沈阳大学 信息工程学院, 辽宁 沈阳 110044
收稿日期：2016-11-06
基金项目：国家自然科学基金资助项目(61174145，61673094)。
作者简介：白鹭(1982-)，男，辽宁沈阳人，东北大学博士研究生;
薛定宇(1963-)，男，辽宁沈阳人，东北大学教授，博士生导师。

摘要：设计了一种计算分数阶微积分的高精度数值算法, 提出了一种构造生成函数的简便方法.分析了基于快速Fourier变换的算法, 该算法误差较大的原因是应用了不准确的生成函数的系数, 而且没有考虑原函数的非零初值条件对计算精度的影响.新算法应用递推公式计算生成函数的系数, 并将原函数分解成零初值条件和非零初值条件两部分, 分别计算它们的分数阶微分和积分, 这样可以减小计算误差.误差分析和计算实例证明新算法具有很高的计算精度.
关键词：分数阶微积分生成函数高精度递推算法
High Precision Recursive Algorithm for Computing Fractional-Order Derivative and Integral
BAI Lu^1,2, XUE Ding-yu¹
1. School of Information Science & Engineering, Northeastern University, Shenyang 110819, China;
2. School of Information Engineering, Shenyang University, Shenyang 110044, China
Corresponding author: XUE Ding-yu, professor, E-mail: xuedingyu@mail.neu.edu.cn
Abstract: A high precision numerical algorithm was designed to compute fractional-order derivative and integral, and a simple method was proposed to construct the generating function. An algorithm based on fast Fourier transform was analyzed. It could be concluded that the reasons of its large computation error were using the inaccurate coefficient of the generating function and no considering the effect of nonzero initial condition of the original function on calculation precision. The recursive formula was used to compute the coefficient of the generating function in the new algorithm, what's more, the original function was decomposed into two parts, i.e., zero initial condition and nonzero initial condition, and their fractional-order derivative and integral were computed to decrease the computation error. The error analysis and the illustrative numerical examples showed that the computation accuracy of the new algorithm was very high.
Key Words: fractional-orderderivative and integralgenerating functionhigh-precisionrecursive algorithm
分数阶微积分可以简洁地描述具有历史记忆的物理过程, 被广泛应用在各工程领域中, 例如在控制理论中出现的一些新算法^[1], 为描述和控制复杂系统提供了方便.这些系统一般被描述成分数阶微分方程, 并且出现了许多求解这种方程的数值算法^[2-5].为了求解这些方程, 首先要解决的问题是如何计算分数阶微分和积分, 因此计算分数阶微积分的数值算法受到越来越多的关注, 常用的算法有分数阶线性多步法^[6]及基于快速Fourier变换(FFT)的算法^[7].当原函数的初值条件不为零时, 上述算法均不能得到高精度的结果.文献[8]提出了一种高精度的数值算法, 但此算法只能计算Caputo分数阶微分, 而且微分阶次只能在(0, 1)区间内.本文提出一种计算分数阶微分和积分的高精度数值算法, 该算法可以计算非零初值条件下任意阶次的分数阶微分和积分.
1 生成函数在分数阶微积分理论中, 有多种分数阶微积分定义, 文献[7]给出了这些定义的具体形式.当原函数的初值条件为零时, RL定义、GL定义和Caputo定义等价, 并可以应用式(1)近似计算^[10].

(1)

其中: 是分数阶微分或积分算子; α是分数阶微分或积分的阶次; 0和t分别是积分区间的下限和上限; y(t)是关于变量t的函数; h是算法的步长; k是离散点的下标; w是生成函数的Taylor展开式的系数.生成函数表达式为

(2)

其中:p是生成函数的阶数.如果原函数的初值条件为零, 式(1)的计算精度为O(h^p)^[9].下面提出一种构造生成函数的新方法.
定理1 ?当式(2)中的α=1时, 生成函数为

(3)

其中, 系数r_i满足

(4)

证明?由式(2)和式(3)得

(5)

令式(5)中的z=1, 得到式(4)中的第一个方程:

(6)

用z乘以式(5)的两边, 再求关于z的一次微分, 得

(7)

令式(7)中的z=1, 得到式(4)中的第二个方程:

(8)

用z乘以式(7)的两边, 再求关于z的一次微分, 得

(9)

令式(9)中的z=1, 得到式(4)中的第三个方程:

(10)

重复以上过程可得矩阵方程(4), 定理得证.
当计算一阶微分时, 可以直接应用定理1构造生成函数.当微分(积分)的阶次α≠1时, 首先应用定理1构造α=1的生成函数, 然后计算此生成函数的α次幂, 这样可以得到计算α阶微分或积分的生成函数.
2 基于FFT算法的局限性由于式(1)中的w是生成函数的Taylor系数, 所以可以计算生成函数在原点的Taylor展开式, 求出对应的Taylor系数, 然后应用式(1)计算分数阶微分或积分, 显然, 这种算法的效率很低.为了提高算法的速度, 文献[7]提出了一种基于FFT的算法, 下面介绍此算法的计算过程.
将生成函数的Taylor展开式写成

(11)

令式(11)中的z=e^－iφ, 则有

(12)

如果将两个函数的内积定义为

(13)

则e^inφ, e^－imφ是正交函数.将下面的内积记为w_k:

(14)

可知w是式(12)的Fourier系数, 因此可用FFT算法计算系数w, 这样可以提高算法的速度.以上的计算过程可以总结成下面的算法.
算法1 ?基于FFT的算法
1) 选择算法的步长, 计算原函数y(t)在离散点上的函数值;
2) 应用定理1构造生成函数, 应用FFT算法计算系数w;
3) 应用式(1)计算分数阶微分或积分.
由于FFT算法是近似计算, 得到的系数w并不是准确值, 因此会造成计算误差.另外, 算法1并没有考虑原函数的非零初值条件对计算精度的影响, 非零初值条件也会造成很大的计算误差.下面的例子将具体说明这一问题.
例1 ?应用基于FFT的算法计算e^-t的0.6阶RL微分, 下面的解析解可以检验数值解的精度.

(15)

其中, ε(·)为Mittag-Leffler函数.分别选择步长h=0.01和h=0.001, 应用定理1构造4阶生成函数:

(16)

应用FFT算法计算式(16)的Taylor系数w, 然后应用式(1)计算RL微分的数值解.将得到的数值解和解析解绘制成曲线, 如图 1所示, 可见数值解和解析解相差很大.正如之前的分析, 造成误差的原因为:应用FFT算法计算的w不准确; 原函数e^-t的非零初值条件也会造成很大的计算误差.
图 1(Fig. 1)

图 1 数值解与解析解Fig.1 Numerical solutions and analytical solutions

3 递推公式算法为了解决系数w不准确的问题, 下面推导计算w的递推公式.
定理2 ?生成函数:

(17)

其Taylor展开式为

(18)

当k=0时, 有w₀=r₀^α; 当k>0时, 计算w的递推公式为

(19)

当k－i < 0时, 有w_k－i=0.
证明?由式(17)和式(18)可得

(20)

令式(20)中的z=0, 可以得到w₀=r₀^α.
对式(20)的两边求关于z的一阶微分, 然后乘以, 得

(21)

将式(20)代入式(21), 得

(22)

在式(22)中, 当j－i+1 < 0时, 有w_j－i+1=0.
在式(22)的两边, z^j项的系数相等, 所以有

(23)

将式(23)中的j+1替换成k, 得

(24)

当k－i < 0时, 有w_k－i=0.整理式(24), 可得递推式(19), 定理得证.
为了消除非零初值条件产生的计算误差, 应用下面的方法将原函数转化成初值条件为零的函数.在算法中如果选择p阶生成函数, 则构造辅助函数:

(25)

为了使u(t)和y(t)在前p+1个离散点上的函数值相等, 式(25)中的系数c需要满足

(26)

其中, h为算法的步长.
将y(t)分解为u(t)+v(t), 由于u(t)和y(t)在前p+1个离散点上的函数值相等, 可以认为u(t)具有和y(t)相同的初值条件, v(t)具有零初值条件, 并且有

(27)

(28)

由于u(t)是幂函数的和, 可直接推导出u(t)的分数阶微分(积分)的解析表达式:

(29)

(30)

这样就可以得到下面的递推公式算法.
算法2 ?递推公式算法
1) 选择算法的步长, 计算原函数y(t)在离散点上的函数值;
2) 用定理1构造生成函数, 用定理2计算系数w;
3) 根据式(25)和式(26)构造u(t), 将y(t)分解为u(t)+v(t);
4) 应用式(1)计算v(t)的分数阶微分或积分;
5) 如果计算RL微积分, 应用式(29)计算u(t)的RL微积分; 如果计算Caputo微分, 应用式(30)计算u(t)的Caputo微分;
6) 应用式(27)计算y(t)的RL微积分, 或者应用式(28)计算y(t)的Caputo微分.
例2 ?应用本文提出的算法重新计算例1中的RL分数阶微分.为了比较两种算法的计算效果, 这里选择与例1相同的步长h=0.01，构造1~5阶的生成函数, 应用本文提出的算法计算数值解, 得到的计算误差如表 1所示.
表 1(Table 1)

表 1 当h=0.01时的计算误差Table 1 Computation errors when h= 0.01 t p=1 p=2 p=3 p=4 p=5 0.5 1.8E-3 1.2E-5 8.9E-8 7.1E-10 5.9E-12 1.0 1.7E-3 1.1E-5 8.6E-8 6.8E-10 5.7E-12 1.5 1.5E-3 1.0E-5 7.5E-8 6.0E-10 5.0E-12 2.0 1.3E-3 8.6E-6 6.4E-8 5.1E-10 4.3E-12 2.5 1.1E-3 7.4E-6 5.5E-8 4.4E-10 3.7E-12 3.0 9.6E-4 6.4E-6 4.8E-8 3.8E-10 3.2E-12 3.5 8.4E-4 5.6E-6 4.2E-8 3.4E-10 2.9E-12 4.0 7.5E-4 5.0E-6 3.7E-8 3.0E-10 2.8E-12 4.5 6.8E-4 4.5E-6 3.4E-8 2.7E-10 3.0E-12 5.0 6.2E-4 4.1E-6 3.1E-8 2.5E-10 3.4E-12

表 1 当h=0.01时的计算误差 Table 1 Computation errors when h= 0.01

由表 1可见, 本文算法的计算精度明显高于基于FFT算法的计算精度.当h=0.01时, 本文算法的计算误差可以控制在10^－12的级别.另外选择大步长h=0.1, 构造6~10阶的生成函数, 然后应用文本算法计算数值解, 得到的计算误差如表 2所示.可见, 即使选择非常大的步长, 本文的算法仍然可以计算出精度很高的数值解.
表 2(Table 2)

表 2 当h=0.1时的计算误差Table 2 Computation errors when h= 0.1 t p=6 p=7 p=8 p=9 p=10 0.5 3.1E-9 8.2E-11 3.0E-12 1.7E-13 4.4E-14 1.0 3.7E-8 2.9E-9 2.5E-10 2.0E-11 1.1E-12 1.5 3.7E-8 3.1E-9 2.5E-10 2.1E-11 1.5E-12 2.0 3.3E-8 2.8E-9 2.5E-10 1.7E-11 1.6E-12 2.5 2.8E-8 2.4E-9 2.5E-10 4.3E-11 3.9E-11 3.0 2.4E-8 2.1E-9 1.5E-10 2.5E-11 1.1E-10 3.5 2.1E-8 1.8E-9 4.5E-11 2.0E-10 6.9E-10 4.0 1.8E-8 1.5E-9 2.3E-10 5.9E-10 8.6E-9 4.5 1.6E-8 1.4E-9 4.9E-10 1.4E-9 2.5E-8 5.0 1.5E-8 1.3E-9 2.2E-10 1.0E-8 2.2E-7

表 2 当h=0.1时的计算误差 Table 2 Computation errors when h= 0.1

为了比较本文算法与文献[8]中的算法, 下面列举一个计算Caputo微分的实例.
例3 ?应用本文算法计算函数e^2t的α阶Caputo微分在t=1处的数值解, 数值解的精度可以用下面的解析解检验:

(31)

此例引用自文献[8], 这里选择与文献[8]相同的α和步长，构造6阶生成函数, 应用本文算法计算数值解, 计算误差如表 3所示.
表 3(Table 3)

表 3 两种算法的计算误差Table 3 Computation errors of two algorithms α h 本文算法 文献[8]算法 0.2 1/10 6.2124E-6 4.9025E-4 1/20 1.5455E-7 4.4094E-5 1/40 2.9006E-9 3.8638E-6 1/80 4.9372E-11 3.4232E-7 1/160 8.0647E-13 3.1440E-8 0.4 1/10 1.6234E-5 1.6156E-3 1/20 3.7989E-7 1.5681E-4 1/40 7.0048E-9 1.4478E-5 1/80 1.1841E-10 1.3103E-6 1/160 1.9256E-12 1.1851E-7 0.6 1/10 3.1264E-5 4.2309E-3 1/20 6.9099E-7 4.5820E-4 1/40 1.2550E-8 4.6839E-5 1/80 2.1093E-10 4.6479E-6 1/160 3.4515E-12 4.5469E-7 0.8 1/10 5.2502E-5 1.0190E-2 1/20 1.1024E-6 1.2525E-3 1/40 1.9780E-8 1.4521E-4 1/80 3.3094E-10 1.6333E-5 1/160 5.6612E-12 1.8089E-6

表 3 两种算法的计算误差 Table 3 Computation errors of two algorithms

如果选择相同的步长, 本文算法与文献[8]的算法相比, 计算误差减小了至少两个数量级.当步长减小时, 本文算法的误差减小得更为明显.比较结果说明, 与文献[8]的算法相比, 本文提出的算法将计算精度提高了两个数量级以上.
文献[8]中的算法不能计算阶次在(0, 1)区间以外的Caputo微分, 但本文算法可以计算这样的Caputo微分.另外, 本文算法还可以计算分数阶积分, 下面给出计算实例.
例4 ?计算函数e^-t的0.6阶RL积分在区间[0, 5]内的数值解, 可以应用下面的解析解检验数值解的精度:

(32)

选择步长h=0.01, 构造1~5阶的生成函数, 应用本文算法计算数值解, 计算误差如表 4所示, 可见计算结果也很精确.
表 4(Table 4)

表 4 计算误差Table 4 Computation errors t p=1 p=2 p=3 p=4 p=5 0.5 5.7E-4 3.6E-6 2.6E-8 2.0E-10 1.6E-12 1.0 1.5E-3 9.5E-6 7.0E-8 5.5E-10 4.5E-12 1.5 2.4E-3 1.6E-5 1.2E-7 9.2E-10 7.6E-12 2.0 3.3E-3 2.2E-5 1.6E-7 1.3E-9 1.1E-11 2.5 4.1E-3 2.7E-5 2.0E-7 1.6E-9 1.3E-11 3.0 4.9E-3 3.3E-5 2.4E-7 1.9E-9 1.6E-11 3.5 5.7E-3 3.7E-5 2.8E-7 2.2E-9 1.8E-11 4.0 6.4E-3 4.2E-5 3.1E-7 2.5E-9 2.1E-11 4.5 7.0E-3 4.6E-5 3.4E-7 2.7E-9 2.3E-11 5.0 7.6E-3 5.0E-5 3.7E-7 3.0E-9 2.5E-11

表 4 计算误差 Table 4 Computation errors

4 误差分析对本文提出的算法进行误差分析, 首先介绍下面的引理, 此引理引自文献[6].
引理1??用式(1)计算函数y(t)=t^q－1(q>0)的分数阶微积分, 则有

(33)

由引理1可知, 式(1)的计算精度不仅与生成函数的阶次有关, 而且与原函数的初值条件有关.这里假设原函数y(t)在原点处解析, 可以写出y(t)在原点的Taylor展开式为

(34)

本文算法将y(t)分解为u(t)+v(t), 由u(t)和v(t)的表达式可得

(35)

应用式(29)和式(30)可计算u(t)的分数阶微积分的精确值.当应用式(1)计算v(t)的分数阶微积分时, v(t)具有零初值条件, 即引理1中的p=q, 因此计算的精度只与生成函数的阶数有关, 而与原函数的初值条件无关.所以递推公式算法的计算精度为O(h^p).
5 结论1) 基于FFT算法误差较大的原因是应用了不准确的生成函数系数, 而且在计算中未考虑原函数的非零初值条件对计算精度的影响.
2) 证明了计算生成函数系数的递推公式, 并设计了递推公式算法, 此算法可以计算非零初值条件下任意阶次的分数阶微分和积分.
3) 在相同的步长下, 本文算法的计算精度高于文献[8]的算法两个数量级以上.即使选择大步长h=0.1, 本文算法的计算误差仍然可以控制在10^－10数量级.
参考文献

[1]	Wei Y H, Tse P W, Du B, et al. An innovative fixed-pole numerical approximation for fractional order systems[J].ISA Transactions, 2016, 62: 94–102.DOI:10.1016/j.isatra.2016.01.010
[2]	Xue D Y, Bai L. Numerical algorithms for Caputo fractional-order differential equations[J].International Journal of Control, 2017, 90(6): 1201–1211.DOI:10.1080/00207179.2016.1158419
[3]	Xue D Y, Bai L. Benchmark problems for Caputo fractional-order ordinary differential equations[J].Fractional Calculus and Applied Analysis, 2017, 20(5): 1305–1312.
[4]	Li M M, Wang J R. Finite time stability of fractional delay differential equations[J].Applied Mathematics Letters, 2017, 64: 170–176.DOI:10.1016/j.aml.2016.09.004
[5]	Li C P, Zeng F H. Numerical methods for fractional calculus[M]. Boca Raton: Chapman & Hall Crc, 2015: 50-100.
[6]	Lubich C. Discretized fractional calculus[J].SIAM Journal on Mathematical Analysis, 1986, 17(3): 704–719.DOI:10.1137/0517050
[7]	Podlubny I. Fractional differential equations[M]. New York: Academic Press, 1999: 41-91.
[8]	Cao J X, Li C P, Chen Y Q. High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (Ⅱ)[J].Fractional Calculus and Applied Analysis, 2015, 18(3): 735–761.