常微分方程初值问题数值解法实用教案

1、2021-12-11 在高等数学中，对于常微分方程在高等数学中，对于常微分方程(wi fn (wi fn fn chn)fn chn)的求解，给出了一些典型方程求解的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程解法等。但能求解的常微分方程(wi fn fn (wi fn fn chn)chn)仍然是有限的，大多数的常微分方程仍然是有限的，大多数的常微分方程(wi fn fn chn)(wi fn fn chn)是不能给出解析

2、解。是不能给出解析解。 譬如譬如 22yxy 这个一阶微分方程就不能用初等这个一阶微分方程就不能用初等(chdng)(chdng)函数来函数来表达它的解。表达它的解。 第1页/共81页第一页，共82页。2021-12-12再如，方程再如，方程(fngchng) (fngchng) 1)0(yyy的解的解 , ,虽然虽然(surn)(surn)有表可查有表可查, ,但但对于表上没有给出对于表上没有给出 的值的值, ,仍需插值方法仍需插值方法来计算来计算xey xe第2页/共81页第二页，共82页。2021-12-13从实际问题当中归纳出来的微分方程从实际问题当中归纳出来的微分方程(wi fn f

3、n (wi fn fn chn)chn)，通常主要依靠数值解法来解决。本章主要讨论一，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程阶常微分方程(wi fn fn chn)(wi fn fn chn)初值问题初值问题 00)(),(yxyyxfy ( 9.1 ) 在区间在区间(q jin)a x b(q jin)a x b上的数值解法。上的数值解法。 可以证明可以证明(zhngmng),(zhngmng),如果函数在带形区域如果函数在带形区域 R: axb, R: axb,-y y内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件

4、，即存在常数条件，即存在常数 L( L(它与它与x,yx,y无关无关) )使使 2121),(),(yyLyxfyxf 对对R内任意内任意x x及两个及两个 都成立都成立, ,则方程则方程( 9.1 )的解的解 在在 a, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 第3页/共81页第三页，共82页。2021-12-14数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题(9.1)式的数值解法，就是要算式的数值解法，就是要算出精确解出精确解y(x)在区间在区间a,b上的一系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值相邻两个节点的

5、间距相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定h为定数为定数(dn sh)，称为定步长，这时节点可表示为，称为定步长，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。散节点的数值解。 bxxxxann110)(,),(),(10nxyxyxy.,10nyyyiiixxh 1niihxxi, 2 , 1,0 第4页/共81页第四页，共82页。2021-12-15 对常微分方程数值解法的基本出发点就是离对常微分方程数值解法的基本出发点就是离散化。

6、其数值解法有两个基本特点，它们散化。其数值解法有两个基本特点，它们(t men)(t men)都采用都采用“步进式步进式”，即求解过程顺着节点排列的，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求次序一步一步地向前推进，描述这类算法，要求给出用已知信息给出用已知信息 计算计算 的递推公式。建的递推公式。建立这类递推公式的基本方法是在这些节点上用数立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对值积分、数值微分、泰勒展开等离散化方法，对初值问题初值问题中的导数中的导数 进行不同的离散化处理。进行不同的离散化处理。 021,yyyyiii1i

7、yy00)(),(yxyyxfy第5页/共81页第五页，共82页。2021-12-16对于初值问题对于初值问题的数值解法，首先要解决的数值解法，首先要解决(jiju)的问题就是如何对微的问题就是如何对微分方程进行离散化，建立求数值解的递推公式。递推分方程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算公式通常有两类，一类是计算yi+1时只用到时只用到xi+1, xi 和和yi，即前一步的值，因此有了初值以后就可以逐步往，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法；其代表是龙格下计算，此类方法称为单步法；其代表是龙格库塔库塔法。另一类是计算法。另一类是

8、计算yi+1时，除用到时，除用到xi+1,xi和和yi以外，以外，还要用到还要用到，即前面，即前面k步的值，此类方法称为多步法；其代表步的值，此类方法称为多步法；其代表是亚当斯法。是亚当斯法。 00)(),(yxyyxfy), 2 , 1( ,kpyxpipi第6页/共81页第六页，共82页。2021-12-179.2 简单简单(jindn)的数值方法与基本概念的数值方法与基本概念9.2.1 Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简单）方法是解初值问题的最简单(jindn)的数值方法。初值问题的数值方法。初值问题的解的解y=y(x)代表通过点代表通过点 的一条称之为微

9、分方的一条称之为微分方程的积分曲线。积分曲线上每一点程的积分曲线。积分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在这点的值。在这点的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf第7页/共81页第七页，共82页。2021-12-18 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 Euler法的求解过程是法的求解过程是:从初始从初始点点P0(即点即点(x0,y0)出发出发,作积分曲线作积分曲线(qxin)y=y(x)在在P0点上切线点上切线 (其斜率为其斜率为 ),与与x=x1直线直线

10、10PP),()(000yxfxy相交于相交于P1P1点点( (即点即点(x1,y1),(x1,y1),得到得到y1y1作为作为y(x1)y(x1)的近似值的近似值, ,如上图所示。过点如上图所示。过点(x0,y0),(x0,y0),以以f(x0,y0)f(x0,y0)为斜率为斜率(xil)(xil)的切线方程为的切线方程为 当当 时时, ,得得 )(,(0000 xxyxfyy1xx )(,(010001xxyxfyy这样这样(zhyng)(zhyng)就获得了就获得了P1P1点的坐标。点的坐标。 第8页/共81页第八页，共82页。2021-12-19 Pi+1 Pn y=y(x) P1 P

11、i Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 同样同样, 过点过点P1(x1,y1),作积分作积分(jfn)曲线曲线y=y(x)的切线的切线交直线交直线x=x2于于P2点点,切线切线 的斜率的斜率 =直线方程为直线方程为21PP)(1xy),(11yxf)(,(1111xxyxfyy)(,(121112xxyxfyy当当 时时, ,得得 2xx 第9页/共81页第九页，共82页。2021-12-110当当 时时, ,得得 Pi+1 Pn y= y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此获得了由此获得了P2P2的坐标。

12、重复的坐标。重复(chngf)(chngf)以上过程以上过程, ,就就可获得一系列的点可获得一系列的点:P1,P1,Pn:P1,P1,Pn。对已求得点。对已求得点以以 = = 为斜率作直线为斜率作直线 ),(nnnyxP)(nxy),(nnyxf)(,(nnnnxxyxfyy1nxx)(,(11nxnnnnxxyxfyynnyxy)(取取第10页/共81页第十页，共82页。2021-12-111 从图形上看从图形上看, ,就获得了一条就获得了一条(y tio)(y tio)近似于曲近似于曲线线y=y(x)y=y(x)的折线的折线 。 Pi+ 1 Pn y= y(x) P1 Pi Pn Pi+

13、1 P0 x0 x1 xi xi+ 1 xn Pi P1 这样这样(zhyng),(zhyng),从从x0 x0逐个算出逐个算出对应的数值解对应的数值解 nxxx,21nyyy,21nPPPP321第11页/共81页第十一页，共82页。2021-12-112通常取通常取 ( (常数常数(chngsh),(chngsh),则则EulerEuler法的法的计算格式计算格式 hhxxiii1)(),(001xyyyxhfyyiiii i=0,1,n ( 9.2 ) 还可用数值微分、数值积分法和泰勒还可用数值微分、数值积分法和泰勒(ti l)(ti l)展开法推导展开法推导EulerEuler格式。以

14、数值积分为例进行推导。格式。以数值积分为例进行推导。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy 选择不同的计算选择不同的计算(j sun)(j sun)方法计算方法计算(j sun)(j sun)上式的积分项上式的积分项 , ,就会得到不同的计算就会得到不同的计算(j (j sun)sun)公式。公式。 1)(,iixxdxxyxf(9.3)第12页/共81页第十二页，共82页。2021-12-113 用左矩形方法计算用

15、左矩形方法计算(j sun)(j sun)积分积分项项 )(,)()(,11iiiixxxyxfxxdxxyxfii代入代入(9.3)(9.3)式式, ,并用并用yiyi近似代替近似代替(dit)(dit)式中式中y(xi)y(xi)即可得到向前欧拉（即可得到向前欧拉（EulerEuler）公式）公式 ),(1iiiiyxhfyy 由于由于(yuy)(yuy)数值积分的矩形方法精度很数值积分的矩形方法精度很低，所以欧拉（低，所以欧拉（EulerEuler）公式当然很粗糙。）公式当然很粗糙。 第13页/共81页第十三页，共82页。2021-12-114例例9.1 用欧拉法解初值问题用欧拉法解初值

16、问题 1)0()6 . 00(2yxxyyy取步长取步长h=0.2 ,h=0.2 ,计算过程保留计算过程保留(boli)4(boli)4位小数位小数 解解: h=0.2, : h=0.2, 欧拉迭代欧拉迭代(di di)(di di)格式格式 2),(xyyyxf21),(iiiiiiiiyhxhyyyxhfyy)2 , 1 , 0()4(2 . 0iyxyiii当当 k=0, x1=0.2时，已知时，已知x0=0，y0=1，有，有 y(0.2) y1=0.21(401)0.8当当 k=1, x2=0.4时，已知时，已知x1 =0.2， y1 =0.8，有，有 y(0.4) y2 =0.20.

17、8(40.20.8)0.6144当当 k=2, x3 =0.6时，已知时，已知x2 =0.4， y2 =0.6144，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 第14页/共81页第十四页，共82页。2021-12-1150.1,2, 01) (0)1, hEulerxyyxyy 取利用公式求解（例例 . .2 . 01.1)2( 1计算及结果如下解：nnnnnnnnyxyyxyhyyclear; y=1, x=0, %初始化for n=1:10 y=1.1*y-0.2*x/y, x=x+0.1,endy = 1 x = 0 y = 1.1000 x

18、 = 0.1000y = 1.1918 x = 0.2000y = 1.2774 x = 0.3000y = 1.3582 x = 0.4000y = 1.4351 x = 0.5000y = 1.5090 x = 0.6000y = 1.5803 x = 0.7000y = 1.6498 x = 0.8000y = 1.7178 x = 0.9000y = 1.7848 x = 1.0000第15页/共81页第十五页，共82页。2021-12-1169.2.2 梯形公式梯形公式为了提高精度为了提高精度,对方程对方程 的两端在区间上的两端在区间上 积分积分(jfn)得，得，改用梯形方法计算其积

19、分改用梯形方法计算其积分(jfn)项，即项，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii(9.4 ) 代入代入(7.4)(7.4)式式, ,并用近似代替式中即可得到并用近似代替式中即可得到(d do)(d do)梯形梯形公式公式 ),(),(2111iiiiiiyxfyxfhyy( 9.5 ) 由于由于(yuy)(yuy)数值积分的梯形公式比矩形公式的精度数值积分的梯形公式比矩形公式的精度高，因此梯形公式（高，因此梯形公式（9.59.5）比欧拉公式）比欧拉公式( 9.2 )

20、( 9.2 )的精度高的精度高一个数值方法。一个数值方法。 第16页/共81页第十六页，共82页。2021-12-117),(),(2111iiiiiiyxfyxfhyy( 9.5 ) (9.5) (9.5)式的右端含有未知的式的右端含有未知的yi+1,yi+1,它是一个关于它是一个关于yi+1yi+1的函数的函数(hnsh)(hnsh)方程方程, ,这类数值方法称为隐式方这类数值方法称为隐式方法。相反地法。相反地, ,欧拉法是关于欧拉法是关于yi+1yi+1的一个直接的计算公的一个直接的计算公式，式， 这类数值方法称为显式方法。这类数值方法称为显式方法。 第17页/共81页第十七页，共82页

21、。2021-12-1189.2.3 两步欧拉公式两步欧拉公式 对方程对方程 的两端在区间的两端在区间(q jin)上上 积分得积分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 9.6 ) 改用中矩形公式改用中矩形公式(gngsh)(gngsh)计算其积分项，即计算其积分项，即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式, ,并用并用yiyi近似代替近似代替(dit)(dit)式中式中y(xi)y(xi)即可得到即可得到两步欧拉公式两步欧拉公式 ),(211iiiiyxhfyy ( 9.7 ) 第18页/共81页第

22、十八页，共82页。2021-12-119 前面介绍过的数值方法前面介绍过的数值方法, ,无论是欧拉方法无论是欧拉方法, ,还是梯形方法，它们都是单步法还是梯形方法，它们都是单步法, ,其特点是在其特点是在计算计算yi+1yi+1时只用到前一步时只用到前一步(y b)(y b)的信息的信息yi;yi;可可是公式是公式(7.7)(7.7)中除了中除了yiyi外外, ,还用到更前一步还用到更前一步(y (y b)b)的信息的信息yi-1,yi-1,即调用了前两步的信息即调用了前两步的信息, ,故称故称其为两步欧拉公式其为两步欧拉公式 第19页/共81页第十九页，共82页。2021-12-1209.2

23、.4 欧拉法的局部截断误差欧拉法的局部截断误差(wch) 衡量求解公式好坏的一个主要标准是求解公式的精度衡量求解公式好坏的一个主要标准是求解公式的精度, 因此引入局部截断误差因此引入局部截断误差(wch)和阶数的概念。和阶数的概念。 定义定义9.1 在在yi准确的前提下准确的前提下, 即即 时时, 用数值方法用数值方法计算计算yi+1的误差的误差(wch) , 称为该数值方法称为该数值方法计算时计算时yi+1的局部截断误差的局部截断误差(wch)。 对于欧拉公式，假定对于欧拉公式，假定 ，则有，则有)(iixyy 11)(iiiyxyR)(iixyy )()()(,()(1iiiiiixyhx

24、yxyxfhxyy而将真解而将真解y(x)在在xi处按二阶泰勒处按二阶泰勒(ti l)展开展开 ),()(! 2)()()(121 iiiiixxyhxyhxyxy)(!2)(211yhyxyii 因此因此(ync)有有 第20页/共81页第二十页，共82页。2021-12-121定义定义9.2 数值方法的局部截断误差为数值方法的局部截断误差为 ,则称这种数则称这种数值方法的阶数是值方法的阶数是P。步长。步长(h N 结束。结束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxx第25页/共81页第二十五页，共82页。2021

25、-12-126（2 2）改进）改进(gijn)(gijn)欧拉法的流程欧拉法的流程图图 开 始 输 入x0, y0,h , N 1 n x0 + h x1 y0+ h f( x0,y0 ) yp y0+ h f( x1,yp) yc ( yp+ yc) /2 y1 输 出x1, y1 n + 1 n n = N ? x1 x0 y1 y0 结 束 n y 第26页/共81页第二十六页，共82页。2021-12-127(3) (3) 程序实现程序实现( (改进欧拉法计算改进欧拉法计算(j sun)(j sun)常微常微 分方程初值问题分方程初值问题 ) ) 例例9.2 9.2 用改进用改进(gi

26、jn)(gijn)欧拉法解初值问题欧拉法解初值问题 1)0(2yyxyy区间区间(q jin)(q jin)为为0,10,1, ,取步长取步长h=0.1 h=0.1 解解: : 改进欧拉法的具体形式改进欧拉法的具体形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyyxyyyyxyyy本题的精确解为本题的精确解为 , ,xxy21)(第27页/共81页第二十七页，共82页。2021-12-128clearx=0,yn=1 %初始化for n=1:10yp=yn+0.1*(yn-2*x/yn); %预测(yc)x=x+0.1;yc=yn+0.1*(yp-2*x/yp) ;y

27、n=(yp+yc)/2 %校正end第28页/共81页第二十八页，共82页。2021-12-129例例9.3 对初值问题对初值问题 1)0(0yyy证明证明(zhngmng)(zhngmng)用梯形公式求得的用梯形公式求得的近似解为近似解为 nnhhynhx 并证明并证明(zhngmng)(zhngmng)当步长当步长h h0 0时时,yn,yn收敛收敛于精确解于精确解证明证明(zhngmng): (zhngmng): 解初值问题的梯形公解初值问题的梯形公式为式为 xe),(),(nnnnnnyxfyxfhyyyyxf),( 211nnnnyyhyy 整理整理(zhngl)(zhngl)成显式

28、成显式 nnyhhy反复迭代反复迭代, ,得到得到 yhhyhhyhhyhhynnnnn.10ynnhhy22 第29页/共81页第二十九页，共82页。2021-12-130由于由于(yuy) (yuy) ，有，有 nhx xxxxhxhhhxhnhhhhhyeee2121lim22limlim222222000nnhhy22xnhy elim0 证毕证毕 第30页/共81页第三十页，共82页。2021-12-1319.3 9.3 龙格龙格- -库塔（库塔（Runge-KuttaRunge-Kutta）法）法9.3.1 9.3.1 龙格龙格- -库塔库塔(Runge-Kutta)(Runge-

29、Kutta)法的基本法的基本(jbn)(jbn)思想思想 Euler Euler公式可改写成公式可改写成 ),(111iiiiyxfKhKyy则则yi+1yi+1的表达式的表达式y(xi+1)y(xi+1)与的与的TaylorTaylor展开式的前两项展开式的前两项完全相同完全相同, ,即局部即局部(jb)(jb)截断误差为截断误差为 。 改进的改进的EulerEuler公式又可改写成公式又可改写成 )(2hO),(),()(21121211hKyxfKyxfKKKhyyiiiiii第31页/共81页第三十一页，共82页。2021-12-132 上述两组公式在形式上有一个共同点上述两组公式在形

30、式上有一个共同点: :都是用都是用f(x,y)f(x,y)在某些在某些(mu xi)(mu xi)点上值的线性组合得出点上值的线性组合得出y(xi+1)y(xi+1)的近似值的近似值yi+1,yi+1,而且增加计算的次数而且增加计算的次数f(x,y)f(x,y)的次数的次数, ,可提高截断误差的阶。如欧拉公式可提高截断误差的阶。如欧拉公式: :每步计每步计算一次算一次f(x,y)f(x,y)的值的值, ,为一阶方法。改进欧拉公式需为一阶方法。改进欧拉公式需计算两次计算两次f(x,y)f(x,y)的值，它是二阶方法。它的局部截的值，它是二阶方法。它的局部截断误差为断误差为 。)(3hO第32页/

31、共81页第三十二页，共82页。2021-12-133 于是可考虑用函数于是可考虑用函数f(x,y)f(x,y)在若干点上的函在若干点上的函数值的线性组合来构造近似公式，构造时要求数值的线性组合来构造近似公式，构造时要求近似公式在近似公式在(xi,yi)(xi,yi)处的处的TaylorTaylor展开式与解展开式与解y(x)y(x)在在xixi处的处的TaylorTaylor展开式的前面几项重合，从而展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免使近似公式达到所需要的阶数。既避免(bmin)(bmin)求偏导求偏导, ,又提高了计算方法精度的阶数。或者说又提高了计算方法精度的阶数

32、。或者说, ,在在 这一步内多预报几个点的斜率值，然后将这一步内多预报几个点的斜率值，然后将其加权平均作为平均斜率，则可构造出更高精其加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格度的计算格式，这就是龙格库塔（库塔（Runge-Runge-KuttaKutta）法的基本思想。）法的基本思想。 1,iixx第33页/共81页第三十三页，共82页。2021-12-1349.3.2 二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi和和 ,以该两点处的斜率值以该两点处的斜率值k1和和k2的加权平均的加权平均(pngjn)(或称为线性组合或称为线性组合)来求取来求取平均平均(pn

33、gjn)斜率斜率k*的近似值的近似值K，即，即 1,iixxphxxipi2211kkK式中式中:k1:k1为为xixi点处的切线斜率点处的切线斜率(xil)(xil)值，值， k2 k2为为 点处的切线斜率点处的切线斜率(xil)(xil)值值, ,比照改比照改进的欧拉法进的欧拉法, ,将将 视为视为 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12phkyphxfkii对常微分方程初值问题对常微分方程初值问题(9.1)(9.1)式的解式的解 y=y(x), y=y(x),根据微分中根据微分中值定理值定理(dngl)(dngl)，存在点，存在点 ，使得，使得 )

34、,(1iixx第34页/共81页第三十四页，共82页。2021-12-135)(,()(yfyK式中式中 K K可看作是可看作是y=y(x)y=y(x)在区间在区间 上的平均斜率上的平均斜率(xil)(xil)。所以可得计算公式为：。所以可得计算公式为： 1,iixxhKxyxyii)()(1)()(2211kkhxyi（9.14） 将将y(xi)y(xi)在在x=xix=xi处进行处进行(jnxng)(jnxng)二阶二阶TaylorTaylor展展开：开： )()(! 2)()()(321hOxyhxyhxyxyiiii （9 9.15） )()()(11iiiixxyxyxy也即也即 h

35、Kxyxyii)()(1（9.13）第35页/共81页第三十五页，共82页。2021-12-136将将 在在x=xix=xi处进行处进行(jnxng)(jnxng)一阶一阶TaylorTaylor展开：展开： ),()(12phkyphxfphxykiii)(),(),(),(),(22hOyxfyxfyxfphyxfkiiyiiiixii)()()(2hOxyphxyii 将以上将以上(yshng)(yshng)结果代入（结果代入（9.149.14）得：）得： )()()(22111kkhxyxyii)()()()()(221hOxyphxyxyhxyiiii )()()()()(32221

36、hOxyphxyhxyiii （9.16） 对式对式(9.15)(9.15)和和(9.16)(9.16)进行比较系数后可知进行比较系数后可知(k zh),(k zh),只要只要 211221p（9.17） 成立成立, ,格式格式(9.14)(9.14)的局部截断误差就等于的局部截断误差就等于)(3hO有有2 2阶阶精度精度第36页/共81页第三十六页，共82页。2021-12-137式式(9.17)(9.17)中具有三个未知量中具有三个未知量, ,但只有两个方程但只有两个方程(fngchng),(fngchng),因而有无穷多解。若取因而有无穷多解。若取 , ,则则p=1p=1，这是无穷多解中

37、的一个解，将以上所解的值代入式这是无穷多解中的一个解，将以上所解的值代入式(9.14)(9.14)并改写可得并改写可得 2121),(),()(21121211hkyxfkyxfkkkhyyiiiiii 不难发现，上面的格式就是改进不难发现，上面的格式就是改进(gijn)(gijn)的欧拉格的欧拉格式。凡满足条件式（式。凡满足条件式（9.179.17）有一簇形如上式的计算格式，）有一簇形如上式的计算格式，这些格式统称为二阶龙格这些格式统称为二阶龙格库塔格式。因此改进库塔格式。因此改进(gijn)(gijn)的欧拉格式是众多的二阶龙格的欧拉格式是众多的二阶龙格库塔法中的一种特殊格库塔法中的一种特

38、殊格式。式。 第37页/共81页第三十七页，共82页。2021-12-138若取若取 , ,则则 ，此时，此时(c sh)(c sh)二阶龙格二阶龙格- -库库塔塔法的计算公式为法的计算公式为 0121, 12p)2,(),(1212121khyxfkyxfkhkyyiiiiii1,2 , 1 , 0ni 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间(q jin) (q jin) 的中点。的中点。 21ix1,iixx第38页/共81页第三十八页，共82页。2021-12-1399.3.3 三阶三阶(sn ji)龙格龙格-库塔法库塔法 为了进一

39、步提高精度，设除为了进一步提高精度，设除 外再增加外再增加(zngji)(zngji)一一点点 pix) 1(qpqhxxiqi并用三个点并用三个点 , , , , 的斜率的斜率(xil)k1,k2,k3(xil)k1,k2,k3加权平均加权平均得出平均斜率得出平均斜率(xil)k(xil)k* *的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式: : ixpixqix11 12233121(,)(,)iiiiiiyyhkkkkf x ykf xph yphk（9.18） 为了预报点为了预报点 的斜率值的斜率值k3 3, ,在区间在区间 内有两内有两个斜率值个斜率值k1 1和和k2

40、 2可以用可以用, ,可将可将k1 1, ,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率, ,从而得到从而得到 的预报值的预报值 qixqiixx,qiixx,)(qixyqiy1122i qiyyqhkk第39页/共81页第三十九页，共82页。2021-12-140于是于是(ysh)(ysh)可得可得 ),(3qiqiyxfk运用运用TaylorTaylor展开展开(zhn ki)(zhn ki)方法选择参数方法选择参数 , ,可可以使格式以使格式(9.18)(9.18)的局部截断误差为的局部截断误差为 , ,即具有三阶精度，即具有三阶精度，这类格式统称为三阶龙格这类格式统称为

41、三阶龙格库塔方法。下列是其中的一种，库塔方法。下列是其中的一种，称为库塔（称为库塔（KuttaKutta）公式。）公式。 12212, ,p q )(4hO)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（9.19） 第40页/共81页第四十页，共82页。2021-12-141121231232221233132611pqpqpq 1141123122666,1,1,2pq 第41页/共81页第四十一页，共82页。2021-12-1429.3.4 四阶龙格四阶龙格库塔法库塔法 如果需要再提高精度，用类似上述的处理方法，如果

42、需要再提高精度，用类似上述的处理方法，只需在区间只需在区间 上用四个点处的斜率加权平均作为上用四个点处的斜率加权平均作为平均斜率平均斜率k*的近似值，构成一系列四阶龙格的近似值，构成一系列四阶龙格库塔公库塔公式。具有式。具有(jyu)四阶精度，即局部截断误差是四阶精度，即局部截断误差是 。 由于推导复杂，这里从略，只介绍最常用的一种四由于推导复杂，这里从略，只介绍最常用的一种四阶经典龙格阶经典龙格库塔公式。库塔公式。 qiixx,)(5hO)22(6),()2,()2,(),(43211314221312121kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（9.

43、20） 第42页/共81页第四十二页，共82页。2021-12-1439.3.5 9.3.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1) 计算步骤计算步骤 输入输入 ,h,N ,h,N 使用龙格使用龙格库塔公式库塔公式(gngsh)(gngsh)（9.209.20）计算）计算出出y1y1 输出输出 ，并使，并使 转到转到 直至直至n N n N 结束。结束。 10,xx11, yx0101,yyxx第43页/共81页第四十三页，共82页。2021-12-144(2) (2) 四阶龙格四阶龙格库塔算法库塔算法(sun f)(sun f)流程图流程图 开 始 输 入 x0, y0,h

44、 , N 1 n x0 + h x1 f(x0,y0 ) k1, f(x0+ h /2 ,y0 + h k1/2 ) k2 f(x0+ h /2 ,y0 + h k2/2 ) k3, f(x1,y0+ h k3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 输 出 x1, y1 n + 1 n n = N ? x1 x0 y1 y0 结 束 n y 第44页/共81页第四十四页，共82页。2021-12-145(3)(3) 程序实现程序实现( (四阶龙格四阶龙格- -库塔法计库塔法计(4)(4) 算常微分方程算常微分方程(wi fn fn chn)(wi fn f

45、n chn)初值问初值问题题) ) 例例9.4 9.4 取步长取步长h=0.2h=0.2，用经典格式，用经典格式(g shi)(g shi)求解初求解初值问题值问题 1)0(2yxyy10 x解解: : 由四阶龙格由四阶龙格- -库塔公式库塔公式(gngsh)(gngsh)可得可得 2 . 0, 1, 0,2),(00hyxxyyxf0),(001yxfk2 . 0) 1 , 1 . 0()2,(10202fkhyxfkh204. 0)02. 1 , 1 . 0()2,(20203fkhyxfkh41632. 0)0408. 1 , 2 . 0(),(3004fhkyxfkh第45页/共81页

46、第四十五页，共82页。2021-12-146040811. 1)22(62 . 0432101kkkkyy可同样进行其余可同样进行其余yiyi的计算。本例方程的解为的计算。本例方程的解为，从表中看到所求的数值，从表中看到所求的数值(shz)(shz)解具有解具有4 4位有效数字。位有效数字。 2xey 龙格龙格库塔方法的推导基于库塔方法的推导基于TaylorTaylor展开展开(zhn ki)(zhn ki)方法，因而它要求所求的解具有方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，较好的光滑性。如果解的光滑性差，那么，使用四阶龙格使用四阶龙格库塔方法求得的数值解，其库塔方法

47、求得的数值解，其精度可能反而不如改进的欧拉方法。在实际精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适计算时，应当针对问题的具体特点选择合适的算法。的算法。 第46页/共81页第四十六页，共82页。2021-12-1479.3.6 变步长的龙格变步长的龙格-库塔法库塔法 在微分方程的数值解中，选择适当的步长是在微分方程的数值解中，选择适当的步长是非常重要的。单从每一步看，步长越小，截断非常重要的。单从每一步看，步长越小，截断误差就越小；但随着步长的缩小，在一定的求误差就越小；但随着步长的缩小，在一定的求解区间内所要完成解区间内所要完成(wn chng)的步数就增加的

48、步数就增加了。这样会引起计算量的增大，并且会引起舍了。这样会引起计算量的增大，并且会引起舍入误差的大量积累与传播。因此微分方程数值入误差的大量积累与传播。因此微分方程数值解法也有选择步长的问题。解法也有选择步长的问题。 以经典的四阶龙格以经典的四阶龙格-库塔法库塔法(9.20)为例。从节为例。从节点点xi出发，先以出发，先以h为步长求出一个近似值，记为步长求出一个近似值，记为为 ，由于局部截断误差为，由于局部截断误差为 ，故有，故有 )(1hiy)(5hO5)(11)(chyxyhii当h值不大时，式中的系数c可近似(jn s)地看作为常数。第47页/共81页第四十七页，共82页。2021-1

49、2-148然后将步长折半然后将步长折半, ,即以为即以为 步长步长, ,从节点从节点(ji (ji din)xidin)xi出发出发, ,跨两步到节点跨两步到节点(ji din)xi+1,(ji din)xi+1,再求得再求得一个近似值一个近似值 , ,每跨一步的截断误差是每跨一步的截断误差是 , ,因此有因此有2h)2(1hiy52 hc5)2(1122)()(hcxyxyhii这样这样(zh(zhyng) yng) 161)()()(11)2(11hiihiiyxyyxy)(151)()(1)2(1)2(11hihihiiyyyxy由此可得由此可得 这表明以这表明以 作为作为 的近似值，其

50、误差可用步的近似值，其误差可用步长折半长折半(zhbn)(zhbn)前后两次计算结果的偏差前后两次计算结果的偏差 )2(1hiy)(1ixy)(1)2(1hihiyy来判断所选步长是否适当来判断所选步长是否适当第48页/共81页第四十八页，共82页。2021-12-149当要求的数值当要求的数值(shz)(shz)精度为精度为时：时： （1 1）如果）如果，反复将步长折半进行计算，直至，反复将步长折半进行计算，直至为止为止, ,并取其最后一次步长的计算结果作为并取其最后一次步长的计算结果作为 （2 2）如果）如果为止，为止，并以上并以上(yshng)(yshng)一次步长的计算结果作为一次步长

51、的计算结果作为 。 这种通过步长加倍或折半来处理步长的方法称为变步这种通过步长加倍或折半来处理步长的方法称为变步长法。表面上看，为了选择步长，每一步都要反复判断长法。表面上看，为了选择步长，每一步都要反复判断，增加了计算工作量，但在方程的解增加了计算工作量，但在方程的解y(x)y(x)变化剧烈的情况下，变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的。总的计算工作量得到减少，结果还是合算的。1iy1iy第49页/共81页第四十九页，共82页。2021-12-150 9.4 算法的稳定性及收敛性算法的稳定性及收敛性9.4.1 稳定性稳定性 稳定性在微分方程的数值解法中是一个非常重要的稳定

52、性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分格式问题。因为微分方程初值问题的数值方法是用差分格式(g shi)进行计算的，而在差分方程的求解过程中，存在进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误差等引起的扰着各种计算误差，这些计算误差如舍入误差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。确性将产生影响。这就涉及到算法稳定性问题。 第50页/共81页第五十页，共82页。2021-12-151 当在某节点上当在某节点上x

53、i的的yi值有大小为值有大小为的扰动时，如果在的扰动时，如果在其后的各节点其后的各节点 上的值上的值yi产生的偏差都不大于产生的偏差都不大于，则，则称这种方法是稳定的。称这种方法是稳定的。 稳定性不仅稳定性不仅(bjn)与算法有关，而且与方程中函数与算法有关，而且与方程中函数f(x,y)也有关，讨论起来比较复杂。为简单起见，通常也有关，讨论起来比较复杂。为简单起见，通常只针对模型方程只针对模型方程)(ijxj)0(yy来讨论来讨论(toln)(toln)。一般方程若局部线性化，也可化为。一般方程若局部线性化，也可化为上述形式。模型方程相对比较简单，若一个数值方法上述形式。模型方程相对比较简单，

54、若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程对模型方程是稳定的，并不能保证该方法对任何方程都稳定，但若某方法对模型方程都不稳定，也就很难都稳定，但若某方法对模型方程都不稳定，也就很难用于其他方程的求解。用于其他方程的求解。第51页/共81页第五十一页，共82页。2021-12-152先考察显式先考察显式EulerEuler方法的稳定性。模型方法的稳定性。模型(mxng)(mxng)方程方程的的EulerEuler公式为公式为 yy)(),(1iiiiiiyhyyxhfyy), 2 , 1 , 0()1 (iyhi将上式反复将上式反复(fnf)(fnf)递推后，可得递推后，可得

55、011)1(yhyii), 2 , 1()1 (00iyyhyiii或或h1式中式中第52页/共81页第五十二页，共82页。2021-12-153 要使要使y yi i有界，其充要条件为有界，其充要条件为 111h即即 由于由于(yuy) (yuy) ，故，故有有 020 h（9.269.26） 可见，如欲保证算法的稳定，显式可见，如欲保证算法的稳定，显式EulerEuler格式格式(g (g shi)shi)的步长的步长h h的选取要受到式（的选取要受到式（9.269.26）的限制。）的限制。 的的绝对值越大，则限制的绝对值越大，则限制的h h值就越小。值就越小。 用隐式用隐式EulerEu

56、ler格式格式(g shi)(g shi)，对模型方程，对模型方程 的计算公式为，可化为的计算公式为，可化为)(11iiiyhyyiiyhy111第53页/共81页第五十三页，共82页。2021-12-154由于由于(yuy) ,(yuy) ,则恒有则恒有 , ,故恒有故恒有 0111hiiyy1 因此，隐式因此，隐式EulerEuler格式是绝对稳定格式是绝对稳定(wndng)(wndng)的（无的（无条件稳定条件稳定(wndng)(wndng)）的（对任何）的（对任何h0h0）。）。 9.4.2 9.4.2 收敛收敛(shulin)(shulin)性性 常微分方程初值问题的求解，是将微分方

57、程转化常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值为差分方程来求解，并用计算值yiyi来近似替代来近似替代y(xi)y(xi)，这种近似替代是否合理，还须看分割区间这种近似替代是否合理，还须看分割区间 的的长度长度h h越来越小时，即越来越小时，即 时，时， 是是否成立。若成立，则称该方法是收敛否成立。若成立，则称该方法是收敛(shulin)(shulin)的，的，否则称为不收敛否则称为不收敛(shulin)(shulin)。 这里仍以这里仍以EulerEuler方法为例，来分析其收敛方法为例，来分析其收敛(shulin)(shulin)性。性。EulerEuler格

58、式格式iixx,101iixxh)(iixyy),(1iiiiyxhfyy第54页/共81页第五十四页，共82页。2021-12-155设设 表示表示(biosh)取取 时时, 按按Euler公式的计算公式的计算结果结果, 即即1iy)(iixyy )(,()(1iiiixyxhfxyyEuler方法方法(fngf)局部截断误差为：局部截断误差为：)()(2)(1211 iiiixxyhyxy设有常数设有常数(chngsh) (chngsh) ，则，则 )(max21xycbxa 211)(chyxyii（9.27） 总体截断误差总体截断误差 1111111)()(iiiiiiiyyyxyyx

59、y)(,()(),(11iiiiiiiixyxhfxyyxhfyyy),()(,()(iiiiiiyxfxyxfhyxy（9.28） 又又 由于由于f(x,y)f(x,y)关于关于y y满足李普希茨条件满足李普希茨条件, ,即即 第55页/共81页第五十五页，共82页。2021-12-156iiiiiiyxyLyxfxyxf)(),()(代入代入(9.28)上式，有上式，有 iiiiyxyhLyy)()1(11ihL)1(再利用再利用(lyng)式（式（9.27），式（），式（9.28） 1111111)()(iiiiiiiyyyxyyxy2)1 (chhLi21)1 (chhLii即即 上式

60、反复上式反复(fnf)递推后，递推后，可得可得 1020)1 ()1 (ikkiihLchhL1)1 ()1 (0iihLLchhL（9.29） 第56页/共81页第五十六页，共82页。2021-12-157设设 （T T为常数为常数(chngsh)(chngsh)） Tihxxi0hLehL 1TLihLieehL)1 (因为因为(yn (yn wi) wi) 所以所以(su(suy) y) 把上式代入式（把上式代入式（9.299.29），得），得 ) 1(0TLTLieLche若不计初值误差，即若不计初值误差，即 ，则有，则有 00heLcTLi)1(（9.30） 式式(9.30)(9.3