分数阶微分方程数值解的一种逼近方法

1、分数阶微分方程数值解的一种逼近方法By:PankajKumar,OmPrakashAgrawal摘要本文提出了一类分数阶微分方程（FDEs）的数值解方案.在这种方法中,FDEs被Caputo型分数阶导数所表现.Caputo型分数阶导数的属性可以让一个分数阶微分方程减弱为一个Volterra型积分方程.这样做了之后，许多研究Volterra型积分方程的数值方法也可以应用于寻找FDEs的数值解.本文总时间被划分为一组小区间,在两个连续区间中，用二次多项式逼近未知函数.这些近似被替换成转化的Volterra型积分方程由此获得一组方程.这些方程的解提供了FDE的解.这种方法被应用于解决两种类型的FDE

2、s，线性和非线性.用这里给出的方法得到的解能与解析解和其他方法的数值解较好的吻合.同时结果说明这种数值方法是稳定的.1.引言本文讨论分数阶微分方程的数值解.分数阶导数和分数阶积分近年来收到了广泛的关注.在许多实际应用中，分数阶导数和分数阶积分为考虑的系统提供了更加精确地模型.比如，分数阶导数已经被成功地运用到模拟许多粘性材料的依赖频率的阻尼行为.1980年之前，Bagley和Torvik提出了这个领域已经被研究的工作的一个回顾，并且说明了半阶导数模型可以非常好地描述阻尼材料的频率以来.另一些学者说明了分数阶导数和分数阶积分在电化学过程，电解质极化，有色噪声，粘性材料和混沌领域的应用.Maina

3、rdi,Rossikhin和Shitikova提出了分数阶导数和分数阶积分在一般固体力学，特定粘弹性阻尼模型中的应用的调查.Magin提出了分数阶微积分在生物工程的三个关键部分的回顾.分数阶导数和分数阶积分在其他领域的应用以及相关的数学工具和技巧还可以在许多其他文献上找到.系统模型中分数阶导数的引进大多会导致分数阶微分方程的出现.对某些特定的分数阶微分方程在通常系统条件下的解,已经有几种方法被找到.这些方法包括,拉普拉斯变换,傅里叶变换,模态综合法和特征向量展开法,数值法以及基于Laguerre积分公式的方法.然而，这些方法中大多数不能被应用到非线性分数阶微分方程.更进一步的，正如Diethe

4、lm等人指出的，这些方法很多只能应用到特定类型的分数阶微分方程，并且人们并不知道他们能否被推广.并且，在很多作者的研究成果中，并没有出现系统性的收敛性分析.最近，对于能被应用到线性和非线性分数阶微分方程的数值稳定数值逼近技巧，人们的兴趣愈发浓厚.这些方法技巧大多利用了分数阶微分方程可以被减弱为Volterra型积分方程的特性.因此，Volterra型积分方程的数值解法也可以应用到分数阶微分方程的解当中.Diethelm等人提出了分数阶微分方程数值解的一种PECE方法，其中P,C,E分别代表预测，校正和估计.这样一来很多学者又推广了应用于常微分方程和分数阶微分方程的Adams-Bashforth

5、-Moulton型预测-校正格式.这种方法的提出也是利用分数阶微分方程可以被转化为Volterra型积分方程的特点.这些作者同时提出了误差分析和用Richardson外推法改善数值精度的延伸.Ford和Simpson提出了一种阶数大于1的分数阶微分方程的数值解法.在该公式中，阶大于1的分数阶微分方程被减弱为阶小于1的分数阶微分方程，然后用相应的数值解法解由此导出的系统.在所有这些方法当中，节点之间的未知函数用线性函数逼近.KumarandAgrawal提出了阶数大于1的分数阶微分方程的数值解法.这种方法要求就y(t)和它的导数在时间节点上连续.本文基于古典分数阶微分方程可以转化为Volterr

6、a型积分方程的特点也提出了一种数值方法来逼近分数阶微分方程的解.特别地，我们用二次逼近函数来建立这种算法，结果说明这种方法可以被应用到寻求分数阶微分方程的数值解.我们还通过两个例子，线性和非线性问题的解决，说明了这种方法的高效和准确，并且这种数值方法是稳定的.2.数值算法关于分数阶导数的定义已经出现有好几种，它们包括Riemann-Liouville,Grun-wald-Letnikov,Weyl,Caputo,Marchaud,和Riesz分数阶导数.这里，我们规定使用Caputo导数.其中，Caputo导数的定义是A3、门八,门亿(n-1an),(1)其中，a是导数的阶数，n是比a大的最小

7、的整数.式(1)早在19世纪就在Liouville的论文中被提出，在Caputo的论文发表前一年它被Rabotnov所用.然而，在文献中，被(1)式所定义的分数阶导数作为Caputo导数被广泛认知.在接下来的讨论中，我们考虑含有Caputo导数的初值问题：D爲|(2)在初始条件：ik=0,1,.,n-1,(3)下的解，其中，f是任意函数，Iii是y的k阶导数，k二0,1,n-1是指定初值条件.假设这个函数关于参数和积分区间都是连续的，并且对于它的第二个参数满足Lipschitz条件.在纯数学中?Riemann-Liouville导数比Caputo导数应用更加广泛.然而，这里考虑Caputo导数

8、是因为以Riemann-Liouville导数为基础的分数阶微分方程要求y(t)在t=0点的导数和积分为0.一般来说，这些条件的物理意义不是已知的，并且在实际应用中，他们是不可用的.LorenzoandHartley讨论了寻找在更一般的情况满足下初始条件的正确格式的问题.在DiethelmandFord的文章中，方程(2)和(3)被证明可以等价描述为：只ti討汀D7/(ry.：r：.iT其中g(t)为：MH工；：川；(5)为了解释以二次多项式为基础的数值方法，我们假设我要求的是由(2)式定义的分数阶微分方程从0到T的积分.为了达到这个目的，我们把时间T等分成N份，令h=T/N，作为时间区间的每

9、一个部分.时间在网格点上被表示为=y敦;同时假设y(t)的数值逼近值被网格点；匚-：其中，L.-v(:i-：l是函数在第k个时间节点的值，山jN八，k=0,1,2是二次插值函数，其中下标(j,k)代表在第j+1,j=0,m步的第k个近似函数.我们首先确定y(t)在；和.处的值.把(7)式带入(4)式，并积分得到：aI工阿屮兀其中，和二打卅7哄汀(9)可以精确计算得到.注意式需要知道f在：和的值(或者间接地说,y的值).所决定.该方法的基本思想是在相邻的两个时间节点乜-I和1.-1上数值地获取函数y(t)的值，然后重复这个过程来接近所求积分直到取到终点T.为了便于接下来的讨论，我们规定如下记号：

10、yltjJ=y0h)心力，官(切=g(jh沁的,f(归)=f(jh$(jh)心耳这里的方法需要对方程(4)每一步求两次积分值这里有两种方法来达成目的.第一种，用一些近似函数逼近y(t),然后用一种数值方法确定式(4和-【作初始的估计.第二种，都用近似函数来显式地逼近y(t)和f(t,y(t)以及确定式(4)的积分.注意在这种情况下，+和一会作为参数出现在函数f当中.本文利用的是第二种逼近方法.现在，我们给出算法的详细思路.首先我们需要确定y(t),，：的值.用二次插值函数可以在区间0，：上逼近y(t)和f(t,y(t):沁说儿(6)以及( )为了得到；，在卩,上把f(t,y(t)近似为:f(坏

11、)=砂理叹(切其中h二;，八I：；：1!/-!是另一个二次插值函数.函数叭丿,k=l,2,3由下式给出,函数：1和a1是成立.(24),(25)的解是：5f、，(26)其中，、：!-!厲：、：I(27)是Mittag-Leffler函数的阶.Timet(sec.)图l.ay(t)的比较(O:解析解，X:本文数值方法)EmirinfuiicdoirL)furxOJS內口MdifYeren(vukvesuf/r!JA=0.4ft=0.2A=0.1=H.U5fr=H.025000000.S9.832E-4-3.OT6E-4一LME-4-4.730E-5-1.423E-5|出-110BE-4-2,11

12、3H-4-W1R-5-il3E-5-6-382E-62.4-2,Jfl3E-4一1.243E一4一39笳E-5-I.ITOF-5-3-537E一63.2-2,IOOE4-K.00JE一5-2.44J6F.-5-713JE-6-2-205E-&4-1.564E-4-5JIJE-5-L6R7E-5-5.021E-6-L487E一&-1.19IE-4-3.998E一5-L2IIE-5-3.592E-6-1.06JE-65.6-9.3S7E-5-3.CHSE-5-9.0E-6-2.685E-6-7.939E-76.4-7.41IE-5-2352E-5-7.W9E-6-2.07HE-6-6.I39C-7

13、表1.a时h在不同值下函数y(t)的误差对不同的a和h可以得到很多结果，这里给出其中一些.在各种情况下，我们另T=6.4s.考虑这个区间是因为它接近a=2的系统的时间.这里图(1)比较了a=0.75时的解析解和二次数值方法.在这种情况下，我们令h=0.1s.注意到这两个结果几乎完全重合.为了强调收敛性，表(1)列出了当a.025时的结果.注意到随着步长的缩小，误差也如期望一样的缩小了在大部分节点误差比R=E(2h)/E(h)(或E(h)=0(龙).图2.a不同时y(t)的比较(0：a=0.25,X：a=0.5,+：a=0.75,:a=0.95,*：a=1.)0246Timet(sec.)图3.

14、a(t)的比较(O:解析解，X:本文数值方法)Table!Errorinbincllohfcr21-5dinertnivaluesdfh1h=0.4k=02丹=CUh=0.05/r=0.0250(00fl4)0.8I.253E-31664E一42.6I5E-54.415E一67.67IE-?I.64.999E-4S.G&9E-5鳴E-52.2-OE一65.267E一72.4-I.346E-4fejOJSE一633E-68.M5E-7I.7I&E-7J.2-3.TO2E-4-3.438E-5-3.707E-6-4.766E-7-7.O61E-84-2.1SE-4一-3-S.523E-*一9.DM

15、E7-1.S54E-74.-I.O20E-4-2J0KI-5-4.0O9E：-6-7.327E-7-1.1171.-75.64.3O2L-5-4275E-6-1572E-7-6.&24E一86.49.386E-55.09讣-66.9CME-B-4.760F-B-J.266E-B表2.a时h在不同值下函数y(t)的误差图(2)展示了h=0.025,a(t)的数值结果.因为解析解和数值结果完全一致，因此图中没用画出解析解.当a=1时，精确解为y(t)=I注意到随着a越来0越接近1,数值解越收敛到解析解y（t）“，即在极限情况下，分数阶微分方程的解接近整数阶导数的解.更进一步地，我们给出了a1的一系

16、列结果.a1的结果是分离的，因为y（t）的斜率在ay（t）在a=1.5,h时的解析解和数值解.两个结果又一次几乎完全重合为了突出收敛性，表2给出了。（或E（h）=O（：）.这样一来，通过观察a1的收敛结果，可以知道误差的收敛阶为1+a（或E（h）=O（一十），即误差的阶不仅依赖于h，还依赖于导数的阶a.图4.a不同时y(t)的比较(0：a=1.25,X：a=1.5,+：a=1.75,:a=1.95,*：a=2.)i.e送.i.e-0.60.4*Increasing0.2.=0.4TimeI(sec.)图5.a（t）的比较（O:解析解，X:本文数值方法）Table3Compaiisanat亡no

17、rinreportedin35andobtainedusingthiswherni!a=0.5x=1.5hRcL35ThisschemeRef(35Thisscheme0.10.050.023-130E-03-WE-04L26E-O4-3,9KE-04-L.4M：-04-5.06E-05-5.46E-U4-I.28E-04-3.04E-052.45E-054.221：-067,40F-07表3.本文数值方法和文献(35)中y(t)的误差的比较.图4展示了h=0.025,a(t)的数值结果.又一次，因为解析解和数值结果完全一致，因此图中没用画出解析解.当a=2时，精确解为y(t)=cos(t).

18、注意到随着a越来越接近2,数值解逐渐收敛到整数阶导数的解.图2和图4展示的收敛结果十分重要，因为他们说明了在极限情况下，分数阶微分方程和他们的解逼近整数阶微分方程以及他们的解析解.表3比较了t=文献(35)的解的误差和用本文数值方法在a=0.5和1.25,h的解的误差.注意到本文的方法得到了更低阶的误差.这是因为，这里的方法是一种高阶方法.当a和h取其他值时这种趋势也能显现出来.3.2在第二个例子中，我们考虑一个如下定义的非线性方程：11)1険-严)-ywi3/2(28)且戸0严山(29)和之前一样，第二个初始条件仅适用于a1.(28)(29)的精确解在文献(35)中已给出，ydir3i:;f

19、(30)注意到当a1,方程的解在t=0处的斜率趋近于无穷.因此，他可能导致在t=0附近出现一个巨大的数值误差.Tble4Irrorin(uiuij的for厘栖1I孔嗣3diHr$g碎lueshseL.5/r=0.1ft=0.05h=0.025Zi=0A=0.05A=0.025D00DU000.1I.9&7E-4-6.I38E-6-5.5H4E-7I.447E-14.OI3E:-62.?IE-70.2-1.374E-4LSitiE5-1.057E-6I.05BE-4WIE-70.33O6tE-4-1552E-52.tWhE-49J88E-65.838E-7-15I9E-4-1.988E-5-I5

20、3IE-6L5WE-49J76E一65,837E一70.51.117E-4-1.621-5-1.2I3E-6JXOiE一57.3CK8E-63.96&E-7Lh一L220E-4-4.75KJ：-6一5I52E-6-6.LfifiE-X0.7-7.S29E-42171E-52.O93L一6-6.05BE一4-L033L-5-K.KyyE-70.H6.yE-46.K61E-55.953L一6-2.92BE一4-2.S79E一5-2J0IE-6-1M3E-31.425E-4l.l4C5i.04E-3-5.650E-5-4.L34E-6i.O2.867E-?2.4BBE-42.O4SE-5L.16IE

21、-39.566E-5-&.872.I-6表4.a=0.75和1.5,h取不同值下y(t)的误差.TableSCoinpariaoncferrorinreportedin35andobtainedusingthisschemeha=0.25.25Ref.351lussclienieRef.35Thisscheme0.I150E-012611：(3-5.53E-03-1.14R-040.05-I.50E-011471；-04-2.S0E-04-XUE-050.0254.09L-024.26E-051.63E-O5-3.12E-06表5.文献(35)中y(t)的误差和用本文数值方法得到的y(t)的误差的比较上面给出了一些在不同a和h下的数值结果.图5表示了h=0.05,a分别等于0.25，0.75，1.25，1.75时解析解和数值解的结果.由它可以得到(1).解析解和数值解基本重合，当a取其他值是，可以得到同样的结果.(2)尽管在t=0处斜率非常大，但是方法给出了非常精确的结果.(3)正如