常微分方程初值问题数值解法.ppt

《常微分方程初值问题数值解法.ppt》由会员分享，可在线阅读，更多相关《常微分方程初值问题数值解法.ppt（81页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、常微分方程初值问题数值解法现在学习的是第1页，共81页 在高等数学中，对于常微分方程的求解，给出了在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。程的解法等。但能求解的常微分方程仍然是有限的，但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。大多数的常微分方程是不可能给出解析解。譬如譬如 22yxy 这个一阶微分方程就不能用初等函数及其积分来表这个一阶微分方程就不能用初等函数及其积分来

2、表达它的解。达它的解。现在学习的是第2页，共81页再如，方程再如，方程 1)0(yyy的解的解 ,虽然有表可查虽然有表可查,但对于表上但对于表上没有给出没有给出 的值的值,仍需插值方法来计算仍需插值方法来计算xey xe现在学习的是第3页，共81页从实际问题当中归纳出来的微分方程，通常主要依靠从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值数值解法来解决。本章主要讨论一阶常微分方程初值问题问题 00)(),(yxyyxfy(9.1)在区间在区间a x ba x b上的数值解法上的数值解法。可以证明可以证明,如果函数在带形区域如果函数在带形区域 R=a

3、xb,R=axb,-y y内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，即存在常数条件，即存在常数L(L(它与它与x,yx,y无关无关)使使 2121),(),(yyLyxfyxf对对R内任意两个内任意两个 都成立都成立,则方程则方程(9.1)的解的解 在在 a,b 上存在且唯一。上存在且唯一。21,yy)(xyy 现在学习的是第4页，共81页数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题(9.1)式的数值解法，就是式的数值解法，就是要算出精确解要算出精确解y(x)y(x)在区间在区间 a,b 上的一

4、系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值。相邻两个节点的间距。相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定h为定数，称为定数，称为为定步长定步长，这时节点可表示为，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。离散节点的数值解。bxxxxann110)(,),(),(10nxyxyxynyyy,10iixxh1niihxxi,2,1 ,0现在学习的是第5页，共81页 对常微分方程数值解法的对常微分方程数值解法的基本

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

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

7、用到x xi+1i+1,x,xi i和和y yi i以外，还要用到以外，还要用到，即前面，即前面k步的值，此类方法称为步的值，此类方法称为多步法多步法；其代表；其代表是亚当斯法。是亚当斯法。00)(),(yxyyxfy),2,1(,kpyxpipi现在学习的是第7页，共81页9.2 简单的数值方法与基本概念简单的数值方法与基本概念9.2.1 Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简单）方法是解初值问题的最简单的数值方法。初值问题的数值方法。初值问题的解的解y=y(x)y=y(x)代表通过点代表通过点 的一条称之为微分的一条称之为微分方程的积分曲线。积分曲线上每一点方程

8、的积分曲线。积分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在这在这点的值。点的值。00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf现在学习的是第8页，共81页 Pi+1 Pn y=y(x)P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 Euler法的求解过程是法的求解过程是:从初始点从初始点P0(即点即点(x(x0 0,y,y0 0)出发出发,作积分曲线作积分曲线y=y(x)y=y(x)在在P0点上切点上切线线 (其斜率为其斜率为 ),),与与x=xx=x1 1直线直线10PP),()(000yxfxy相交于相交于P1点点

9、(即点即点(x(x1 1,y,y1 1),),得到得到y y1 1作为作为y(xy(x1 1)的近似值的近似值,如上如上图所示。过点图所示。过点(x(x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0)为为斜率的切线方程为斜率的切线方程为 当当 时时,得得 )(,(0000 xxyxfyy1xx)(,(010001xxyxfyy这样就获得了这样就获得了P P1 1点的坐标。点的坐标。现在学习的是第9页，共81页 Pi+1 Pn y=y(x)P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 同样同样,过过点点P1(x x1 1,y,y1 1),),作

10、积分曲线作积分曲线y=y(x)y=y(x)的切线的切线交直线交直线x=xx=x2 2于于P2点点,切线切线 的斜率的斜率 =直线方程为直线方程为21PP)(1xy),(11yxf)(,(1111xxyxfyy)(,(121112xxyxfyy当当 时时,得得 2xx 现在学习的是第10页，共81页当当 时时,得得 Pi+1 Pn y=y(x)P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程,就可获得一系列的点就可获得一系列的点:P P1 1,P P1 1,P Pn n。对已求得点对已求得点以

11、以 =为斜率作直线为斜率作直线 ),(nnnyxP)(nxy),(nnyxf)(,(nnnnxxyxfyy1nxx)(,(11nxnnnnxxyxfyynnyxy)(取取现在学习的是第11页，共81页 从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)y=y(x)的折线的折线 。Pi+1 Pn y=y(x)P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 这样这样,从从x x0 0逐个算出逐个算出对应的数值解对应的数值解 nxxx,21nyyy,21nPPPP321现在学习的是第12页，共81页通常取通常取 (常数常数),),则则Eul

12、er法的计算格式法的计算格式 hhxxiii1)(),(001xyyyxhfyyiiii i=0,1,n (9.2)还可用还可用数值微分数值微分、数值积分法数值积分法和和泰勒展开法泰勒展开法推导推导EulerEuler格格式。以数值积分为例进行推导。式。以数值积分为例进行推导。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得，),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy 选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项 ,就会得到不同的计算公式。就会得到

13、不同的计算公式。1)(,iixxdxxyxf(9.3)现在学习的是第13页，共81页 用左矩形方法计算积分项用左矩形方法计算积分项 )(,)()(,11iiiixxxyxfxxdxxyxfii代入代入(9.3)(9.3)式式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到向前欧即可得到向前欧拉（拉（EulerEuler）公式）公式 ),(1iiiiyxhfyy 由于数值积分的矩形方法精度很低，所以欧拉（由于数值积分的矩形方法精度很低，所以欧拉（EulerEuler）公式当然很粗糙。）公式当然很粗糙。现在学习的是第14页，共81页例例9.1 用欧拉法解初值问题用欧拉法解

14、初值问题 1)0()6.00(2yxxyyy取步长取步长h=0.2,h=0.2,计算过程保留计算过程保留4 4位小数位小数 解解:h=0.2,:h=0.2,欧拉迭代格式欧拉迭代格式 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.8(40.20.8)0.6144当当 k=2,x3=0.6时，已知时，已知x2=0.4，y

15、2=0.6144，有，有 y(0.6)y3=0.20.6144(4-0.40.6144)=0.4613 现在学习的是第15页，共81页0.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=0.1000y=1.1918 x=0.2000y=1.2774 x=0.3000y=1.3582 x=0.4000y=1.4351 x=0.5000y=1.5090

16、x=0.6000y=1.5803 x=0.7000y=1.6498 x=0.8000y=1.7178 x=0.9000y=1.7848 x=1.0000现在学习的是第16页，共81页9.2.2 梯形公式梯形公式为了提高精度为了提高精度,对方程对方程 的两端在区间上的两端在区间上 积分得，积分得，改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii(9.4)代入代入(7.4)(7.4)式式,并用近似代替式中即可得到梯形公式并用近似

17、代替式中即可得到梯形公式 ),(),(2111iiiiiiyxfyxfhyy(9.5)由于数值积分的梯形公式比矩形公式的精度高，因由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（此梯形公式（9.59.5）比欧拉公式）比欧拉公式(9.2)(9.2)的精度高一个数的精度高一个数值方法。值方法。现在学习的是第17页，共81页),(),(2111iiiiiiyxfyxfhyy(9.5)(9.5)式的右端含有未知的式的右端含有未知的y yi+1i+1,它是一个关于它是一个关于y yi+1i+1的函的函数方程数方程,这类数值方法称为这类数值方法称为隐式方法隐式方法。相反地。相反地,欧拉法是关欧拉法

18、是关于于y yi+1i+1的一个直接的计算公式，的一个直接的计算公式，这类数值方法称为这类数值方法称为显式方显式方法。法。现在学习的是第18页，共81页9.2.3 两步欧拉公式两步欧拉公式 对方程对方程 的两端在区间上的两端在区间上 积分得积分得),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy(9.6)改用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到两步欧拉即可得到两步欧拉公式公式 ),(211i

19、iiiyxhfyy(9.7)现在学习的是第19页，共81页 前面介绍过的数值方法前面介绍过的数值方法,无论是欧拉方法无论是欧拉方法,还是梯形方法，它们都是单步法还是梯形方法，它们都是单步法,其特点是在其特点是在计算计算y yi+1i+1时只用到前一步的信息时只用到前一步的信息y yi i;可是公式可是公式(7.7)中除了中除了y yi i外外,还用到更前一步的信息还用到更前一步的信息y yi-1i-1,即调用了前两步的信息即调用了前两步的信息,故称其为两步欧拉公故称其为两步欧拉公式式 现在学习的是第20页，共81页9.2.4 欧拉法的局部截断误差欧拉法的局部截断误差 衡量求解公式好坏的一个主要

20、标准是求解公式的精度衡量求解公式好坏的一个主要标准是求解公式的精度,因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。定义定义9.1 在在yi准确准确的前提下的前提下,即即 时时,用数值方法计算用数值方法计算yi+1的误差的误差 ,称为该数值方法计算时称为该数值方法计算时yi+1的的局部截断误差。局部截断误差。对于欧拉公式，假定对于欧拉公式，假定 ，则有，则有)(iixyy 11)(iiiyxyR)(iixyy)()()(,()(1iiiiiixyhxyxyxfhxyy而将真解而将真解y(x)在在xi处按二阶泰勒展开处按二阶泰勒展开),()(!2)()()(121 iiiii

21、xxyhxyhxyxy)(!2)(211yhyxyii 因此有因此有 现在学习的是第21页，共81页定义定义9.2 数值方法的局部截断误差为数值方法的局部截断误差为 ,则称这种数值则称这种数值方法的阶数是方法的阶数是P。步长步长(h N 结束。结束。10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11,yx0101,yyxx现在学习的是第26页，共81页（2）改进欧拉法的流程图）改进欧拉法的流程图 开 始 输 入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

22、,y1 n+1 n n=N?x1 x0 y1 y0 结 束 n y 现在学习的是第27页，共81页(3)3)程序实现程序实现(改进欧拉法计算常微改进欧拉法计算常微 分方程初值问题分方程初值问题)例例9.2 9.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 1)0(2yyxyy区间为区间为 0,10,1,取步长取步长h=0.1h=0.1 解解:改进欧拉法的具体形式改进欧拉法的具体形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyyxyyyyxyyy本题的精确解为本题的精确解为 ,xxy21)(现在学习的是第28页，共81页clearx=0,yn=1%初始化for n=

23、1:10yp=yn+0.1*(yn-2*x/yn);%预测x=x+0.1;yc=yn+0.1*(yp-2*x/yp);yn=(yp+yc)/2%校正end现在学习的是第29页，共81页例例9.3 对初值问题对初值问题 1)0(0yyy证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为 nnhhynhx 并证明当步长并证明当步长h h0 0时时,y,yn n收敛于精确解收敛于精确解证明证明:解初值问题的梯形公式为解初值问题的梯形公式为 xe),(),(nnnnnnyxfyxfhyyyyxf),(211nnnnyyhyy 整理成显式整理成显式 nnyhhy反复迭代反复迭代,得到得到 yhhy

24、hhyhhyhhynnnnn.10ynnhhy22 现在学习的是第30页，共81页由于由于 ，有，有 nhx xxxxhxhhhxhnhhhhhyeee2121lim22limlim222222000nnhhy22xnhy elim0 证毕证毕 现在学习的是第31页，共81页9.3 9.3 龙格龙格-库塔（库塔（Runge-KuttaRunge-Kutta）法）法9.3.1 9.3.1 龙格龙格-库塔库塔(Runge-Kutta)(Runge-Kutta)法的基本思想法的基本思想 Euler Euler公式可改写成公式可改写成 ),(111iiiiyxfKhKyy则则yi+1i+1的表达式的表

25、达式y(xi+1i+1)与的与的TaylorTaylor展开式的前两项完全展开式的前两项完全相同相同,即局部截断误差为即局部截断误差为 。改进的改进的EulerEuler公式又可改写成公式又可改写成 )(2hO),(),()(21121211hKyxfKyxfKKKhyyiiiiii现在学习的是第32页，共81页 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点:都是用都是用f(x,y)f(x,y)在某些点上值的线性组合得出在某些点上值的线性组合得出y(xy(xi+1i+1)的近似值的近似值y yi+1i+1,而且增而且增加计算的次数加计算的次数f f(x x,y y)的次数的

26、次数,可提高截断误差的阶。如可提高截断误差的阶。如欧拉公式欧拉公式:每步计算一次每步计算一次f f(x x,y y)的值的值,为一阶方法。改为一阶方法。改进欧拉公式需计算两次进欧拉公式需计算两次f f(x x,y y)的值，它是二阶方法。它的的值，它是二阶方法。它的局部截断误差为局部截断误差为 。)(3hO现在学习的是第33页，共81页 于是可考虑用函数于是可考虑用函数f(x,y)在若干点上的函数值的在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在线性组合来构造近似公式，构造时要求近似公式在(xi i,yi i)处的处的Taylor展开式与解展开式与解y(x)在在xi i处的处

27、的Taylor展展开式的前面几项重合，从而使近似公式达到所需要的开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导阶数。既避免求偏导,又提高了计算方法精度的阶数又提高了计算方法精度的阶数。或者说。或者说,在在 这一步内多这一步内多预报几个点的斜率值，然后将其加预报几个点的斜率值，然后将其加权平均作为平均斜率权平均作为平均斜率，则可构造出更高精度的计算，则可构造出更高精度的计算格式，这就是龙格格式，这就是龙格库塔（库塔（Runge-Kutta）法的基）法的基本思想。本思想。1,iixx现在学习的是第34页，共81页9.3.2 二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi

28、 i和和 ,以该两点处的斜以该两点处的斜率值率值k1 1和和k2 2的加权平均的加权平均(或称为线性组合或称为线性组合)来求取平均斜来求取平均斜率率k*的近似值的近似值K，即，即 1,iixxphxxipi2211kkK式中式中:k1 1为为xi i点处的切线斜率值，点处的切线斜率值，k2 2为为 点处的切线斜率值点处的切线斜率值,比照改进的欧拉法比照改进的欧拉法,将将 视为视为 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12phkyphxfkii对常微分方程初值问题对常微分方程初值问题(9.1)(9.1)式的解式的解 y=y(x),),根据微分中值定根据微分中

29、值定理，存在点理，存在点 ，使得，使得 ),(1iixx现在学习的是第35页，共81页)(,()(yfyK式中式中 K K可看作是可看作是y=y(x)y=y(x)在区间在区间 上的平均斜率。所以可上的平均斜率。所以可得计算公式为：得计算公式为：1,iixxhKxyxyii)()(1)()(2211kkhxyi（9.14）将将y(xy(xi i)在在x=xx=xi i处进行二阶处进行二阶TaylorTaylor展开：展开：)()(!2)()()(321hOxyhxyhxyxyiiii（9 9.15）)()()(11iiiixxyxyxy也即也即 hKxyxyii)()(1（9.13）现在学习的是

30、第36页，共81页将将 在在x=xx=xi i处进行一阶处进行一阶TaylorTaylor展开：展开：),()(12phkyphxfphxykiii)(),(),(),(),(22hOyxfyxfyxfphyxfkiiyiiiixii)()()(2hOxyphxyii 将以上结果代入（将以上结果代入（9.149.14）得：）得：)()()(22111kkhxyxyii)()()()()(221hOxyphxyxyhxyiiii)()()()()(32221hOxyphxyhxyiii（9.16）对式对式(9.15)(9.15)和和(9.16)(9.16)进行比较系数后可知进行比较系数后可知,只

31、要只要 211221p（9.17）成立成立,格式格式(9.14)(9.14)的局部截断误差就等于的局部截断误差就等于)(3hO有有2 2阶阶精度精度现在学习的是第37页，共81页式式(9.17)(9.17)中具有三个未知量中具有三个未知量,但只有两个方程但只有两个方程,因而有无因而有无穷多解。若取穷多解。若取 ,则则p p=1=1，这是无穷多解中的一个，这是无穷多解中的一个解，将以上所解的值代入式解，将以上所解的值代入式(9.14)(9.14)并改写可得并改写可得 2121),(),()(21121211hkyxfkyxfkkkhyyiiiiii 不难发现，上面的格式就是改进的欧拉格式。凡不难

32、发现，上面的格式就是改进的欧拉格式。凡满足条件式（满足条件式（9.179.17）有一簇形如上式的计算格式，这）有一簇形如上式的计算格式，这些格式统称为二阶龙格些格式统称为二阶龙格库塔格式。因此改进的欧拉格库塔格式。因此改进的欧拉格式是众多的二阶龙格式是众多的二阶龙格库塔法中的一种特殊格式。库塔法中的一种特殊格式。现在学习的是第38页，共81页若取若取 ,则则 ，此时二阶龙格，此时二阶龙格-库塔库塔法的计算公式为法的计算公式为 0121,12p)2,(),(1212121khyxfkyxfkhkyyiiiiii1,2,1,0ni 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式

33、中库塔法。式中 为区间为区间 的中点。的中点。21ix1,iixx现在学习的是第39页，共81页9.3.3 三阶龙格三阶龙格-库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 pix)1(qpqhxxiqi并用三个点并用三个点 ,的斜率的斜率k1 1,k2 2,k3 3加权平均加权平均得出平均斜率得出平均斜率k*的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式:ixpixqix11 12233121(,)(,)iiiiiiyyhkkkkf x ykf xph yphk（9.18）为了预报点为了预报点 的斜率值的斜率值k3 3,在区间在区间 内

34、有两内有两个斜率值个斜率值k1 1和和k2 2可以用可以用,可将可将k1 1,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率,从而得到从而得到 的预报值的预报值 qixqiixx,qiixx,)(qixyqiy1122i qiyyqhkk现在学习的是第40页，共81页于是可得于是可得),(3qiqiyxfk运用运用Taylor展开方法选择参数展开方法选择参数 ,可以使格式可以使格式(9.18)的局部截断误差为的局部截断误差为 ,即具有三阶精度，这类格式即具有三阶精度，这类格式统称为统称为三阶龙格三阶龙格库塔方法库塔方法。下列是其中的一种，称为。下列是其中的一种，称为库库塔（塔（Ku

35、tta）公式。）公式。12212,p q )(4hO)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（9.19）现在学习的是第41页，共81页121231232221233132611pqpqpq 1141123122666,1,1,2pq 现在学习的是第42页，共81页9.3.4 四阶龙格四阶龙格库塔法库塔法 如果需要再提高精度，用类似上述的处理方法，只需在如果需要再提高精度，用类似上述的处理方法，只需在区间区间 上用四个点处的斜率加权平均作为平均斜率上用四个点处的斜率加权平均作为平均斜率k*的近似的近似值，构成一系列四

36、阶龙格值，构成一系列四阶龙格库塔公式。具有四阶精度，即局部库塔公式。具有四阶精度，即局部截断误差是截断误差是 。由于推导复杂，这里从略，只介绍最常用的一种由于推导复杂，这里从略，只介绍最常用的一种四阶经典四阶经典龙格龙格库塔公式库塔公式。qiixx,)(5hO)22(6),()2,()2,(),(43211314221312121kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（9.20）现在学习的是第43页，共81页9.3.5 9.3.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1)计算步骤计算步骤 输入输入 ,h,Nh,N 使用龙格使用龙格库塔公式

37、（库塔公式（9.209.20）计算出）计算出y y1 1 输出输出 ，并使，并使 转到转到 直至直至n n N N 结束。结束。10,xx11,yx0101,yyxx现在学习的是第44页，共81页(2)(2)四阶龙格四阶龙格库塔算法流程图库塔算法流程图 开 始 输 入 x0,y0,h,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 现在

38、学习的是第45页，共81页 程序实现程序实现(四阶龙格四阶龙格-库塔法计库塔法计 算常微分方程初值问题算常微分方程初值问题)例例9.4 9.4 取步长取步长h=0.2h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 1)0(2yxyy10 x解解:由四阶龙格由四阶龙格-库塔公式可得库塔公式可得 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现在学习的是第46页，共81页040811

39、.1)22(62.0432101kkkkyy可同样进行其余可同样进行其余y yi i的计算。本例方程的解为的计算。本例方程的解为，从表中看到所求的数值解具有，从表中看到所求的数值解具有4位有效数字。位有效数字。2xey 龙格龙格库塔方法的库塔方法的推导基于推导基于Taylor展开方法，因展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格性差，那么，使用四阶龙格库塔方法求得的数值解，库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法其精度可能反而不如改进的欧拉方法。在实际计算时，。在实际计算时，应当针对问题的具体特

40、点选择合适的算法。应当针对问题的具体特点选择合适的算法。现在学习的是第47页，共81页9.3.6 变步长的龙格变步长的龙格-库塔法库塔法 在微分方程的数值解中，选择适当的步长是非常重要在微分方程的数值解中，选择适当的步长是非常重要的。单从每一步看，的。单从每一步看，步长越小，截断误差就越小；但随着步长越小，截断误差就越小；但随着步长的缩小，在一定的求解区间内所要完成的步数就增加了步长的缩小，在一定的求解区间内所要完成的步数就增加了。这样会引起计算量的增大，并且会引起舍入误差的大量积。这样会引起计算量的增大，并且会引起舍入误差的大量积累与传播。因此微分方程数值解法也有选择步长的问题。累与传播。因

41、此微分方程数值解法也有选择步长的问题。以经典的四阶龙格以经典的四阶龙格-库塔法库塔法(9.20)为例。从节点为例。从节点x xi i出出发，先以发，先以h为步长求出一个近似值，记为为步长求出一个近似值，记为 ，由于局部截，由于局部截断误差为断误差为 ，故有，故有)(1hiy)(5hO5)(11)(chyxyhii当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。现在学习的是第48页，共81页然后将步长折半然后将步长折半,即以为即以为 步长步长,从节点从节点x xi i出发出发,跨两跨两步到节点步到节点x xi+1i+1,再求得一个近似值再求得一个近似

42、值 ,每跨一步的截断误每跨一步的截断误差是差是 ,因此有因此有2h)2(1hiy52 hc5)2(1122)()(hcxyxyhii这样这样 161)()()(11)2(11hiihiiyxyyxy)(151)()(1)2(1)2(11hihihiiyyyxy由此可得由此可得 这表明以这表明以 作为作为 的近似值，其误差可用步长的近似值，其误差可用步长折半前后两次计算结果的偏差折半前后两次计算结果的偏差 )2(1hiy)(1ixy)(1)2(1hihiyy来判断所选步长是否适当来判断所选步长是否适当现在学习的是第49页，共81页当要求的数值精度为当要求的数值精度为时：时：（1 1）如果）如果，

43、反复将步长折半进行计算，直至，反复将步长折半进行计算，直至为止为止,并取其最后一次步长的计算结果作为并取其最后一次步长的计算结果作为 （2 2）如果）如果为为止，并以上一次步长的计算结果作为止，并以上一次步长的计算结果作为 。这种通过步长加倍或折半来处理步长的方法称为变步长法这种通过步长加倍或折半来处理步长的方法称为变步长法。表面上看，为了选择步长，每一步都要反复判断。表面上看，为了选择步长，每一步都要反复判断，增加，增加了计算工作量，但在方程的解了计算工作量，但在方程的解y(x)y(x)变化剧烈的情况下，总的变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的。计算工作量得到减少，结果还

44、是合算的。1iy1iy现在学习的是第50页，共81页 9.4 算法的稳定性及收敛性算法的稳定性及收敛性9.4.1 稳定性稳定性 稳定性在微分方程的数值解法中是一个非常重要稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程中，存在格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误差等引起的着各种计算误差，这些计算误差如舍入误差等引起的扰动，在传播过程中，可能会大量积累，对计算结果扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。

45、这就涉及到算法稳定性问题。的准确性将产生影响。这就涉及到算法稳定性问题。现在学习的是第51页，共81页 当在某节点上当在某节点上x xi i的的y yi i值有大小为值有大小为的扰动时，如果在的扰动时，如果在其后的各节点其后的各节点 上的值上的值y yi i产生的偏差都不大于产生的偏差都不大于，则称这种方法是稳定的。，则称这种方法是稳定的。稳定性不仅与算法有关，而且与方程中函数稳定性不仅与算法有关，而且与方程中函数f(x,y)f(x,y)也有关，讨论起来比较复杂。也有关，讨论起来比较复杂。为简单起见，通常只针对模为简单起见，通常只针对模型方程型方程)(ijxj)0(yy来讨论。一般方程若局部线

46、性化，也可化为上述形式。来讨论。一般方程若局部线性化，也可化为上述形式。模型方程相对比较简单，若一个数值方法对模型方程是模型方程相对比较简单，若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程都稳定，但若某稳定的，并不能保证该方法对任何方程都稳定，但若某方法对模型方程都不稳定，也就很难用于其他方程的求方法对模型方程都不稳定，也就很难用于其他方程的求解。解。现在学习的是第52页，共81页先考察显式先考察显式EulerEuler方法的稳定性。模型方程方法的稳定性。模型方程的的EulerEuler公式为公式为 yy)(),(1iiiiiiyhyyxhfyy),2,1,0()1(iyhi将上

47、式反复递推后，可得将上式反复递推后，可得 011)1(yhyii),2,1()1(00iyyhyiii或或h1式中式中现在学习的是第53页，共81页 要使要使y yi i有界，其充要条件为有界，其充要条件为 111h即即 由于由于 ，故有，故有 020 h（9.269.26）可见，如欲保证算法的稳定，显式可见，如欲保证算法的稳定，显式EulerEuler格式的步格式的步长长h h的选取要受到式（的选取要受到式（9.269.26）的限制。）的限制。的绝对值越大的绝对值越大，则限制的，则限制的h h值就越小。值就越小。用隐式用隐式EulerEuler格式，对模型方程格式，对模型方程 的计算公式为，

48、可化为的计算公式为，可化为)(11iiiyhyyiiyhy111现在学习的是第54页，共81页由于由于 ,则恒有则恒有 ,故恒有故恒有 0111hiiyy1 因此，隐式因此，隐式EulerEuler格式是绝对稳定的（无条件稳定）的格式是绝对稳定的（无条件稳定）的（对任何（对任何h0h0）。）。9.4.2 9.4.2 收敛性收敛性 常微分方程初值问题的求解，是将微分方程转化为差分方常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值程来求解，并用计算值y yi i来近似替代来近似替代y(xy(xi i)，这种近似替代是否合理这种近似替代是否合理，还 须 看 分 割 区 间，还

49、须 看 分 割 区 间 的 长 度的 长 度 h h 越 来 越 小 时，即越 来 越 小 时，即 时，时，是否成立。若成立，则称该方法是收敛的是否成立。若成立，则称该方法是收敛的，否则称为不收敛。，否则称为不收敛。这里仍以这里仍以EulerEuler方法为例，来分析其收敛性。方法为例，来分析其收敛性。EulerEuler格式格式iixx,101iixxh)(iixyy),(1iiiiyxhfyy现在学习的是第55页，共81页设设 表示取表示取 时时,按按Euler公式的计算结果公式的计算结果,即即1iy)(iixyy)(,()(1iiiixyxhfxyyEuler方法局部截断误差为：方法局部

50、截断误差为：)()(2)(1211 iiiixxyhyxy设有常数设有常数 ，则，则 )(max21xycbxa 211)(chyxyii（9.27）总体截断误差总体截断误差 1111111)()(iiiiiiiyyyxyyxy)(,()(),(11iiiiiiiixyxhfxyyxhfyyy),()(,()(iiiiiiyxfxyxfhyxy（9.28）又又 由于由于f(x,y)f(x,y)关于关于y y满足李普希茨条件满足李普希茨条件,即即 现在学习的是第56页，共81页iiiiiiyxyLyxfxyxf)(),()(代入代入(9.28)上式，有上式，有 iiiiyxyhLyy)()1(1