2018, Vol. 39
2. 沈阳大学 信息工程学院, 辽宁 沈阳 110044
2. School of Information Engineering, Shenyang University, Shenyang 110044, China
分数阶微积分可以简洁地描述具有历史记忆的物理过程, 被广泛应用在各工程领域中, 例如在控制理论中出现的一些新算法[1], 为描述和控制复杂系统提供了方便.这些系统一般被描述成分数阶微分方程, 并且出现了许多求解这种方程的数值算法[2-5].为了求解这些方程, 首先要解决的问题是如何计算分数阶微分和积分, 因此计算分数阶微积分的数值算法受到越来越多的关注, 常用的算法有分数阶线性多步法[6]及基于快速Fourier变换(FFT)的算法[7].当原函数的初值条件不为零时, 上述算法均不能得到高精度的结果.文献[8]提出了一种高精度的数值算法, 但此算法只能计算Caputo分数阶微分, 而且微分阶次只能在(0, 1)区间内.本文提出一种计算分数阶微分和积分的高精度数值算法, 该算法可以计算非零初值条件下任意阶次的分数阶微分和积分.
1 生成函数在分数阶微积分理论中, 有多种分数阶微积分定义, 文献[7]给出了这些定义的具体形式.当原函数的初值条件为零时, RL定义、GL定义和Caputo定义等价, 并可以应用式(1)近似计算[10].
|
(1) |
其中: 
|
(2) |
其中:p是生成函数的阶数.如果原函数的初值条件为零, 式(1)的计算精度为O(hp)[9].下面提出一种构造生成函数的新方法.
定理1 当式(2)中的α=1时, 生成函数为
|
(3) |
其中, 系数ri满足
|
(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) |
则einφ, e－imφ是正交函数.将下面的内积记为wk:
|
(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 Numerical solutions and analytical solutions |
为了解决系数w不准确的问题, 下面推导计算w的递推公式.
定理2 生成函数:
|
(17) |
其Taylor展开式为
|
(18) |
当k=0时, 有w0=r0α; 当k>0时, 计算w的递推公式为
|
(19) |
当k－i < 0时, 有wk－i=0.
证明 由式(17)和式(18)可得
|
(20) |
令式(20)中的z=0, 可以得到w0=r0α.
对式(20)的两边求关于z的一阶微分, 然后乘以
|
(21) |
将式(20)代入式(21), 得
|
(22) |
在式(22)中, 当j－i+1 < 0时, 有wj－i+1=0.
在式(22)的两边, zj项的系数相等, 所以有
|
(23) |
将式(23)中的j+1替换成k, 得
|
(24) |
当k－i < 0时, 有wk－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 当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 当h=0.1时的计算误差 Table 2 Computation errors when h= 0.1 |
为了比较本文算法与文献[8]中的算法, 下面列举一个计算Caputo微分的实例.
例3 应用本文算法计算函数e2t的α阶Caputo微分在t=1处的数值解, 数值解的精度可以用下面的解析解检验:
|
(31) |
此例引用自文献[8], 这里选择与文献[8]相同的α和步长，构造6阶生成函数, 应用本文算法计算数值解, 计算误差如表 3所示.
|
|
表 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 Computation errors |
对本文提出的算法进行误差分析, 首先介绍下面的引理, 此引理引自文献[6].
引理1 用式(1)计算函数y(t)=tq－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(hp).
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.
|