常微分方程数值解法的误差分析教材Word格式.docx

1、（一） 、欧拉法 5（二） 、向后Euler方法 6（三） 、0-法 7（四） 、改进欧拉法 7（五） Rung* Kutta 方法 92、 线性多步法 14总结 16参考文献： 17自然界中很多事物的运动规律可用微分方程来刻画。常微分方程是研究自然 科学和社会科学中的事物、物体和现象运动、演化和变化规律的最为基本的数学 理论和方法。物理、化学、生物、工程、航空航天、医学、经济和金融领域中的 许多原理和规律都可以描述成适当的常微分方程，因此，常微分方程的理论和方 法不仅广泛应用丁自然科学，而且越来越多的应用丁社会科学的各个领域。由丁 该问题比较复杂且涉及的面广，使得有些问题的解析解很难求出，而

2、对丁一些典 型的微分方程（如线性方程、某些特殊的一阶非线性方程等）可以运用基本方法求 出其解析解，并在理论上可以根据初值问题的条件把其中的任意常数完全确定下 来。然而，在生产实际和科学研究中所遇到的微分方程往往很复杂，在很多情况 下都不可能给出解的解析表达式，有时即使能求出形式的解，也往往因计算量太 大而不实用，而且高次代数方程求根也并不容易，所以用求解析解的方法来计算 微分方程的数值解往往是不适宜的。从实际意义来讲我们更关心的是某些特定的 自变量在某一个定义范围内的一系列离散点上的近似值。本文研究的主要是针对 常微分方程各种数值解法的误差进行分析。一、常微分方程1、定义首先，我们在这部分给出

3、所需的一些基本概念和基本知识。我们已经知道微分方程就是联系着自变量、未知函数以及其导数的关系式。 如果在微分方程中，自变量的个数只有一个，我们称这种微分方程为常微分方程。方程变y b曳 cy=f t dt2 dt就是常微分方程的例子，这里y是未知函数，t是自变量。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。2、常微分方程初值问题描述在自然科学和经济的许多领域中，常常会遇到一阶常微分方程的初值问题dy=f x, ydx yj,a4xb （1）yg ） = y（）这里f（x, y）是充分光滑，即关丁 x或y满足李普希茨条件的二元函数，y0是 给定的初值，y（xo ）=y称为初始条件

4、。3、数值解法的基本思想与途径一阶微分方程的初值问题（1）的解y（x）是区间Ia,b上的连续变量x的函数，因 而问题（1）实际上是一个连续性的问题，求这个问题的数值解，就是要求在区间 la,b上的若干个离散点处的函数近似值，例如：ax0 为 . xnb ,然后计算出解y（x ）的近似值y x。,y xi ,.,y xn .股常取xo,xi,.,xn为等距离的点，即为一 x。=x2 -为=.=xn xn= h或称h为步长xi = a i,h 0,1, . n,建立数值方法的第1步，就是把连续性问题（1）通过一定的方法化为在给定的 n+1个点上的近似的差分方程的初值问题，称这个过程为离散化。常用离

5、散化的 方法如下：（1）用差商替代导数在点x处的导数yxi）可以近似地表示成差商y（x 片从而把初值问题（1）化为差分问题y Xo =ys其中表小解 y（x）在点x处的近似解，即yi = y（ x ）。当然，用差商来近似地表示导数，方法不是唯一的，这里所用的是所谓的向 前差商。（二）Taylor展开法在一点（例如点x）的附近，y（x）的同次数的近似多项式中的 Taylor多项式py（Xi +h ）女 y（x ）+hy（为）+.+土 yp xi p!为最好。其中p为一正整数。通过微分方程 y=f（x,y）,便可以逐次把各阶 导数y；y,IH在x处的值表示出来。（三）数值积分法对微分方程y=f（x

6、,y侄区问仅,X上求积分，得xi 1y（x. ）y（xi ）=七 f（x,y（x）px, j = 0 ,“丁是，初值I可题（1）便可以近似地化为y i = yy xo = y。,xi f x,y x dxi=0,1川 n.这样，关丁上式右端的积分，可以用数值积分方法计算其近似值。4、数值解的分类常微分方程初值问题的数值解法一般分为两大类：单步法：所谓单步法是指这类方法在计算yn中时，只用到前一步的值 x 4i, xn, yn然后逐步往下计算。这个算法的代表是龙格-库塔算法，简称R k 方法。四阶显示Runge-Kutta方法是求解普通常微分方程初值问题数值解法中的重 要方法，而隐式Runge-

7、Kutta公式是求解刚性常微分方程初值问题的重要方法。多步法：这类方法在计算yn+时，除了用到前一步的值xn书,xn,yn，之外，还要 用到xn* y p =1,2,HI,k;k 0这前面k步的值，这个算法的代表就是阿达姆斯（Adams）方法。5、问题（1）解的存在惟一性定理一个常微分方程是不是有特解呢？如果有，乂有几个呢？这是微分方程论中一个 基本的问题，数学家把它归纳成基本定理，叫做存在和唯一性定理。因为如果没 有解，而我们要去求解，那是没有意义的；如果有解而乂不是唯一的，那乂不好 确定。因此，存在和唯一性定理对丁微分方程的求解是十分重要的。这个重要的 存在和唯一性就是下面列出的著名的存在

8、惟一性定理。定理如果f （ x, y ）在带形区域R = （x, y）|a壬xb, -* y+=中连 续，且关丁 y满足Lipchiz条件：即存在正常数L,使得f（ x y） f ,x2 0，然后在x0x x0 + h情况下用解函数的切线l x = y x- x f ,0 x 0 y代替解函数。这样对丁点 x1 = x0 + h就可以得到y1 = y 0 hf x y 0在点（xi, y1 ）重复如上的程序再次计算新的方向，就会得到所谓的递推公式：xm41 = xm + h , ymF = ym + hf（xm,ym）,这就是Euler方法。由此，再通过连接所有这些切线得到的函数被称为 Eul

9、er折线。如果我们令h- 0,这些折线就会越来越接近解函数。（一）、欧拉法Euler方法是最简单的一步法，它是一阶的，精度较差，但公式很简单，即yn 1 = yn hy Xn,yn n = 0 , 1 j。 （3）Euler方法的几何意义在数值计算思想中已经体现出来了， 实际上就是用过已知点的折线来近似代替过此点的积分曲线。 因此，这种方法乂称为折线法。在Euler法中，数值解的误差首先是由差商代替导数引起的，这种近似替代所产生的误差 称为截断误差。另外，计算过程中还会由丁数值的舍入产生另一种误差一一舍入 误差。显然只有当初产生的误差在以后各步的计算中不会无限制扩大时，即当初 始误差充分小时，

10、以后各步的误差也可以充分小， Euler法才具有实用价值。收敛性、截断误差估计与稳定性闷题是常微分方程各种数值解法研究中必须考虑的基 本问题。显然这些问题在 Euler法中是得到验证的，详见下面例子分析。在xy平面上，微分方程（1）的解y=y（x）称作它的积分曲线。积分曲线上一点（x,y） 线斜率等丁函数f（x,y）的值，如果按函数f（x, y）在xy平面上建立一个方向场，那 么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致，基丁上述几 何解释，从初始点Po（x0,yo）出发，先依方向场在该点的方向推进到x = x上一点臼， 然后再从p1依方向场的方向推进到 x = x2上一点p2

11、,循此前进推出一条折线P0P1P2，一般地，设已做出该折线的顶点 Pn，过Pn（X,yn）依方向场的方向再推 进到Pn41（xn+,yn+），显然两个顶点Pn , Pn*的坐标有关系yn 1 = yn hf xn, yn这就是著名的欧拉（Euler）公式。若初值y已知，则依公式可逐步算出V1=* hf x,yo ,V2=V hf *,乂 ,III例1求解初值问题y = y _2x/ y, （0 ： x : 1）y（0） 二1,解欧拉公式的具体形式为取步长h =0.1 ,计算结果如下表:表1计算结果对比初值问题（2.1.2）有解y =5/芯，按这个解析式子算出的准确值 y（xn）同近似值 yn

12、一起歹0在表1,两者相比较可以看出欧拉方法的精度很差。（二）、向后EuIer方法向后Euler方法和Euler方法差不多，只是把y（x）用y Xn 1 - y Xnh去代替，这时计算公式为（6）yn 1 = y hf xn1,yn1 , NJ（n =0,1,2 川）y x = y向后Euler方法的总体截断误差也是一阶的，因此向后 Euler方法是收敛的。这 里需要指出它与Euler方法的一个很大不同之处，Euler方法是显式方法，即由xn41,xn,yn,h明显地表示出来了，而向后 Euler方法是隐式方法，计算yn41时 要解隐式方程（6）。通常解此方程用迭代法。因此计算较为麻烦，但比显式

13、 Euler 方法精度要高。（三）、8-法将Euler方法公式与向后Euler方法公式作加权平均，得到如下公式：fyJ = yn +h 丐 f（xn, yn 广（1 顼）f （x忡，）。刃巾 =0,|,2 ）y 冷 ）=%,称式 为初值问题 的8 -法公式。其中如y（x0 ）= y0为初值条件。 1 -在此法中当6 =-时，即2 hyn = yn +2（xn, Yn ）+ f （ x” 甲，Yn* ） （8）此时的e一法称为梯形公式法。梯形公式也隐式格式，用起来要进行迭代，其计算 公式为,y）=mn,yn）+f（XmynQ）k = 0,i,2 川） ynJ）= yn +hf （xn, yn ）

14、,这里在应用本迭代法时，是先用 Euler方法求初值yn+1的近似值yn+i（）即： yn+i（Lyn + hf （Xn,yn ）然后将yn+1（陪代梯形公式（8）中的y n +1得到的式（9）。式乂称为预测校正公式。换言之，由 Euler方法给出预测值，再用梯形法 予以校正。很显然，当步长办取得适当小时，由 Euler方法算出的值已是较好的 近似。格式（9）收敛很快，通常只需一两次迭代即可满足精度要求， 若需多次迭代,则应缩小步长h后再行计算。梯形公式法比用向后Euler方法的迭代步长可以放宽一倍，它的总体截断误 差为O（h2 ），比Euler方法高一阶。但它每积分一步要计算二次函数值，这说

15、明 的精度的提高是以增加计算量为代价的。（4）xn中积分，得（10）、改进欧拉法为得到比欧拉法精度高的计算公式，在等式如果对方程 （1）从xn到xn 1y xn 1 = y xn f t,y t dtxn右端积分中若用梯形求积公式近似，并用 代替y（Xn） , yn书代替y（XnG，则得 h .yn 1 = yn 2 IIXn,yn f Xn 1, y” 1 （11）称为改进欧拉法.改进欧拉方法是隐式单步法，可用迭代法求解.用欧拉方法提供迭代初值，则改 进欧拉法的迭代公式为 乂拦=人 +；/（吊）+了0心丈2）（12）侬=。,成为了分析迭代过程的收敛性，将（11）式与（10）相减，得= 5 仇

16、 f m+l） -/（w，照） M于是有Am-欢倍竽*照式中L为f （x, y）对y满足Lipschitz常数，如果选取h充分小，使得hL *12 ,则当kT 8时有y：%T外十，这说明迭代过程（12）是收敛的.例2用改进的欧拉方法求解初值问题（1）.解改进的欧拉公式为* =人 + 第（L）+/g璀）2（上=。,12 ）.仍取h =0.1，计算结果见下表.同例1中欧拉法的计算结果比较，改进欧拉法明 显改善了精度.表2计算结果对比%旭）所）|0.10.20.30.40.51.09591.18411.26621.34341.41641.09541.18321.26491.34161.41420.6

17、0.70.80.91.01.48601.55251.61531.67821.73791.48321.54921.61251.67331.7321（五）RungKutta 方法德国数学家C. D. T. Runge（185 1927）是数值方法发展史上具有里程碑作用的人物。1895年，他在Hanover发表了关丁微分方程数值解法的经典论文 常微分方程数值解法。此文成为常微分方程RungaKutta方法的发端。此后，Runge 结合教学活动积极投身丁发展一般的数值分析特别是各种实际应用中的 Runge-Kutta 方法（严格来说，此方法在Kutta作出工作后才能称作 RungbKutta方法）。

18、RungbKutta方法是一种特殊的单步方法，事实上，这个方法可以看作在（Xm,XmG 上取若干条积分曲线的若干个点的切线斜率，再进行一次 （或多次）算术（或加权）平均后产生的新斜率，再按这个斜率从 （Xm,ym）出发，以直线带曲线向前推进一步的过程。与Taylor展示法相比，RungbKutta方法不用增加微商f（x,y）的次数就 可以得到较高的阶。Runge Kutta方法除了在微分方程求解中扮演的传统角色外， 人们发现相关类型的初值问题可以用 Runge 一 Kutta方法或适合更一般问题的RungbKutta方法求解，比如 RungbKutta方法被应用到了 Hamilton系统中。前

19、面提到的几种数值解法的精度是很低的，下面给出高阶一步法一一 Runge- Kutta 方法。它是最常用的一种数值解法，因为它相当精确、稳定、容易编程。Runge Kutta方法至今仍然得到广泛地应用。1、二级二阶 RungbKutta方法由RungbKutta方法的思想，我们得到二阶 RungbKutta公式为=yn+4（Ki+3K2）Ki =hf Xn,yn ,（ 2 2 ）（13）K2 =hf . Xn + h, yn + Ki .3 3yn* = yn +*Ki + K2 ）Ki =hf （Xn,% ）, K2 =hf （Xn +h,yn +Ki ）（i4）RungbKutta公式是在计

20、算两次函数值的情况下，局部截断误差的阶最高是 3,式（i3）是允许函数f（x, y）任意变化情况下截断误差最小的二阶方法。要再提高阶 就必须增加计算函数值的次数。上述式（i4） 乂称为欧拉预估一一校正公式。2、三级三阶 RungbKuuta方法两个常用的三阶RungbKuuta方法分别为：i yn* = yn +&h（2Ki+3K2 +4K3 ）Ki = f （Xn,% ）, i iK2 = f Xn +，h,yn +，hKi ,I 2 2 ）/工3 工3 ）（i5）K3 = f Xn + h, yn +hK2 , i 4 4i , 、yn = yn +h（Ki +4K2 + K3 ）, 6K

21、i = f Xn，yn ,Ji i ）K2 = f Xn + h,yn + hKi ,I 2 2 J（i6）K3 = f Xn h,yn - hKi 2hK2 .这两个常用方法在解决实际问题中能够达到较低的精度要求。但是要更高精度要求的，我们必须了解更高阶的方法一四级四阶 Rung* Kuuta方法。3、四级四阶RungbKutta方法这种方法在解决实际问题中常用。在这里，我们将公式进行导出，以熟悉其方法的实用过程。由丁 Runges Kutta公式对初值问题（i）中的一般f （x, y ）都适用, 则它必然对特殊的f （x, y）=Xy , f （x, y）=X或f （x, y）=y也适用。

22、从而可以定出特定的参数来，因此，这里采用的是一种待定系数法，来导出四级四阶显式Kutta公式，过程十分简单明了。为了计算简单，令=Xn h,y（xn）近似值为yn,这样求解问题（1）的四阶显示Runge Kutta公式为=忘 +顺”+小2+ +cK ）,Ki = f （Xn，yn ）,K2 = f （Xn +ha2,yn +hb2iKi）,K3 = f （Xn +ha3,yn +hb3iKi +hb43K2 ）,K4 = f （Xn +ha4,yn +hb4iKi + hb43L ）,这里要求a2 = b2i ,a b3i 炫2 , a4 = b4i b42 b43 -由丁 RungaKutt

23、a公式的思想与数值积分Xn if X,y dxXi相似。所以希望四阶显示 RungbKutta公式对y = X3, y 0 =0是精确的。这样，将（i7）式应用于（i9）式，就有Ki =Xn3,Runge（17）（18）（19）（20）K2 a?）3 = Xn3 3用端 Mx a23,_ 3K4 =（Xn a4）K3 = （Xn W3）3 = Xn3 3Xn2a3 3x2 a33, =x3 3xn2a4 3xna42 a43,yn i = yn （Cii 2-c2a 2 c ac2 c3 c4）yn h c2a2 c3& c4& y h%+c a4 2）*% 3-（c a /32C a fca ）y34h（.4）6与Taylor公式（21）相比较，有