1、常微分方程初值问题数值解法综述常微分方程初值问题数值解法朱欲辉(浙江海洋学院 数理信息学院 , 浙江 舟山 316004)摘要 ：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法 . 然而在生产实际和科学研究中所遇到的问题往往很复杂 , 在很多情况下都不可能给出解的解析表达式 . 本篇文章详细介绍了常微分方程初值问题的一些数值方法 , 导出了若干种数值方法 , 如Euler法、改进的Euler法、Runge Kutta法以及线性多步法中的 Adams显隐式公式和预测校正公式 , 并且对其稳定性及收敛性作了理论分析 . 最后给出了数值例子 , 分别用不同的方法 计算出近似解 , 从得出
2、的结果对比各种方法的优缺点 .关键词：常微分方程 ; 初值问题; 数值方法 ; 收敛性; 稳定性; 误差估计Numerical Method for Initial-Value ProblemsZhu Yuhui(School of Mathematics, Physics, and Information Science,Zhejiang Ocean University, Zhoushan, Zhejiang 316004)Abstract: In the course about ordinary differential equations, the methods for analy
3、tic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can bte expressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods
4、 include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of the
5、oretical analysis.Keywords: Ordinary differential equation; Initial-value problem; Numerical method;Convergence; Stability ; Error estimate1刖言自然界和工程技术中的很多现象，例如自动控制系统的运行、电力系统的运行、飞行器 的运动、化学反应的过程、生态平衡的某些问题等，都可以抽象成为一个常微分方程初值问 题其真解通常难以通过解析的方法来获得 ，至今有许多类型的微分方程还不能给出解的解析表达式，一般只能用数值的方法进行计算 有关这一问题的研究早在十八世纪就已经
6、开始了，特别是计算机的普遍应用 ，许多微分方程问题都获得了数值解 ，从而能使人们认识解的种种性质及其数值特征 ，为工程技术等实际问题提供了定量的依据 关于常微分方程初值问题的数值计算方法 ，许多学者己经做了大量的工作 1768年,Euler提出了关于常微分方程初值问题的方法 ，1840年，Cauchy第一次对初值问题进行了仔细的分析，早期的常微分方程数值解的问题来源于天体力学 在1846年，当Adams还是一个学生的时候，和Le Verrier 起根据天王星轨道中出现的己知位置 ，预测了它下一次出现的位置.1883年，Adams 提出了 Adams Bashforth 和 Adams Moul
7、ton 方法.Rull （1895 年）、 Heun（1900年）和 Kutta （1901 年）提出 Runge.Kutta 方法.二十世纪五十年代，Dahlquist建立了常微分方程数值解法的稳定性理论 ，线性多步法是常微分方程初值问题的一种数值方法 .由于通常的数值方法，其绝对稳定区域是有限的，不适用于求解刚性常微分的初值问题 .刚性微分方程常常出现于航空、 航天、热核反应、自动控制、电子网络及化学动力学等一系列与国防和现代化建设密切相关的高科技领域 ，具有无容置疑的重要性.因此，刚性微分方程的研究工作早在二十世纪五十年代就开始了 ，1965年，在爱丁堡举行的IFIP会议后，更进一步地认
8、识刚性方程的普遍性和重要性 .自从六十年代初，许多数值分析家致力于探讨刚性问题的数值方法及其理论 ，注意到刚性问题对传统数值积分方法所带来的挑战 .这一时期，人们的研究主要集中在算法的线性稳定性上,就是基于试验方程 y，- C）数值解的稳定性研究.在此领域发表了大量的 论文，取得了许多重要的理论成果 .例如，1963年，Dahlquist给出A稳定性理论，1967年, Widlund给出A（）一稳定性理论，1969年，Gear将 A 稳定性减弱，给出刚性（Stiff）稳定性 理论，并找到了当k空6的k步k阶的刚性稳定方法，1969年Dill找到刚性稳定的7阶和8阶以 及1970年Jain找到刚
9、性稳定的9阶到11阶，但可用性没有检验.这些稳定性理论和概念都是 在线性试验方程的框架下推导出的 ，从严格的数学意义上来说，这些理论只适用于常系数线性自治系统.但从实用的观点来说，这些理论无疑是合理和必要的 ，对刚性问题的算法设计具有重要的指导意义 在八十至九十年代，国内也有一些学者研究线性理论 ，主要有匡蛟勋、陈果良、项家祥、李寿佛、黄乘明、李庆扬和费景高等线性理论虽然对一般问题具有指导作用 ，但其不能作为非线性刚性问题算法的稳定性理论研究基础为了将线性理论推广到非线性问题中 ，人们开始对非线性模型问题进行研究但是,早期文献主要致力于数值方法基于经典 Lipschitz条件下的经典收敛理论,
10、即认为良好的稳定性加上经典相容性和经典相容阶就足以描述方法的整体误差性态 直到1974年,Prothero和Robi nson首先注意到算法的经典误差估计由于受刚性问题巨大参数的影响而严重 失真，产生阶降低现象，这时人们认识到经典收敛理论对于非线性刚性问题以及线性模型 的不足于是，1975年，Dahlquist和Butcher分别提出了单支方法和线性多步法的 G一急定概念和B稳定概念.这两个概念填补了非线性稳定性分析理论 ，引起了计算数学家们的极大关注，在上述理论的基础上，1975年至1979年，Burrage和Butcher提出了 AN 一稳定性与BN 一 稳定性概念，并相应地建立了基本的
11、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）弱
12、代数稳定准则和多步Runge- Kutta法的一系列代数准则.此外，Dahquist, Butcher和Hairer分别深刻地揭示了单支 方法、一般线性方法和 Runge.Kutta方法线性与非线性稳定性之间的内在关系 .为了求解刚性微分方程，不少文献中构造含有稳定参数的线性多步方法 ，利用适当选择稳定参数来扩大方法的稳定区域.所有改进的思想，都是通过构造一些特殊的显式或隐式线性多步法 ，使其具有增大的稳定域，或使A（） 稳定的角增大.八十年代，就成为国内外学者所 研究的一个课题，学者主要有Rodabaugh和Thompson、Feinberg、李旺尧、李寿佛、包雪松、 徐洪义、刘发旺、匡蛟勋
13、、项家祥、蒋立新、李庆扬、谢敬东和李林忠等 .当前国内外研究刚性问题的一个主要趋势就是在 E一理论指导下寻找更为有效的新算法.另一个发展趋势就是力图突破单边lipschitz常数和内积范数的局限，建立比B理论更为普遍的定量分析收敛 理论.近年来，刚性延迟系统的算法研究成为刚性问题的另一个热点研究领域 ，张诚坚将Burrage等人创立的针对刚性常微分系统的 B理论拓展到非线性刚性延迟系统起来的单支方法、分块方法、循环方法、外推法、混合方法、二阶导数法以及各种常用的估校正算法.其中经常用到的线性多步法公式有 Euler公式、Heun公式、中点公式、Milne公式、Adams公式、simpson公式
14、、Hamming公式,Gear方法、Adams预估一校正法和 Mile预估 一 Hamming校正法公式等，此外还包含许多迄今尚末探明的新公式 .Burage曾将线性多步法和Runge Kutta法比作大海中的两座小岛，在浩瀚的汪洋之中，还有许多到现在没有发 现的新方法.本篇文章详细介绍了常微分方程初值问题的一些数值方法 ，导出了若干种数值方法 ，如Euler法、改进的欧拉法、显式龙格 -库塔法、隐式龙格-库塔法以及线性多步法中的 Adams显隐式公式和预测校正公式，并且对其稳定性及收敛性作了理论分析 .最后给出了数值例子，进行了计算机程序算法的分析与实现 ，以计算机的速度优势来弥补计算量大的
15、不足 .从得出的结果对比各种方法的优缺点 .2常微分方程初值问题的数值解法2.1数值方法的基本思想考虑一阶常微分方程初值问题dxf(x,y),x a,b,dyy(xc)=yo (2.1)的数值解法，其中f是x和y的己知函数，y0是给定的初始值.对于常微分方程初值问题 (2.1)数值解法，就是要算出精确解 y(x)在区间a,b上的一系列离散节点 x0 x I : xn4 xn = b的函数值y(x) , y(xj, Hl,y(xn)的近似值y0, %, ，yn.相邻两个节点的间距 xi称为步长，这时节点也可以表示为X ih (i =1,2, III,n).数值解法需要把连续性的问题加以离散化 ,
16、从而求出离散节点的数值解.通常微分方程初值问题(2.1)的数值方法可以分为两类(2)多步法-计算y(x)在X二Xn 1处的值需要应变量及其导数在 Xn 1之前的多个网个节点出的值.2.2 Euler 方法2.2.1 Euler 公式2.2.2 梯形公式欧拉法形式简单精度低,为了提高精度，对方程y = f (x, y)的两端在区间Xn,Xn 1上积分得改用梯形方法计算其积分项代入式(2.3),并用yn近似代替式中y(xn)即可得到梯形公式h(2.4)yn yn - f (Xn,n) f (Xn 1, Yn 1 )由于数值积分的梯形公式比矩形公式精度高 ，所以梯形公式(2.4)比欧拉公式(2.2)
17、的精度高一个数值方法.式(2.4)的右端含有未知的 yn1,它是关于yn d的函数方程，这类方法称为隐式方法223改进的欧拉公式梯形公式实际计算时要进行多次迭代 ，因而计算量较大在实用上，对于梯形公式(2.4)只迭代一次,即先用欧拉公式算出 yn1的预估值了宀，再用梯形公式(2.4)进行一次 迭代得到校正值yn勺，即+ yn 1 =yn hf(Xn,yn) h (25)| yn 1 = yn * f (Xn , yn ) f (xn 1, yn -1)2.2.4欧拉法的局部截断误差衡量求解公式好坏的一个主要标准是求解公式的精度 ，因此引入局部截断误差和阶数概念.定义2.1 在yn准确的前提下，
18、即yn二y(xn)时，用数值方法计算yn1的误差 尺1二y(Xn 1)- yn 1，称为该数值方法计算 yn 1时的局部截断误差对于欧拉公式，假定yn二y(xn)，则有yn 1 =y(Xn) hf(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)，则称这种数值方法的阶数是长(h1)越小，P越高，则局部截断误差越小，计算精度越高2.3 Runge Kutta 方法2.3.1Runge Kutta方法的基本思想欧拉公式可改写成
19、yn 1 = ynK f (Xn, yn)则yn 1的表达式与y(xn1)的泰勒展开式的前两项完全相同 ,即局部截断误差为 0().改进的欧拉公式又可改写成hy = yn十2(心+心)& = f (Xn 朴 yn +hkj1 L上述两组公式在形式上有一个共同点 ：都是用f(x,y)在某些点上值的线形组合得出y(X.,)的近似值yn d,而且增加计算f(x,y)的次数，可提高截断误差的阶如欧拉公式每步计算一次f(x, y)的值，为一阶方法改进的欧拉公式需计算两次 f (x, y)的值，它是二阶方法它的局部截断误差为 0(h3).于是可考虑用函数 f (x, y)在若干点上的函数值的线形组合来构造
20、近似公式 ，构造时要求近似公式在(xn, yn)处的泰勒展开式与解 y( x)在xn处的泰勒展开式的前面几项重合 ， 从而使近似公式达到所需要的阶数 既避免求偏导，又提高了计算方法精度的阶数 或者 说，在xx这一步内多预报几个点的斜率值 ,然后将其加权平均作为平均斜率 则可 构造出更高精度的计算格式，这就是Runge Kutta方法的基本思想232 Runge Kutta 方法的构造一般地，Runge Kutta方法设近似公式为 (下面的公式修改了) p%卅= yn+h瓦 qKi! 7“ K f (xn，yn) (1 =2,3，|，p)i 4Kj = f(焉 +ahyn + 运 bjKj)I.
21、 j其中aj，q，G都是参数，确定它们的原则是使近似公式在 ym处的泰勒展开式与 y(x)在Runge Kutta 公式xn处的泰勒展开式的前面的项尽可能多的重合 ，这样就使近似公式有尽可能高的精度 以此我们可以通过一个复杂的计算过程得出常用的的三阶和四阶yn = yn +：(心 +4K2 +K3)6Ki=f(Xn,yn)1 h hK2 = f(Xn +q，yn +君1)、K3 = f(Xn+h,ynhKi+2hK2)和 h*卑7.+6(心+2心+2心+心)Ki = f (Xn, yn)h hK2 = f (Xn + yn + Ki) (2.6)2 2h hK3 = f (Xn yn K2)2
22、 2K f (Xn h, yn hK3)式(2.6)称为经典 Runge Kutta方法.2.4线性多步法在逐步推进的求解过程中，计算yn1之前事实上已经求出了一系列的近似值 y0, %,yn.如果充分利用前面多步的信息来预测 ym,则可期望获得较高的精度 ，这就是构造多步法的基本思想.线性k步方法的一般公式为k 4 k4yn1=ajyn4 h bj f(Xn4,yn_j) (2.7)j=e j=J其中aj, bj均为与n无关的常数，|ak+bkj HO.当b_ =0时为显格式；当时1为隐格式.特别当 k = 1, a0 = b0 = 1,b = 0时为 Euler 公式；当 k = 1,a0
23、 =1,b0 =b_i = 2 时为梯形公式.定义2.3称kJ kJRn 1二y(Xn 1)yn 1八2 bj f (人小)j j =J为k步公式(2.7)在xn d处的局部截断误差当Rn l =0(hP )时称式(2.7)是p阶的.应用方程y (x)二f (x, y(x)可知局部截断误差也可写成为k A kJI Rn+ = y(Xn 卅)yn 卑二瓦 ajyn+h 瓦 bjy (Xn) j=0 j=J定义2.4如果线性k步方法(2.7)至少是1阶的，则称是相容的；如果线性k步法(2.7)是 p阶的，贝U称是p阶相容的.2.4.1Adams 外插法将微分方程生=f (x, y),x a,b,d
24、y的两端从xn到xn 1进行积分，得到xn +y(xnJ=y(x) x f(x,y(x)dx (2.8)xn我们用插值多项式代替右端的被积函数 .Adams外插法选取k个点xn, xn斗=111, xn乂十作为插值基点构造f (x, y)的k-1阶多项式Adams外插法的计算公式为kAYn Yn h a jfm7其中aj满足如下代数递推式j =0,1,hi根据此递推公式，可逐个的计算q (j =0,1,|),表2.1给出了 aj的部分数值表2.1j012345书* $aj11/25/123/8251/72095/288书* $2.4.2Adams 内插法根据插值理论知道，插值节点的选择直接影响
25、着插值公式的精度 ，同样次数的内插公式的精度要比外插公式的高 .仍假定已按某种方法求得问题 (2.1)的解y(x)在xXi =0,1,l|n)处的数值yi,并选取插 值节点xn4, III, XniH ,Xn+p P是正整数，用Lagrange型插值多项式LPiL(x)构造 y= f( x, y)可以导出解初值问题(2.1)的Adams内插公式：冷1丫风十)=y(Xn) + J Lnku(x)dx (2.9)Xn当p = 0时上式就退化为内插公式.内插公式(2.9)除了包含f在Xn1,lH,Xn处的已知值外，还包含了在点Xn,川，Xn.p,处的未知值.因此内插公式(2.9)只给出了未知量 yn
26、,IILyn.p的关系式，实际计算时仍需要解联立方程组 p =1的内插公式是最适用的 ，Lk j(x)采用Newt on向后插值公式得到Adams内插公式kyn yn aj j fmj=0其中系数a定义为j =0,1,HI其中a j满足如下代数递推式：_ 1,H0_ 0,j 0根据此递推公式，可逐个的计算a (j =0,1,1),表2.2给出了 a的部分数值表2.2j012345* * *a j1-1/2-1/12-1/24-19/720-3/160* * 2.5算法的收敛性和稳定性2.5.1相容性初值问题(2.1):f(x,y),x a,b, dyy() = y。的显式单步法的一般形式为yn
27、 .1 h ：(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)成立，则称单
28、步法(2.10)与微分方程初值问题(2.1)相容.容易验证，Euler法与改进Euler法均满足相容性条件.事实上，对Euler法，增量函数 为(x，y，h)二 f(x，y)自然满足相容性条件，改进Euler法的增量函数为1(x，y，h) f(x，y) f(x h，y hf(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阶精 度(P 0),则其局部误差为y(x h)-y(x) h(x,y(x),h)HO(hp1)上式两端同除以h,得；n勺
29、0令h 0 ,如果(x,y,h),则有y(x)_ ：(x,y,h)=O即(x.y.hH f(x, y)所以p 0的显示单步法均与初值问题 (2.1)相容由此各阶RK方法与初值问题(2.1)相容.2.5.2收敛性(x 二 a nh)定义2.6 一种数值方法称为是收敛的,如果对于任意初值 y0及任意(a,b,都有忸 yn = y(x)其中y(x)为初值问题(2.1)的准确解.按定义，数值方法的收敛性需根据该方法的整体截断误差 ；n -来判定.已知Euler方法的整体截断误差有估计式当h r 0 , ；n q 0,故Euler方法收敛.定理2.1设显式单步法具有p阶精度，其增量函数 (x, y, h
30、)关于y满足Lipschitz条件,问题(2.1)是精确的，既y() ，则显式单步法的整体截断误差为；n 1 =y(Xn 1) - 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)对于xa,b, h (0,h0以及一切实数y,关于y满足Lipschitz条件，则单步法(2.7)是稳定的.定义2.8对给定的微分方程和给定的步长 h,如果由单步法计算 yn时有大小为:.的误差，即计算得出yn二yn.，而引起其后之值 ym(mn)的变化小于,则说该单步法 是绝对稳定的一般只限于典型微分方
