常微分方程初值问题数值解法综述.docx
《常微分方程初值问题数值解法综述.docx》由会员分享，可在线阅读，更多相关《常微分方程初值问题数值解法综述.docx（25页珍藏版）》请在冰豆网上搜索。
常微分方程初值问题数值解法综述
常微分方程初值问题数值解法
朱欲辉
（浙江海洋学院数理信息学院,浙江舟山316004）
[摘要]：
在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实
际和科学研究中所遇到的问题往往很复杂,在很多情况下都不可能给出解的解析表达式.本篇文章详细介绍了常微分方程初值问题的一些数值方法,导出了若干种数值方法,如
Euler法、改进的Euler法、Runge—Kutta法以及线性多步法中的Adams显隐式公式和预测校正
公式,并且对其稳定性及收敛性作了理论分析.最后给出了数值例子,分别用不同的方法计算出近似解,从得出的结果对比各种方法的优缺点.
[关键词]：
常微分方程;初值问题;数值方法;收敛性;稳定性;误差估计
NumericalMethodforInitial-ValueProblems
ZhuYuhui
（SchoolofMathematics,Physics,andInformationScience,
ZhejiangOceanUniversity,Zhoushan,Zhejiang316004）
[Abstract]:
Inthecourseaboutordinarydifferentialequations,themethodsforanalyticsolutionsofsometypicalequationsareoftendiscussed.However,inscientificresearch,theproblemsareverycomplexandtheanalyticsolutionsabouttheseproblemscan'bteexpressedexplicitly.Inthispaper,somenumericalmethodsfortheinitial-valueproblemsareintroduced.thesemethodsincludeEulermethod,improvedEulermethod,Runge-Kuttamethodandsomelinearmultistepmethod（e.g.Adamsformulaandpredicted-correctedformula）.Thestabilityandconvergenceaboutthemethodsarepresented.Somenumericalexamplesaregivetodemonstratetheeffectivenessandaccuracyoftheoreticalanalysis.
[Keywords]:
Ordinarydifferentialequation;Initial-valueproblem;Numericalmethod;
Convergence;Stability;Errorestimate
1刖言
自然界和工程技术中的很多现象，例如自动控制系统的运行、电力系统的运行、飞行器的运动、化学反应的过程、生态平衡的某些问题等，都可以抽象成为一个常微分方程初值问题•其真解通常难以通过解析的方法来获得，至今有许多类型的微分方程还不能给出解的
解析表达式，一般只能用数值的方法进行计算•有关这一问题的研究早在十八世纪就已经
开始了，特别是计算机的普遍应用，许多微分方程问题都获得了数值解，从而能使人们认
识解的种种性质及其数值特征，为工程技术等实际问题提供了定量的依据•
关于常微分方程初值问题的数值计算方法，许多学者己经做了大量的工作•1768年,
Euler提出了关于常微分方程初值问题的方法，1840年，Cauchy第一次对初值问题进行了仔
细的分析，早期的常微分方程数值解的问题来源于天体力学•在1846年，当Adams还是一
个学生的时候，和LeVerrier—起根据天王星轨道中出现的己知位置，预测了它下一次出现
的位置.1883年，Adams提出了Adams—Bashforth和Adams—Moulton方法.Rull（1895年）、Heun（1900年）和Kutta（1901年）提出Runge.Kutta方法.
二十世纪五十年代，Dahlquist建立了常微分方程数值解法的稳定性理论，线性多步法
是常微分方程初值问题的一种数值方法.由于通常的数值方法，其绝对稳定区域是有限的，
不适用于求解刚性常微分的初值问题.刚性微分方程常常出现于航空、航天、热核反应、自
动控制、电子网络及化学动力学等一系列与国防和现代化建设密切相关的高科技领域，具有
无容置疑的重要性.因此，刚性微分方程的研究工作早在二十世纪五十年代就开始了，
1965年，在爱丁堡举行的IFIP会议后，更进一步地认识刚性方程的普遍性和重要性.自从
六十年代初，许多数值分析家致力于探讨刚性问题的数值方法及其理论，注意到刚性问题
对传统数值积分方法所带来的挑战.这一时期，人们的研究主要集中在算法的线性稳定性
上,就是基于试验方程y，-C）数值解的稳定性研究.在此领域发表了大量的论文，取得了许多重要的理论成果.例如，1963年，Dahlquist给出A稳定性理论，1967年,Widlund给出A（〉）一稳定性理论，1969年，Gear将A—稳定性减弱，给出刚性（Stiff）稳定性理论，并找到了当k空6的k步k阶的刚性稳定方法，1969年Dill找到刚性稳定的7阶和8阶以及1970年Jain找到刚性稳定的9阶到11阶，但可用性没有检验.这些稳定性理论和概念都是在线性试验方程的框架下推导出的，从严格的数学意义上来说，这些理论只适用于常系数
线性自治系统.但从实用的观点来说，这些理论无疑是合理和必要的，对刚性问题的算法
设计具有重要的指导意义•在八十至九十年代，国内也有一些学者研究线性理论，主要有
匡蛟勋、陈果良、项家祥、李寿佛、黄乘明、李庆扬和费景高等
线性理论虽然对一般问题具有指导作用，但其不能作为非线性刚性问题算法的稳定性
理论研究基础•为了将线性理论推广到非线性问题中，人们开始对非线性模型问题进行研
究•但是,早期文献主要致力于数值方法基于经典Lipschitz条件下的经典收敛理论,即认为
良好的稳定性加上经典相容性和经典相容阶就足以描述方法的整体误差性态•直到1974年,
Prothero和Robinson首先注意到算法的经典误差估计由于受刚性问题巨大参数的影响而严重失真，产生阶降低现象，这时人们认识到经典收敛理论对于非线性刚性问题以及线性模型的不足•于是，1975年，Dahlquist和Butcher分别提出了单支方法和线性多步法的G一急定概
念和B稳定概念.这两个概念填补了非线性稳定性分析理论，引起了计算数学家们的极大关
注，在上述理论的基础上，1975年至1979年，Burrage和Butcher提出了AN一稳定性与BN一稳定性概念，并相应地建立了基本的B一稳定及代数稳定理论.1981至1985年，Frank,
schneid和ueberhuber建立Runge一kutta方法的旷收敛理论.B稳定与B一收敛理论统称B—理论，它是常微分方程数值解法研究领域的巨大成就之一，是刚性问题算法理论的突破性进
展，标志着刚性问题研究从线性向非线性情形深入发展.国内也有众多学者致力于B一理
论的研究，如李寿佛、曹学年等.1989年，李寿佛将Dahlquist的G—稳定概念推广到更一般的（C,P,Q）代数稳定，克服了G急定的线性多步法不能超过二阶的限制.对于一般线性方法，
李寿佛建立了一般线性方法的（K,P,Q）稳定性理论及（K,P,Q）弱代数稳定准则和多步
Runge-Kutta法的一系列代数准则.此外，Dahquist,Butcher和Hairer分别深刻地揭示了单支方法、一般线性方法和Runge.Kutta方法线性与非线性稳定性之间的内在关系.为了求解刚
性微分方程，不少文献中构造含有稳定参数的线性多步方法，利用适当选择稳定参数来扩
大方法的稳定区域.所有改进的思想，都是通过构造一些特殊的显式或隐式线性多步法，
使其具有增大的稳定域，或使A（＞）—稳定的〉角增大.八十年代，就成为国内外学者所研究的一个课题，学者主要有Rodabaugh和Thompson、Feinberg、李旺尧、李寿佛、包雪松、徐洪义、刘发旺、匡蛟勋、项家祥、蒋立新、李庆扬、谢敬东和李林忠等.当前国内外研究
刚性问题的一个主要趋势就是在E一理论指导下寻找更为有效的新算法.另一个发展趋势就
是力图突破单边lipschitz常数和内积范数的局限，建立比B—理论更为普遍的定量分析收敛理论.近年来，刚性延迟系统的算法研究成为刚性问题的另一个热点研究领域，张诚坚将
Burrage等人创立的针对刚性常微分系统的B—理论拓展到非线性刚性延迟系统
起来的单支方法、分块方法、循环方法、外推法、混合方法、二阶导数法以及各种常用的估
校正算法.其中经常用到的线性多步法公式有Euler公式、Heun公式、中点公式、Milne公
式、Adams公式、simpson公式、Hamming公式,Gear方法、Adams预估一校正法和Mile预估一Hamming校正法公式等，此外还包含许多迄今尚末探明的新公式.Burage曾将线性多步
法和Runge—Kutta法比作大海中的两座小岛，在浩瀚的汪洋之中，还有许多到现在没有发现的新方法.
本篇文章详细介绍了常微分方程初值问题的一些数值方法，导出了若干种数值方法，
如Euler法、改进的欧拉法、显式龙格-库塔法、隐式龙格-库塔法以及线性多步法中的Adams
显隐式公式和预测校正公式，并且对其稳定性及收敛性作了理论分析.最后给出了数值例
子，进行了计算机程序算法的分析与实现，以计算机的速度优势来弥补计算量大的不足.
从得出的结果对比各种方法的优缺点.
2常微分方程初值问题的数值解法
2.1数值方法的基本思想
考虑一阶常微分方程初值问题
dx
f（x,y）,x[a,b],
dy
y（xc）=yo（2.1）
的数值解法，其中f是x和y的己知函数，y0是给定的初始值.
对于常微分方程初值问题（2.1）数值解法，就是要算出精确解y（x）在区间[a,b]上的一
系列离散节点x0x^I:
xn4xn=b的函数值y（x））,y（xj,Hl,y（xn）的近似
值y0,%,•…，yn.相邻两个节点的间距xi称为步长，这时节点也可以表示为
Xih（i=1,2,III,n）.数值解法需要把连续性的问题加以离散化,从而求出离散节
点的数值解.
通常微分方程初值问题（2.1）的数值方法可以分为两类
（2）多步法-----计算y（x）在X二Xn1处的值需要应变量及其导数在Xn1之前的多个网个节
点出的值.
2.2Euler方法
2.2.1Euler公式
2.2.2梯形公式
欧拉法形式简单精度低,为了提高精度，对方程y=f（x,y）的两端在区间
[Xn,Xn1]上积分得
改用梯形方法计算其积分项
代入式（2.3）,并用yn近似代替式中y（xn）即可得到梯形公式
h
（2.4）
ynyn-[f（Xn,『n）f（Xn1,Yn1）]
由于数值积分的梯形公式比矩形公式精度高，所以梯形公式（2.4）比欧拉公式（2.2）的精
度高一个数值方法.
式（2.4）的右端含有未知的yn1,它是关于ynd的函数方程，这类方法称为隐式方法
223改进的欧拉公式
梯形公式实际计算时要进行多次迭代，因而计算量较大•在实用上，对于梯形公式
（2.4）只迭代一次,即先用欧拉公式算出yn1的预估值了宀，再用梯形公式（2.4）进行一次迭代得到校正值yn勺，即
+yn1=ynhf（Xn,yn）
|yn1=yn*^[f（Xn,yn）f（xn1,yn-1）]
2.2.4欧拉法的局部截断误差
衡量求解公式好坏的一个主要标准是求解公式的精度，因此引入局部截断误差和阶数
概念.
定义2.1在yn准确的前提下，即yn二y（xn）时，用数值方法计算yn1的误差尺1二y（Xn1）-yn1，称为该数值方法计算yn1时的局部截断误差•
对于欧拉公式，假定yn二y（xn），则有
yn1=y（Xn）h[f（Xn，y（Xn））Hy（Xn）hy'（Xn）
而将真解y（X）在Xn处按二阶泰勒展开式有
h2
Yn^y（Xn）hy'（Xn）=2（），X，X.1）
2!
因此有
P.步
定义2.2若数值方法的局部截断误差为0（hp1），则称这种数值方法的阶数是
长（h<1）越小，P越高，则局部截断误差越小，计算精度越高
2.3Runge—Kutta方法
2.3.1Runge—Kutta方法的基本思想
欧拉公式可改写成
yn1=yn
K^f（Xn,yn）
则yn1的表达式与y（xn1）的泰勒展开式的前两项完全相同,即局部截断误差为0（『）.
改进的欧拉公式又可改写成
h
y^=yn十2（心+心）
&=f（Xn朴yn+hkj
1L
上述两组公式在形式上有一个共同点：
都是用f（x,y）在某些点上值的线形组合得出
y（X.,）的近似值ynd,而且增加计算f（x,y）的次数，可提高截断误差的阶•如欧拉公式
每步计算一次f（x,y）的值，为一阶方法•改进的欧拉公式需计算两次f（x,y）的值，它是
二阶方法•它的局部截断误差为0（h3）.
于是可考虑用函数f（x,y）在若干点上的函数值的线形组合来构造近似公式，构造时
要求近似公式在（xn,yn）处的泰勒展开式与解y（x）在xn处的泰勒展开式的前面几项重合，从而使近似公式达到所需要的阶数•既避免求偏导，又提高了计算方法精度的阶数•或者说，在［x^x^］这一步内多预报几个点的斜率值,然后将其加权平均作为平均斜率•则可构造出更高精度的计算格式，这就是Runge—Kutta方法的基本思想•
2•3•2Runge—Kutta方法的构造
一般地，Runge—Kutta方法设近似公式为（下面的公式修改了）
'p
%卅=yn+h瓦qKi
!
7
“K^f（xn，yn）（1=2,3，|||，p）
i4
Kj=f（焉+ahyn+运bjKj）
I.j^
其中aj，q，G都是参数，确定它们的原则是使近似公式在ym处的泰勒展开式与y（x）在
Runge—Kutta公式
xn处的泰勒展开式的前面的项尽可能多的重合，这样就使近似公式有尽可能高的精度•以
此我们可以通过一个复杂的计算过程得出常用的的三阶和四阶
yn^=yn+：
（心+4K2+K3）
6
Ki=f（Xn,yn）
1hh
K2=f（Xn+q，yn+君1）
、K3=f（Xn+h,yn—hKi+2hK2）
和
•h
*卑7.+6（心+2心+2心+心）
Ki=f（Xn,yn）
hh
K2=f（Xn+yn+Ki）（2.6）
22
hh
K3=f（XnynK2）
22
K^f（Xnh,ynhK3）
式（2.6）称为经典Runge—Kutta方法.
2.4线性多步法
在逐步推进的求解过程中，计算yn1之前事实上已经求出了一系列的近似值y0,%,
yn.如果充分利用前面多步的信息来预测ym,则可期望获得较高的精度，这就是
构造多步法的基本思想.
线性k步方法的一般公式为
k4k4
yn1=»ajyn4h'bjf（Xn4,yn_j）（2.7）
j=ej=J
其中aj,bj均为与n无关的常数，|ak』+bkjHO.当b_^=0时为显格式；当时
1
为隐格式.特别当k=1,a0=b0=1,b^=0时为Euler公式；当k=1,a0=1,b0=b_i=2时
为梯形公式.
定义2.3称
kJkJ
Rn1二y（Xn1）—[yn1八2bjf（人」小」）]
j」j=J
为k步公式（2.7）在xnd处的局部截断误差•当
Rnl=0（hP）
时称式（2.7）是p阶的.
应用方程y（x）二f（x,y（x））可知局部截断误差也可写成为
kAkJ
IRn+=y（Xn卅）—[yn卑二瓦ajyn」+h瓦bjy（Xn」）]j=0j=J
定义2.4如果线性k步方法（2.7）至少是1阶的，则称是相容的；如果线性k
步法（2.7）是p阶的，贝U称是p阶相容的.
2.4.1Adams外插法
将微分方程
生=f（x,y）,x[a,b],
dy
的两端从xn到xn1进行积分，得到
xn+
y（xnJ=y（x）xf（x,y（x））dx（2.8）
xn
我们用插值多项式代替右端的被积函数.
Adams外插法选取k个点xn,xn斗=111,xn乂十作为插值基点构造f（x,y）的k-1阶多项式
Adams外插法的计算公式为
kA
Yn^Ynha^jfm
7
其中aj满足如下代数递推式
j=0,1,hi
根据此递推公式，可逐个的计算q（j=0,1,|||）,表2.1给出了aj的部分数值
表2.1
j
0
1
2
3
4
5
书*$
aj
1
1/2
5/12
3/8
251/720
95/288
书*$
2.4.2Adams内插法
根据插值理论知道，插值节点的选择直接影响着插值公式的精度，同样次数的内插公
式的精度要比外插公式的高.
仍假定已按某种方法求得问题（2.1）的解y（x）在xXi=0,1,l||n）处的数值yi,并选取插值节点xn」4,III,XniH,Xn+pP是正整数，用Lagrange型插值多项式LPiL（x）构造y=f（x,y）可以导出解初值问题（2.1）的Adams内插公式：
冷1
丫风十）=y（Xn）+JLnku（x）dx（2.9）
Xn
当p=0时上式就退化为内插公式.
内插公式（2.9）除了包含f在Xn»1,lH,Xn处的已知值外，还包含了在点Xn,川，Xn.p,
处的未知值.因此内插公式（2.9）只给出了未知量yn,IILyn.p的关系式，实际计算时仍需要
解联立方程组•p=1的内插公式是最适用的，L^kj（x）采用Newton向后插值公式得到
Adams内插公式
k
ynynajjfm
j=0
其中系数a］定义为
j=0,1,HI
其中aj满足如下代数递推式：
_1,H0
_0,j0
根据此递推公式，可逐个的计算a］（j=0,1,1"）,表2.2给出了a］的部分数值
表2.2
j
0
1
2
3
4
5
**■
*
aj
1
-1/2
-1/12
-1/24
-19/720
-3/160
**■
2.5算法的收敛性和稳定性
2.5.1相容性
初值问题（2.1）:
f（x,y）,x[a,b],dy
y（«）=y。
的显式单步法的一般形式为
yn.1h：
（Xn,yn,h）,n=0,1川|,N-1（2.10）
b_a
其中h二〒，人"nh.这样我们用差分方程初值问题⑵1。
）的解yn作为问题⑵1）的解y（xj在x=xn处的近似值，即yX）:
“yn•因此，只有在问题（2.1）的解使得
y（X[胁-心恥）」）
h
逼近
y‘—f（x,y（x））=0
时，才有可能使（2.10）的解逼近于问题（2.1）的解.从而，我们期望对任一固定的x[a,b],都有
hmJy（x+hhW（x,y（x）,h）]=0
假设增量函数®（x,y,h）关于h连续，则有
（x,y，h）二f（x，y）
定义2.5若关系式
（x，y，h）=f（x，y）
成立，则称单步法（2.10）与微分方程初值问题（2.1）相容.
容易验证，Euler法与改进Euler法均满足相容性条件.事实上，对Euler法，增量函数为
（x，y，h）二f（x，y）
自然满足相容性条件，改进Euler法的增量函数为
1
（x，y，h）[f（x，y）f（xh，yhf（x，y））]
2
因为f（x,y）,从而有
1
（x,y,h^-[f（x,y）f（x,y）]=f（x,y）
2
所以Euler法与改进Euler法均与初值问题（2.1）相容.一般的，如果显示单步法有p阶精度（P0）,则其局部误差为
y（xh）-[y（x）h（x,y（x）,h）HO（hp1）
上式两端同除以h,得；n勺0令h—0,如果（x,y,h）,则有
y（x）_：
（x,y,h）=O
即
（x.y.hHf（x,y）
所以p0的显示单步法均与初值问题（2.1）相容•由此各阶RK方法与初值问题（2.1）相容.
2.5.2收敛性
（x二anh）
定义2.6一种数值方法称为是收敛的,如果对于任意初值y0及任意（a,b],都有
忸yn=y（x）
其中y（x）为初值问题（2.1）的准确解.
按定义，数值方法的收敛性需根据该方法的整体截断误差；n-来判定.已知Euler方法
的整体截断误差有估计式
当hr0,；nq0,故Euler方法收敛.
定理2.1设显式单步法具有p阶精度，其增量函数（x,y,h）关于y满足Lipschitz条件,
问题（2.1）是精确的，既y（^），则显式单步法的整体截断误差为
；n1=y（Xn1）-yn产。
伸卩）.
定理2.2设增量函数（x,y,h）在区域S上连续，且关于y满足Lipschitz条件时，则
显式单步法收敛的充分必要条件是相容性条件成立，即
：
（x,y,h）二f（x,y）.
2.5.3稳定性
定义2.7如果存在正常数h0及C,使得对任意的初始出发值y°,y。
，单步法（2.10）的
相应精确解yn,yn,对所有的0”：
h空h0,恒有
y。
一yo
则称单步法是稳定的•
定理2.3若（x,y,h）对于x・[a,b],h（0,h0]以及一切实数y,关于y满足Lipschitz
条件，则单步法（2.7）是稳定的.
定义2.8对给定的微分方程和给定的步长h,如果由单步法计算yn时有大小为:
.的误
差，即计算得出yn二yn〔.，而引起其后之值ym（m・n）的变化小于,则说该单步法是绝对稳定的•
一般只限于典型微分方
升级会员 