第七章常微分方程数值解

1、计算方法河北大学电信学院第七章第七章 常微分方程的数值解法常微分方程的数值解法 7.1 7.1 引言引言 包含自变量、未知函数及未知函数的导数或微包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中分的方程称为微分方程。在微分方程中, , 自变量的自变量的个数只有一个个数只有一个, , 称为常微分方程称为常微分方程. .。自变量的个数。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数方程的阶数。如果未知函数y

2、y及其各阶导数及其各阶导数都是一次的都是一次的, ,则称它是线性的则称它是线性的, ,否则称为非线性的。否则称为非线性的。 )(,nyyy 在高等数学中，对于常微分方程的求解，给出在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。有限的，大多数的常微分方程是不可能给出解析解。 譬如譬如 22yxy

3、这个一阶微分方程就不能用初等函数及其积这个一阶微分方程就不能用初等函数及其积分来表达它的解。分来表达它的解。 再如，方程再如，方程 1)0(yyy的解的解 , ,虽然有表可查虽然有表可查, ,但对于表但对于表上没有给出上没有给出 的值的值, ,仍需插值方法来仍需插值方法来计算计算xey xe从实际问题当中归纳出来的微分方程，通常主要依从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程靠数值解法来解决。本章主要讨论一阶常微分方程初值问题初值问题 00)(),(yxyyxfy ( 7.1 ) 在区间在区间a x ba x b上的数值解法上的数值解法。 可以证明

4、可以证明, ,如果函数在带形区域如果函数在带形区域 r=axb,r=axb,-y y内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹( (lipschitzlipschitz) )条件，即存在常数条件，即存在常数l(l(它与它与x,yx,y无关无关) )使使 2121),(),(yylyxfyxf对对r内任意两个内任意两个 都成立都成立, ,则方程则方程( 7.1 )的解的解 在在 a, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 7.2 数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题( (7.1)式的数值解法，就式的数值解法，就是要算出

5、精确解是要算出精确解y(x)y(x)在区间在区间 a,b 上的一系列离散节上的一系列离散节点点 处的函数值处的函数值 的近似值的近似值。相邻两个节点的间距。相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定h为定数，为定数，称为称为定步长定步长，这时节点可表示为，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。出离散节点的数值解。 bxxxxann110)(,),(),(10nxyxyxynyyy,10iixxh1niihxxi, 2 , 1 ,0 对常

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

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

8、+1时，除用到时，除用到x xi i+1+1,x,xi i和和y yi i以外，以外，还要用到还要用到，即前面，即前面k步的值，此类方法称为步的值，此类方法称为多步法多步法；其代表；其代表是亚当斯法。是亚当斯法。 00)(),(yxyyxfy), 2 , 1( ,kpyxpipi7.3 欧拉（欧拉（euler）法）法7.3.1 euler公式公式 欧拉（欧拉（euler）方法是解初值问题的最）方法是解初值问题的最简单的数值方法。初值问题简单的数值方法。初值问题的解的解y=y(x)y=y(x)代表通过点代表通过点 的一条称之为的一条称之为微分方程的积分曲线。积分曲线上每一点微分方程的积分曲线。积

9、分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在在这点的值。这点的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf 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点点( (即点即点( (x

10、 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点的坐标。点的坐标。 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),),作积分曲线作积分曲

11、线y=y(x)y=y(x)的切线的切线交直线交直线x=xx=x2 2于于p2点点, ,切线切线 的斜率的斜率 = =直线方程为直线方程为21pp)(1xy),(11yxf)(,(1111xxyxfyy)(,(121112xxyxfyy当当 时时, ,得得 2xx 当当 时时, ,得得 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。对已求得点对已求得点以以 =

12、 = 为斜率作直线为斜率作直线 ),(nnnyxp)(nxy),(nnyxf)(,(nnnnxxyxfyy1nxx)(,(11nxnnnnxxyxfyynnyxy)(取取 从图形上看从图形上看, ,就获得了一条近似于曲线就获得了一条近似于曲线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通常取通常取 ( (常数常数),),则则euler法的计算格式法的计算格式 hhx

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

14、(,iixxdxxyxf(7.3) 用左矩形方法计算积分项用左矩形方法计算积分项 )(,)()(,11iiiixxxyxfxxdxxyxfii代入代入(7.3)(7.3)式式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到即可得到向前欧拉（向前欧拉（eulereuler）公式）公式 ),(1iiiiyxhfyy 由于数值积分的矩形方法精度很低，所以由于数值积分的矩形方法精度很低，所以欧拉（欧拉（eulereuler）公式当然很粗糙。）公式当然很粗糙。 例例7.1 用欧拉法解初值问题用欧拉法解初值问题 1)0()6 . 00(2yxxyyy取步长取步长h=0.2

15、 ,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， y2 =0.6144

16、，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 7.3.2 梯形公式梯形公式为了提高精度为了提高精度, ,对方程对方程 的两端在区间上的两端在区间上 积分得，积分得，改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii( 7.4 ) 代入代入(7.4)(7.4)式式, ,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 ),(),(2111iiiiiiyxfyxfhyy(

17、 7.5 ) 由于数值积分的梯形公式比矩形公式的精度高，由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（因此梯形公式（7.57.5）比欧拉公式）比欧拉公式( 7.2 )( 7.2 )的精度高的精度高一个数值方法。一个数值方法。 ),(),(2111iiiiiiyxfyxfhyy( 7.5 ) ( (7.5)式的右端含有未知的式的右端含有未知的y yi i+1+1, ,它是一个关于它是一个关于y yi i+1+1的函数方程的函数方程, ,这类数值方法称为这类数值方法称为隐式方法隐式方法。相。相反地反地, ,欧拉法是关于欧拉法是关于y yi i+1+1的一个直接的计算公式，的一个直接的计算

18、公式， 这类数值方法称为这类数值方法称为显式方法。显式方法。 7.3.3 两步欧拉公式两步欧拉公式 对方程对方程 的两端在区间上的两端在区间上 积分得积分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 7.6 ) 改用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到两即可得到两步欧拉公式步欧拉公式 ),(211iiiiyxhfyy ( 7.7 ) 前面介绍过的数值方法前面介绍过的数值

19、方法, ,无论是欧拉无论是欧拉方法方法, ,还是梯形方法，它们都是单步法还是梯形方法，它们都是单步法, ,其其特点是在计算特点是在计算y yi i+1+1时只用到前一步的信息时只用到前一步的信息y yi i; ;可是公式可是公式( (7.7)中除了中除了y yi i外外, ,还用到更前一步还用到更前一步的信息的信息y yi i-1-1, ,即调用了前两步的信息即调用了前两步的信息, ,故称其故称其为两步欧拉公式为两步欧拉公式 7.3.4 欧拉法的局部截断误差欧拉法的局部截断误差 衡量求解公式好坏的一个主要标准是求解公式的衡量求解公式好坏的一个主要标准是求解公式的精度精度, 因此引入局部截断误差

20、和阶数的概念。因此引入局部截断误差和阶数的概念。 定义定义7.1 在在yi准确的前提下准确的前提下, 即即 时时, 用数值用数值方法计算方法计算yi+1的误差的误差 , 称为该数值方法称为该数值方法计算时计算时yi+1的局部截断误差。的局部截断误差。 对于欧拉公式，假定对于欧拉公式，假定 ，则有，则有)(iixyy 11)(iiiyxyr)(iixyy )()()(,()(1iiiiiixyhxyxyxfhxyy而将真解而将真解y(x)在在xi处按二阶泰勒展开处按二阶泰勒展开 ),()(! 2)()()(121 iiiiixxyhxyhxyxy)(!2)(211yhyxyii 因此有因此有 定

21、义定义7.2 数值方法的局部截断误差为数值方法的局部截断误差为 ,则称这则称这种数值方法的阶数是种数值方法的阶数是p。步长。步长(h n 结束。结束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxx（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, y1 n + 1 n n = n ? x1 x0 y1 y0 结 束 n y ( (3)3)

22、程序实现程序实现( (见附录见附录a a-15 改进欧拉法计算常微改进欧拉法计算常微 分方程初值问题分方程初值问题 ) 例例7.2 7.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 1)0(2yyxyy区间为区间为 0,10,1 , ,取步长取步长h=0.1h=0.1 解解: : 改进欧拉法的具体形式改进欧拉法的具体形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyyxyyyyxyyy本题的精确解为本题的精确解为 , ,计算见计算见p p158158列表所示列表所示 xxy21)(例例7.3 对初值问题对初值问题 1)0(0yyy证明用梯形公式求得的近似解为证明用梯

23、形公式求得的近似解为 nnhhynhx 并证明当步长并证明当步长h h0 0时时, ,y yn n收敛于精确解收敛于精确解证明证明: : 解初值问题的梯形公式为解初值问题的梯形公式为 xe),(),(nnnnnnyxfyxfhyyyyxf),( 211nnnnyyhyy 整理成显式整理成显式 nnyhhy反复迭代反复迭代, ,得到得到 yhhyhhyhhyhhynnnnn.10ynnhhy22 由于由于 ，有，有 nhx xxxxhxhhhxhnhhhhhyeee2121lim22limlim222222000nnhhy22xnhy elim0 证毕证毕 7.4 7.4 龙格龙格- -库塔（库

24、塔（runge-kuttarunge-kutta）法）法7.4.1 7.4.1 龙格龙格- -库塔库塔( (runge-kuttarunge-kutta) )法的基本思想法的基本思想 eulereuler公式可改写成公式可改写成 ),(111iiiiyxfkhkyy则则yi i+1+1的表达式的表达式y( (xi i+1+1) )与的与的taylortaylor展开式的前两项展开式的前两项完全相同完全相同, ,即局部截断误差为即局部截断误差为 。 改进的改进的eulereuler公式又可改写成公式又可改写成 )(2ho),(),()(21121211hkyxfkyxfkkkhyyiiiiii

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

26、差为的值，它是二阶方法。它的局部截断误差为 。)(3ho 于是可考虑用函数于是可考虑用函数f( (x, ,y) )在若干点上的函在若干点上的函数值的线性组合来构造近似公式，构造时要求数值的线性组合来构造近似公式，构造时要求近似公式在近似公式在( (xi i, ,yi i) )处的处的taylor展开式与解展开式与解y( (x) )在在xi i处的处的taylor展开式的前面几项重合，从而展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导使近似公式达到所需要的阶数。既避免求偏导, ,又提高了计算方法精度的阶数。或者说又提高了计算方法精度的阶数。或者说, ,在在 这一步内多预报几个

27、点的斜率值，然后这一步内多预报几个点的斜率值，然后将其加权平均作为平均斜率，则可构造出更高将其加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格精度的计算格式，这就是龙格库塔（库塔（runge-kutta）法的基本思想。）法的基本思想。 1,iixx7.4.2 二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi i和和 , ,以该两点处的以该两点处的斜率值斜率值k1 1和和k2 2的加权平均的加权平均( (或称为线性组合或称为线性组合) )来求取平来求取平均斜率均斜率k* *的近似值的近似值k，即，即 1,iixxphxxipi2211kkk式中式中: :k1 1为为xi i

28、点处的切线斜率值，点处的切线斜率值， k2 2为为 点处的切线斜率值点处的切线斜率值, ,比照改进的欧比照改进的欧拉法拉法, ,将将 视为视为 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12phkyphxfkii对常微分方程初值问题对常微分方程初值问题(7.1)(7.1)式的解式的解 y= =y( (x),),根据微根据微分中值定理，存在点分中值定理，存在点 ，使得，使得 ),(1iixx)(,()(yfyk式中式中 k k可看作是可看作是y=y(x)y=y(x)在区间在区间 上的平均斜率。上的平均斜率。所以可得计算公式为：所以可得计算公式为： 1,iixxhk

29、xyxyii)()(1)()(2211kkhxyi（7.14） 将将y(xy(xi i) )在在x=xx=xi i处进行二阶处进行二阶taylortaylor展开：展开： )()(! 2)()()(321hoxyhxyhxyxyiiii （7.15） )()()(11iiiixxyxyxy也即也即 hkxyxyii)()(1（7.13）将将 在在x=xx=xi i处进行一阶处进行一阶taylortaylor展开：展开： ),()(12phkyphxfphxykiii)(),(),(),(),(22hoyxfyxfyxfphyxfkiiyiiiixii)()()(2hoxyphxyii 将以上结

30、果代入（将以上结果代入（7.147.14）得：）得： )()()(22111kkhxyxyii)()()()()(221hoxyphxyxyhxyiiii )()()()()(32221hoxyphxyhxyiii （7.16） 对式对式(7.15)(7.15)和和(7.16)(7.16)进行比较系数后可知进行比较系数后可知, ,只要只要 211221p（7.17） 成立成立, ,格式格式(7.14)(7.14)的局部截断误差就等于的局部截断误差就等于)(3ho有有2 2阶阶精度精度式式(7.17)(7.17)中具有三个未知量中具有三个未知量, ,但只有两个方程但只有两个方程, ,因而因而有无

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

32、种特殊库塔法中的一种特殊格式。格式。 若取若取 , ,则则 ，此时二阶龙格，此时二阶龙格- -库塔库塔法的计算公式为法的计算公式为 0121, 12p)2,(),(1212121khyxfkyxfkhkyyiiiiii1,2 , 1 , 0ni 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间 的中点。的中点。 21ix1,iixx7.4.3 三阶龙格三阶龙格- -库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 pix) 1(qpqhxxiqi并用三个点并用三个点 , , , , 的斜率的斜率k1 1, ,

33、k2 2, ,k3 3加权平均加权平均得出平均斜率得出平均斜率k* *的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式: : ixpixqix),(),()1 (1213211phkyphxfkyxfkkkkhyyiiiiii（7.18） 为了预报点为了预报点 的斜率值的斜率值k3 3, ,在区间在区间 内有两内有两个斜率值个斜率值k1 1和和k2 2可以用可以用, ,可将可将k1 1, ,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率, ,从而得到从而得到 的预报值的预报值 qixqiixx,qiixx,)(qixyqiy21)1 (kkqhyyiqi于是可得于是可得

34、 ),(3qiqiyxfk运用运用taylor展开方法选择参数展开方法选择参数 , ,可以使可以使格式格式( (7.18)的局部截断误差为的局部截断误差为 , ,即具有三阶精度，即具有三阶精度，这类格式统称为这类格式统称为三阶龙格三阶龙格库塔方法库塔方法。下列是其中。下列是其中的一种，称为的一种，称为库塔（库塔（kutta）公式）公式。 ,qp)(4ho)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（7.19） 7.4.4 四阶龙格四阶龙格库塔法库塔法 如果需要再提高精度，用类似上述的处理方法，如果需要再提高精度，用类似

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

36、阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1) 计算步骤计算步骤 输入输入 , ,h,nh,n 使用龙格使用龙格库塔公式（库塔公式（7.207.20）计算出）计算出y y1 1 输出输出 ，并使，并使 转到转到 直至直至n n n n 结束。结束。 10,xx11, yx0101,yyxx(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 k

37、3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 输 出 x1, y1 n + 1 n n = n ? x1 x0 y1 y0 结 束 n y (3)(3) 程序实现程序实现( (见附录见附录a a-16 四阶龙格四阶龙格- -库塔法计库塔法计 算常微分方程初值问题算常微分方程初值问题) ) 例例7.4 7.4 取步长取步长h=0.2h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 1)0(2yxyy10 x解解: : 由四阶龙格由四阶龙格- -库塔公式可得库塔公式可得 2 . 0, 1, 0,2),(00hyxxyyxf0),(001yxfk2 . 0

38、) 1 , 1 . 0()2,(10202fkhyxfkh204. 0)02. 1 , 1 . 0()2,(20203fkhyxfkh41632. 0)0408. 1 , 2 . 0(),(3004fhkyxfkh040811. 1)22(62 . 0432101kkkkyy可同样进行其余可同样进行其余y yi i的计算。本例方程的解为的计算。本例方程的解为，数值解，数值解y yi i与准确解与准确解y(xy(xi i) )的对照表见教材的对照表见教材p p163163所示，所示，从表中看到所求的数值解具有从表中看到所求的数值解具有4位有效数字。位有效数字。 2xey 龙格龙格库塔方法的推导基

39、于库塔方法的推导基于taylor展开方法，展开方法，因而它要求所求的解具有较好的光滑性。如果解的因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格光滑性差，那么，使用四阶龙格库塔方法求得的库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应当针对问题的具体特点选择合适的实际计算时，应当针对问题的具体特点选择合适的算法。算法。 7.4.6 变步长的龙格变步长的龙格-库塔法库塔法 在微分方程的数值解中，选择适当的步长是非常在微分方程的数值解中，选择适当的步长是非常重要的。单从每一步看，步长越小，截断误差就越重要的

40、。单从每一步看，步长越小，截断误差就越小；但随着步长的缩小，在一定的求解区间内所要小；但随着步长的缩小，在一定的求解区间内所要完成的步数就增加了。这样会引起计算量的增大，完成的步数就增加了。这样会引起计算量的增大，并且会引起舍入误差的大量积累与传播。因此微分并且会引起舍入误差的大量积累与传播。因此微分方程数值解法也有选择步长的问题。方程数值解法也有选择步长的问题。 以经典的四阶龙格以经典的四阶龙格- -库塔法库塔法( (7.20)为例。从节点为例。从节点x xi i出发，先以出发，先以h为步长求出一个近似值，记为为步长求出一个近似值，记为 ，由于局部截断误差为由于局部截断误差为 ，故有，故有

41、)(1hiy)(5ho5)(11)(chyxyhii当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。然后将步长折半然后将步长折半, ,即以为即以为 步长步长, ,从节点从节点x xi i出发出发, ,跨跨两步到节点两步到节点x xi i+1+1, ,再求得一个近似值再求得一个近似值 , ,每跨一步的每跨一步的截断误差是截断误差是 , ,因此有因此有2h)2(1hiy52 hc5)2(1122)()(hcxyxyhii这样这样 161)()()(11)2(11hiihiiyxyyxy)(151)()(1)2(1)2(11hihihiiyyyxy由此可

42、得由此可得 这表明以这表明以 作为作为 的近似值，其误差可用步的近似值，其误差可用步长折半前后两次计算结果的偏差长折半前后两次计算结果的偏差 )2(1hiy)(1ixy)(1)2(1hihiyy来判断所选步长是否适当来判断所选步长是否适当当要求的数值精度为当要求的数值精度为时：时： （1 1）如果）如果，反复将步长折半进行计算，直，反复将步长折半进行计算，直至至为止为止, ,并取其最后一次步长的计算结果作为并取其最后一次步长的计算结果作为 （2 2）如果）如果为止，并以上一次步长的计算结果作为为止，并以上一次步长的计算结果作为 。 这种通过步长加倍或折半来处理步长的方法称为这种通过步长加倍或折

43、半来处理步长的方法称为变步长法。表面上看，为了选择步长，每一步都要变步长法。表面上看，为了选择步长，每一步都要反复判断反复判断，增加了计算工作量，但在方程的解，增加了计算工作量，但在方程的解y(x)y(x)变化剧烈的情况下，总的计算工作量得到减少，结变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的。果还是合算的。1iy1iy7.5 亚当姆斯方法亚当姆斯方法7.5.1 亚当姆斯格式亚当姆斯格式 龙格龙格-库塔方法是一类重要算法，但这类库塔方法是一类重要算法，但这类算法在每一步都需要先预报几个点上的斜率算法在每一步都需要先预报几个点上的斜率值，计算量比较大。考虑到计算值，计算量比较大。考虑

44、到计算yi+1之前已得之前已得出一系列节点上出一系列节点上 的斜率值，的斜率值，能否利用这些已知值来减少计算量呢？能否利用这些已知值来减少计算量呢？ 这就是亚当姆斯（这就是亚当姆斯（adams）方法的设计）方法的设计思想。思想。 11,xxxii 设用设用x xi i,x,xi i+1+1两点的斜率值加权平均作为区间两点的斜率值加权平均作为区间 上的平均斜率，有计算格式上的平均斜率，有计算格式 1,iixx),(),()1 (11111iiiiiiiiiiyxfyyxfyyyhyy（7.21） 选取参数选取参数，使格式（，使格式（7.217.21）具有二阶精度。）具有二阶精度。 将将 在在x

45、xi i处处taylor展开展开 1iy)()(! 21)(321hohyhyyyiii 代入计算格式代入计算格式(7.21)(7.21)化简化简, ,并假设并假设, , 因此有因此有 )(),(11iiiixyyxyy)()()()(321hoxyhxyhxyyiiii 与与y(xi+1)在在xi处的处的taylor展开式展开式 )()(21)()()(321hoxyhxyhxyxyiiii 相比较相比较, 需取需取 21才使格式才使格式(7.21)具有二阶精度。这样导出的计算格式具有二阶精度。这样导出的计算格式 )3(211iiiiyyhyy称之为二阶亚当姆斯格式。类似地可以导出三阶亚称之

46、为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式当姆斯格式。 )51623(12211iiiiiyyyhyy和和四阶亚当姆斯格式。四阶亚当姆斯格式。 )9375955(243211iiiiiiyyyyhyy( (7.22) 这里和下面均记这里和下面均记 ),(kikikiyxfy 上述几种亚当姆斯格式都是显式的，算法比较上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点简单，但用节点 的斜率值来预报的斜率值来预报区间区间 上的平均斜率是个外推过程，效果不够上的平均斜率是个外推过程，效果不够理想。为了进一步改善精度，变外推为内插，即理想。为了进一步改善精度，变外推为内插，即增加节点增加节点

47、xi+1的斜率值来得出的斜率值来得出 上的平均斜率。上的平均斜率。譬如考察形如譬如考察形如 11,xxxii1,iixx1,iixx),(),()1 (11111iiiiiiiiiiyxfyyxfyyyhyy（7.23） 的隐式格式的隐式格式, ,设设( (7.23)右端的右端的taylor展开有展开有 )(),(11iiiixyyxyy)()()1 ()()(321hoxyhxyhxyyiiii 可见要使格式可见要使格式(7.23)(7.23)具有二阶精度具有二阶精度, ,需令需令 , ,这样就可构造二阶隐式亚当姆斯格式这样就可构造二阶隐式亚当姆斯格式 21)(211iiiiyyhyy其实是

48、梯形格式。类似可导出三阶隐式亚当姆斯格其实是梯形格式。类似可导出三阶隐式亚当姆斯格式式 )85(12111iiiiiyyyhyy和四阶隐式亚当姆斯格式和四阶隐式亚当姆斯格式 )5199(242111iiiiiiyyyyhyy（7.247.24） 7.5.2 亚当姆斯预报亚当姆斯预报-校正格式校正格式 参照改进的欧拉格式的构造方法，以四阶亚当参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式（姆斯为例，将显式（7.22）和隐式（）和隐式（7.24）相结合，）相结合，用显式公式做预报，再用隐式公式做校正，可构成用显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报亚当姆斯预报-校正格式校正

49、格式 )9375955(243211iiiiiiyyyyhyy),(111iiiyxfy)5199 (242111iiiiiiyyyyhyy),(111iiiyxfy（7.257.25） 预报：预报： 校正：校正： 这种预报这种预报- -校正格式是四步法，它在计校正格式是四步法，它在计算算y yi i+1+1时不但用到前一步的信息时不但用到前一步的信息 ，而且，而且要用到再前面三步的信息要用到再前面三步的信息 ，因，因此它不能自行启动。在实际计算时，可借助此它不能自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格于某种单步法，譬如四阶龙格库塔法提供库塔法提供开始值开始值 。 iiyy,3

50、21,iiiyyy321,yyy例例7.5 取步长取步长h=0.1,h=0.1,用亚当姆斯预报用亚当姆斯预报- -校正公式求解校正公式求解 初值问题初值问题 1)0(yyxy10 x的数值解。的数值解。 解解: : 用四阶龙格用四阶龙格- -库塔公式求出发值库塔公式求出发值 , ,计算得：计算得：1 . 0, 1, 0,),(00hyxyxyxf321,yyy399717. 1,242805. 1,110342. 1321yyy表中的表中的 , ,y yi i和和y(xy(xi i) )分别为预报值、校正值和准确分别为预报值、校正值和准确解解( ),( ),以比较计算结果的精度。以比较计算结果

51、的精度。 再使用亚当姆斯预报再使用亚当姆斯预报- -校正公式校正公式(7.25),(7.25),见教材见教材p p166166列表算得其余的计算结果列表算得其余的计算结果iy12xeyx 7.6 算法的稳定性及收敛性算法的稳定性及收敛性 7.6.1 稳定性稳定性 稳定性在微分方程的数值解法中是一个非常稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是重要的问题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误中，存在着各种计算误差，这些计算误差如舍

52、入误差等引起的扰动，在传播过程中，可能会大量积累，差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。稳定性问题。 当在某节点上当在某节点上x xi i的的y yi i值有大小为值有大小为的扰动时，的扰动时，如果在其后的各节点如果在其后的各节点 上的值上的值y yi i产生的偏差产生的偏差都不大于都不大于，则称这种方法是稳定的。，则称这种方法是稳定的。 稳定性不仅与算法有关，而且与方程中函数稳定性不仅与算法有关，而且与方程中函数f(x,y)f(x,y)也有关，讨论起来比较复杂。也有关，讨论起来比较复杂

53、。为简单起见，为简单起见，通常只针对模型方程通常只针对模型方程)(ijxj)0(yy来讨论。一般方程若局部线性化，也可化为上述形来讨论。一般方程若局部线性化，也可化为上述形式。模型方程相对比较简单，若一个数值方法对模式。模型方程相对比较简单，若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程都型方程是稳定的，并不能保证该方法对任何方程都稳定，但若某方法对模型方程都不稳定，也就很难稳定，但若某方法对模型方程都不稳定，也就很难用于其他方程的求解。用于其他方程的求解。先考察显式先考察显式eulereuler方法的稳定性。模型方程方法的稳定性。模型方程的的eulereuler公式为公式为 y

54、y)(),(1iiiiiiyhyyxhfyy), 2 , 1 , 0()1 (iyhi将上式反复递推后，可得将上式反复递推后，可得 011)1(yhyii), 2 , 1()1 (00iyyhyiii或或h1式中式中 要使要使y yi i有界，其充要条件为有界，其充要条件为 111h即即 由于由于 ，故有，故有 020 h（7.267.26） 可见，如欲保证算法的稳定，显式可见，如欲保证算法的稳定，显式eulereuler格式格式的步长的步长h h的选取要受到式（的选取要受到式（7.267.26）的限制。）的限制。 的绝的绝对值越大，则限制的对值越大，则限制的h h值就越小。值就越小。 用隐式

55、用隐式eulereuler格式，对模型方程格式，对模型方程 的计算公式为，可化为的计算公式为，可化为)(11iiiyhyyiiyhy111由于由于 , ,则恒有则恒有 , ,故恒有故恒有 0111hiiyy1 因此，隐式因此，隐式eulereuler格式是绝对稳定的（无条件稳格式是绝对稳定的（无条件稳定）的（对任何定）的（对任何h0h0）。）。 7.6.2 7.6.2 收敛性收敛性 常微分方程初值问题的求解，是将微分方程转化常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值为差分方程来求解，并用计算值y yi i来近似替代来近似替代y(xy(xi i) )，这种近似替代是否

56、合理，还须看分割区间这种近似替代是否合理，还须看分割区间 的的长度长度h h越来越小时，即越来越小时，即 时，时， 是是否成立。若成立，则称该方法是收敛的，否则称为否成立。若成立，则称该方法是收敛的，否则称为不收敛。不收敛。 这里仍以这里仍以eulereuler方法为例，来分析其收敛性。方法为例，来分析其收敛性。eulereuler格式格式iixx,101iixxh)(iixyy),(1iiiiyxhfyy设设 表示取表示取 时时, 按按euler公式的计算结果公式的计算结果, 即即1iy)(iixyy )(,()(1iiiixyxhfxyyeuler方法局部截断误差为：方法局部截断误差为：)

57、()(2)(1211 iiiixxyhyxy设有常数设有常数 ，则，则 )(max21xycbxa 211)(chyxyii（7.27） 总体截断误差总体截断误差 1111111)()(iiiiiiiyyyxyyxy)(,()(),(11iiiiiiiixyxhfxyyxhfyyy),()(,()(iiiiiiyxfxyxfhyxy（7.28） 又又 由于由于f(x,y)f(x,y)关于关于y y满足李普希茨条件满足李普希茨条件, ,即即 iiiiiiyxylyxfxyxf)(),()(代入代入(7.28)上式，有上式，有 iiiiyxyhlyy)()1(11ihl)1( 再利用式（再利用式（

58、7.27），式（），式（7.28） 1111111)()(iiiiiiiyyyxyyxy2)1 (chhli21)1 (chhlii即即 上式反复递推后，可得上式反复递推后，可得 1020)1 ()1 (ikkiihlchhl1)1 ()1 (0iihllchhl（7.29） 设设 （t t为常数）为常数） tihxxi0hlehl 1tlihlieehl)1 (因为因为 所以所以 把上式代入式（把上式代入式（7.297.29），得），得 ) 1(0tltlielche若不计初值误差，即若不计初值误差，即 ，则有，则有 00helctli)1(（7.30） 式式(7.30)(7.30)说明说明

59、, ,当时当时 , , ,即即 , ,所以所以eulereuler方法是收敛的，且其收敛速度为方法是收敛的，且其收敛速度为 ，即具，即具有一阶收敛速度。同时还说明有一阶收敛速度。同时还说明eulereuler方法的整体截断方法的整体截断误差为误差为 ，因此算法的精度为一阶。，因此算法的精度为一阶。 0h0i)(iixyy )(ho)(ho7.7 一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程 我们已介绍了一阶常微分方程初值问题的各种我们已介绍了一阶常微分方程初值问题的各种数值解法，对于一阶常微分方程组，可类似得到各数值解法，对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化

60、为一阶常微分方种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。程组来求解。 7.7.1 一阶常微分方程组一阶常微分方程组对于一阶常微分方程组的初值问题对于一阶常微分方程组的初值问题 0000)(),()(),(zxzzyxgzyxyzyxfy（7.31） 可以把单个方程可以把单个方程 中的中的f 和和y看作向量看作向量来处理，这样就可把前面介绍的各种差分算法推广来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。到求一阶方程组初值问题中来。 ),(yxfy 设设 为节点上的近似解，为节点上的近似解，则有改进的则有改进的eulereuler格式为格式为 iiizyii