分数阶微分方程初值问题的Runge-Kutta型法
2018-08-06樊孝仁
王 颖,樊孝仁
(太原工业学院 理学系,山西 太原 030008)
0 引言
分数阶导数在描述不同物质的记忆与遗传性质方面提供了有力的工具,在科学和工程的不同领域都用分数阶方程(组)来描述动力系统.分数阶微积分在数学、物理、化学、控制等领域都有广泛的应用,已经吸引了许多科学家和工程学家的关注和研究[1].同整数阶微分方程一样,许多分数阶微分方程的解都是理论上存在,其解析解并不能明确给出.然而,除了模型,对于方程出现突然发散,收敛或分歧的情况下,对于抓住其关键点是至关重要的.因此,高精度数值解仍然是需要的[2].因此,许多分数阶常微分方程的数值解法被提出,最常见的技巧是区域分解法[3-5],变分法[6-7],分数阶微分变换法[8-12],分数阶差分法[13]及幂级数解法[14].同时,整数阶微分方程的经典方法也经常用到.如 :Laplace变换法[15].然而,常微分方程初值问题的数值解法的经典方法之一Runge-Kutta方法却很少被人们所提到.
本文将把Runge-Kutta方法应用到分数阶常微分方程的初值问题中,给出其算法格式.
1 定义
以下是该文可能用到的基本定义.
定义1[16]对t>0,函数x(t)的阶数为α∈R+的分数阶积分定义为
其中,Γ(·)是伽马函数.
定义2[16]对t>0,函数x(t)的阶数为α∈[n-1,n),n∈Z+的Riemann-Liouville分数阶导数定义为
定义3[16]如果f(t,x(t))是连续函数,那么分数阶微分方程
(1)
等价于以下的第二类非线性Volterra积分方程
(2)
2 算法
2.1 Runge-Kutta型方法简介
Runge-Kutta型方法为单步方法,针对方程
(3)
给出其算法.一般而言,若xn已知,则
xn+1-xn=hΦ(tn,xn,h),n=0,1,….
(4)
算法的核心在于构造增量函数Φ. Taylor级数展开法是用函数f在同一点(tn,xn)的高阶导数表示Φ(tn,xn,h),这不便于数值计算.Runge-Kutta型方法是用函数f在一些点上的值来表示Φ(tn,xn,h),使单步法的局部截断误差与Taylor级数展开法相等.我们针对方程(3),在[0,h]上取m个点0=t1≤t2≤t3≤…≤tm≤h,若知道ki=f(ti,x(ti)),i=1,2,…,m,则用它们的一次组合去近似f,
(5)
问题是如何计算ki,(因x(ti)未知).一个直观的想法是:设已知(t1,k1)=(t1,f(t1,x(t1)),由Eurler方法,x(t2)=x(t1)+(t2-t1)f(t1,x(t1))=x(t1)+(t2-t1)k1,于是
k2=f(t2,x(t1)+(t2-t1)k1).
(6)
同样,利用Eurler方法,又可以从(t2,k2)算出,
k3=f(t3,x(t1)+(t2-t1)k1+(t3-t2)k2).
(7)
如此可以继续进行下去,节点{ti}和系数ci可如此选择,使得近似式(5)有尽可能的逼近阶.
为便于推导,先引进若干记号.首先令ti=t1+aih,i=1,2,…,m.其中,ai与h无关.再引进如下的三角系数阵:
其中,bij与h无关,
(8)
又ci≥0,
(9)
{ai},{bij}和{ci}已给定,则Runge-Kutta法的计算过程如下:
un+1=un+hΦ(tn,un,h).
(10)
其中,
(11)
(12)
系数{ai},{bij}和{ci}按如下原则确定:将ki关于h展开,以之代到(11),使得l次幂hl(l=0,1,…,p-1)的系数和
(13)
中的同次幂的系数相等.如此得到的算法(10)称为m级p阶Runge-Kutta算法.
常用的Runge-Kutta算法为如下的4级4阶经典的Runge-Kutta算法.
(14)
2.2 Runge-Kutta型方法的稳定性与收敛性
定理1[17]设Φ(t,x,h)(t0≤t≤T,0≤h≤h0,x∈(-,+))关于u满足Lipschitz条件(Lipschitz常数为L),则单步法(10)稳定.
定理2[17]设Φ(t,x,h)满足定理1的条件,又单步法相容,则当h→0时,它的数值解xn→x(t),只要t0+nh→t,初值x0→x(t0).更设单步法(10)是p阶单步法,则还有收敛估计|en·≤eLT[(1+Lh)|e0·+cThp]. 设f(t,x)(t0≤t≤T,x∈(-,+))上连续,且关于u满足Lipschitz条件.由Φ(t,x(t),h)的表达式(12)知Φ(t,x,0)=f(t,x),即Runge-Kutta法相容,所以定理1和定理2对Runge-Kutta法恒成立.
2.3 分数阶微分方程的Runge-Kutta法构造
考虑如下的分数阶微分方程:
(15)
的Runge-Kutta方法的构造.其中,f(t,x(t))是[0,1]×R上的连续函数.
由引理3,问题(15)可以写成如下的积分形式:
(16)
步骤2:计算x1,
设xn-1已经计算出,下一步xn.
这样,我们就可以计算出方程(15)的数值解{xn}.