第9章常微分方程初值问题数值解法

1、2022-4-251第第9章章 常微分方程初值问题数值解法常微分方程初值问题数值解法1. Euler 公式公式2. 改进的欧拉公式改进的欧拉公式3. 龙格龙格库塔法库塔法4. 亚当斯法亚当斯法5. 算法的稳定性及收敛性算法的稳定性及收敛性2022-4-2529.1 9.1 引言引言 包含自变量、未知函数及未知函数的导数或微包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中分的方程称为微分方程。在微分方程中, , 自变量的自变量的个数只有一个个数只有一个, , 称为常微分方程。自变量的个数为称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方两个或两

2、个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数程的阶数。如果未知函数 y 及其各阶导数及其各阶导数都是一次的都是一次的, ,则称它是线性微分方程则称它是线性微分方程, ,否则称为非线否则称为非线性的。性的。 )(,nyyy 2022-4-253 在高等数学中，对于常微分方程的求解，给出了在高等数学中，对于常微分方程的求解，给出了一些典型方程求一些典型方程求解析解解析解的基本方法，如分离变量法、的基本方法，如分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程常系数齐次线性方程的解法、常系数

3、非齐次线性方程的解法等。但能求解析解的方程是有限的的解法等。但能求解析解的方程是有限的，大多数，大多数的的常微分方程常微分方程不能给出解析解。不能给出解析解。 譬如：譬如： 22yxy 这个一阶微分方程就不能用初等函数来表达这个一阶微分方程就不能用初等函数来表达它的解。它的解。 2022-4-254再如，方程再如，方程 1)0(yyy虽然有解虽然有解 , ,仍需插值方法来计算各点的值仍需插值方法来计算各点的值xey 2022-4-255 从实际问题当中归纳出来的微分方程，通常主从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分要依靠数值解法来解决。本章主要讨论

4、一阶常微分方程初值问题方程初值问题 00)(),(yxyyxfy ( 9.1 ) 在区间在区间 a x b 上的数值解法上的数值解法。 可以证明可以证明, ,如果函数在带形区域如果函数在带形区域 R: aR: axb,b,-y 内连续，且关于内连续，且关于 y 满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，即存在常数条件，即存在常数 L(L(它与它与x,yx,y无关无关) )使使 2121),(),(yyLyxfyxf 对对R内任意内任意x及两个及两个 都成立都成立, ,则方程则方程( 9.1 )的的解解 在在 a, b 上上存在且唯一存在且唯一。 21, yy)(

5、xyy 2022-4-256数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题( (9.1)式的数值解法，就式的数值解法，就是要算出精确解是要算出精确解 y( (x) )在区间在区间 a,b 上的一系列离散节上的一系列离散节点点 处的函数值处的函数值 的近似值的近似值相邻两个节点的间距相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定 h 为定数，为定数，称为称为定步长定步长，这时节点可表示为，这时节点可表示为 数值解法需要把连续性的问题加以离散化，数值解法需要把连续性的问题加以离散化，从而求出离散节

6、点的数值解。从而求出离散节点的数值解。 bxxxxann 110)(,),(),(10nxyxyxy.,10nyyyiiixxh 1niihxxi, 2 , 1,0 2022-4-257 对常微分方程数值解法的基本出发点就是离散对常微分方程数值解法的基本出发点就是离散化。其数值解法有化。其数值解法有两个基本特点，两个基本特点，11它们都采用它们都采用“步进式步进式”，即求解过程顺着节点排列的次序一步，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知一步地向前推进，描述这类算法，要求给出用已知信息信息 计算计算 的递推公式。的递推公式。22建建立这类递推公式的基本方法

7、是在这些节点上用数值立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等方法，对初值问题中积分、数值微分、泰勒展开等方法，对初值问题中的的导数导数 进行不同的离散化进行不同的离散化处理处理。 iyyy,201 iyy 00)(),(yxyyxfy11(,)kkkkxxxxy dxfx y dx 11(,)kkxkkxyyfx y dx 11(,)kkkkkkyyfxyxx 2022-4-258对于初值问题对于初值问题的数值解法，首先要解决的问题就是的数值解法，首先要解决的问题就是如何如何对微分方对微分方程中的程中的导数导数进行离散化，建立求数值解的递推公式。进行离散化，建立求

8、数值解的递推公式。递推公式通常有两类，一类是计算递推公式通常有两类，一类是计算yi+1时只用到时只用到xi+1, xi 和和yi，即前一步的值，因此有了初值以后就可以逐步，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为往下计算，此类方法称为单步法单步法；其代表是；其代表是龙格龙格库塔法库塔法。另一类是计算。另一类是计算yi+1i+1时，除用到时，除用到xi+1i+1, ,xi i和和yi i以以外，还要用到外，还要用到 ，即前面，即前面k步的值，此类方法称为步的值，此类方法称为多步法多步法其代表是亚当斯法。其代表是亚当斯法。 00)(),(yxyyxfy), 2 , 1(,kp

9、yxpipi 2022-4-259 一阶微分方程的数值解可用一阶微分方程的数值解可用数值微分数值微分、数值积数值积分分和和泰勒展开泰勒展开等方法得到。等方法得到。hxyxyxykkk)()()(1 ),(yxfy 微微分分方方程程),()()(1kkkkyxfhxyxy ),(1kkkkyxhfyy 9.2 简单的数值方法与基本概念简单的数值方法与基本概念此公式称为欧拉公式此公式称为欧拉公式的的导导数数用用差差商商代代替替节节点点kx2022-4-2510 ),(1kkkkyxhfyy 欧拉公式欧拉公式 )(),(001xyyyxhfyyiiiiEuler 法的计算格式为：法的计算格式为： )

10、2 . 9(1, 1 , 0 ni )(),(00 xyyyxfy微微分分方方程程2022-4-2511 Euler 格式格式用用泰勒展开泰勒展开方法得到。方法得到。，微分方程),(yxfy 展开处在节点将Taylorxxykk)(1 )(2)()()()(21 yhxyhxyhxyxykkkk),(1kkkkyxhfyy )(22 yhh的高次项略去)()()(1kkkxyhxyxy 2022-4-2512 Euler 格式格式用用数值积分数值积分方法得到。方法得到。 11),(iiiixxxxdxyxfdxy)3 . 9()(,),()()(111 iiiixxxxiidxxyxfdxyx

11、fxyxy 选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项就会得到不同的计算公式。就会得到不同的计算公式。 1)(,iixxdxxyxf),(yxfy 1, iixx将方程将方程 的两端在区间的两端在区间 上积分得，上积分得，2022-4-2513 用左矩形方法计算积分项用左矩形方法计算积分项 )(,)()(,11iiiixxxyxfxxdxxyxfii 代入代入(9.3)(9.3)式式, ,并用并用 yi i近似代替式中近似代替式中 y( (xi i) )即可得即可得到向前欧拉（到向前欧拉（Euler）公式）公式 ),(1iiiiyxhfyy 由于数值积分的矩形方法精度

12、很低，所以由于数值积分的矩形方法精度很低，所以欧拉（欧拉（Euler）公式当然很粗糙。）公式当然很粗糙。 )3 . 9()(,)()(11 iixxiidxxyxfxyxy2022-4-2514Euler 公式的几何解释公式的几何解释 欧拉（欧拉（Euler）方法是解初值问题的最简单）方法是解初值问题的最简单的数值方法。初值问题的数值方法。初值问题的解的解 y=y(x)为通过点为通过点 的一条积分曲线。的一条积分曲线。积分曲线上每一点积分曲线上每一点 的切线的斜率的切线的斜率 等等于函数于函数 在这点的值。在这点的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf2

13、022-4-2515 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) 在在P0点上点上切线切线 ( (其斜率为其斜率为 ),),与与x=x1 1直线直线10PP),()(000yxfxy 相交于相交于P1点点( (即点即点( (x1 1, ,y1 1),),得到得到y1 1作为作为y( (x1 1) )的近似值的近似值, ,如上图所示。过点如上图所示。过点( (x0 0, ,y0

14、0),),以以f( (x0 0, ,y0 0) )为为斜率的切线斜率的切线方程为方程为 当当 时时, ,得得 )(,(0000 xxyxfyy 1xx )(,(010001xxyxfyy 这样就获得了这样就获得了P1 1点的坐标点的坐标: : hyxfyy),(0001 2022-4-2516 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 同样同样, 过过点点P1(x1 1, ,y1 1),),作积分曲线作积分曲线 y= =y( (x) )的切线的切线交直线交直线 x=x2 2于于P2点点, ,切线切线 的斜率的斜率直线方程为直线

15、方程为21PP),()(111yxfxy )(,(1111xxyxfyy )(,(121112xxyxfyy 当当 时时, ,得得 2xx ),(1112yxhfyy 2022-4-2517当当 时时, ,得得 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程, ,就可获得一系就可获得一系列的点列的点: : 。对已求得点对已求得点以以 为斜率作直线为斜率作直线 ),(kkkyxP)(,(kkkkxxyxfyy 1 kxx)(,(11kkkkkkxxyxfyy .

16、 1, 1 , 0,)(11 nkyxykk取取nPPP,21),()(kkkyxfxy 2022-4-2518 从图形上看从图形上看, ,就获得了一条近似于曲线就获得了一条近似于曲线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,21nPPPP3212022-4-2519通常取通常取 ( (常数常数),),则则hhxxiii 1)2 . 9(1, 1 , 0)(),(001 nixyyyxh

17、fyyiiii 00)(),(yxyyxfy微分方程的初值问题微分方程的初值问题Euler 法的计算格式为：法的计算格式为： 用折线近似于曲线得到的计算公式称为欧拉公式用折线近似于曲线得到的计算公式称为欧拉公式2022-4-2520例例9.1 用欧拉法解初值问题用欧拉法解初值问题 )6 . 00( ,1)0(2 xyxyyy取步长取步长 h=0.2 ,=0.2 ,计算过程保留计算过程保留4 4位小数位小数 解解: :)2 , 1 , 0(1)4(2 . 001 kyyxyykkkk欧拉迭代格式欧拉迭代格式),(1kkkkyxhfyy 的的近近似似值值计计算算)6 . 0(),4 . 0(),2

18、 . 0(yyy)( 2 . 02kkkkyxyy ihxxhi020,.232106040200 xyyyxfxxxx),(,.,.,.,2022-4-2521例例9.1 用欧拉法解初值问题用欧拉法解初值问题 )6 . 00( ,1)0(2 xyxyyy取步长取步长 h=0.2 ,=0.2 ,计算过程保留计算过程保留4 4位小数。位小数。 解解: :当当 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.2

19、0.8)0.6144)2 , 1 , 0(1)4(2 . 001 kyyxyykkkk当当 k=2, x3 =0.6时，已知时，已知x2 =0.4， y2 =0.6144，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 2022-4-25221(,)nnnnyyhfxy clear; 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.358

20、2 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公公式式求求解解利利用用取取例例Eulerh, 1 . 02 )10( ,1)0(2 xyyxyy. . 20.1()0.21.1nnnnnnnxyyyxyy 解：计算结果如下：2022-4-2523对方程对方程 的两端在区间的两端在区间 上积分得，上积分得，),(yxfy 1, iixx) 2 . 9(),(iiiyxh

21、fy 2)(,()(,()(,1111 iiiiiixxxyxfxyxfxxdxxyxfii代入代入(9.(9.2 2) )式式, ,即可得到梯形公式即可得到梯形公式 为了提高精度为了提高精度,改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 9.2.2 梯形公式梯形公式 由于由于欧拉（欧拉（Euler）公式使用的是）公式使用的是数值积分的左数值积分的左矩形方法精度很低，所以欧拉公式很粗糙。矩形方法精度很低，所以欧拉公式很粗糙。 1)(,)()(1iixxiidxxyxfxyxy2022-4-2524)2 . 9()(,)()(11 iixxiidxxyxfxyxy111(, ()(,

22、 () , ( )2iixiiiixf x y xf xy xf x y x dxh )5 . 9(),(),(2111 iiiiiiyxfyxfhyy 由于数值积分的梯形公式比矩形公式的精度高，由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（因此梯形公式（9.59.5）是比欧拉公式）是比欧拉公式( 9.2 )( 9.2 )精度高精度高的一个数值方法。的一个数值方法。 9.2.2 梯形公式梯形公式梯形公式梯形公式2022-4-2525111( ,)(,)2iiiiiihyyf x yf xy ( 9.5 ) ( (9.2) 和和(9.5)式粗看差不多式粗看差不多， )(),(001xy

23、yyxhfyyiiii( 9.2 ) Euler 公式梯形公式但但( (9.5)式的右端也含有未知的式的右端也含有未知的 yi+1, ,它是一个关于它是一个关于 yi+1的方程的方程, ,这类数值方法称为这类数值方法称为隐式方法隐式方法。相反地。相反地, ,欧欧拉公式是关于拉公式是关于 yi+1 的一个直接的计算公式，这类数的一个直接的计算公式，这类数值方法称为值方法称为显式方法显式方法。隐式公式隐式公式,需用迭代法求解。需用迭代法求解。2022-4-25269.2.3 两步欧拉公式两步欧拉公式 对方程对方程 的两端在区间上的两端在区间上 积分得积分得 ),(yxfy 11, iixx 11)

24、(,)()(11iixxiidxxyxfxyxy ( 9.6 ) 改用改用中矩形中矩形公式计算其积分项，即公式计算其积分项，即 )(,)()(,1111iiiixxxyxfxxdxxyxfii 代入上式代入上式, ,并用并用 yi 近似代替式中近似代替式中y(xi)即可得到两步的即可得到两步的欧拉公式欧拉公式 )7 .9(),(211iiiiyxhfyy 111111,),(),( iiiiiiiiiyxPyxPyxP1 iy2022-4-2527 前面介绍的欧拉方法、梯形方法，它们都前面介绍的欧拉方法、梯形方法，它们都是单步法是单步法, ,其特点是在计算其特点是在计算 yi+1i+1时只用到

25、前一步的时只用到前一步的信息信息 yi i; ; 无论是单步法无论是单步法, ,还是两步法只要是直接法就还是两步法只要是直接法就需估计它的误差。需估计它的误差。 而中矩形公式而中矩形公式( (9.7)中除了中除了yi i外外, ,还用到更前一还用到更前一步的信息步的信息 yi-1i-1, ,即调用了前两步的信息即调用了前两步的信息, ,故称其为故称其为两步欧拉公式。两步欧拉公式。 2022-4-2528 定义定义9.1 在在 yi准确准确的前提下的前提下, 即即 时时, 用数值方用数值方法计算法计算yi+1的误差的误差 , 称为该数值方法计称为该数值方法计算时算时 yi+1 的局部截断误差。的

26、局部截断误差。)(iixyy 111)( iiiyxyR9.2.4 欧拉法的局部截断误差欧拉法的局部截断误差 衡量求解公式好坏的一个主要标准是求解公式的衡量求解公式好坏的一个主要标准是求解公式的精度精度, 因此引入局部截断误差和收敛阶的概念。因此引入局部截断误差和收敛阶的概念。()iiiey xy 整整体体截截断断误误差差。 在假设前在假设前 i-1 步精确的前提下，得到的第步精确的前提下，得到的第 i 步的步的误差称为局部截断误差误差称为局部截断误差处的精确值在iixyxy)(处的近似值在iixyy2022-4-2529假定假定 ，由欧拉公式则有由欧拉公式则有)(iixyy )()(1iii

27、xyhxyy ),()(! 2)()()(121 iiiiixxyhxyhxyxy)(! 2)(211 yhyxyii因此有因此有 展展开开处处二二阶阶在在将将精精确确解解Taylorxhxyxyiii)()(1 近似解局部截断误差局部截断误差 1ke)(2hO 2022-4-2530定义定义9.2 如果数值方法的局部截断误差为如果数值方法的局部截断误差为 ,则称这种数值方法的则称这种数值方法的阶数是阶数是 P。)(1 phO 显然，步长显然，步长(h N 结束。结束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxx2022

28、-4-2539（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 2022-4-2540例例9.2 9.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 区间为区间为 0,10,1 , ,取步长取步长h=0.1h=0.1 解解: : 改进欧拉法的具体形式改进欧拉法的具体形式 本题的精确解为本题的精确解为xxy21)( 0,1,9i 0

29、00,1xy ,)(21)2( 1 . 0)2( 1 . 011 cpipipiciiiipyyyyxyyyyxyyy1)0(,2 yyxyy),(1iiiiyxhfyyEuler 公公式式2022-4-2541clearx=0,yn=1 %初始化for n=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 %校正end2022-4-2542例例9.3 对初值问题对初值问题 1)0(0yyy求得的近似解为求得的近似解为 nnhhy 22并证明当步长并证明当步长h h0 0时时, ,yn n收敛于精

30、确解收敛于精确解xe ),(),(2111 nnnnnnyxfyxfhyyyyxf ),(211 nnnnyyhyy整理成显式整理成显式 nnyhhy 221反复迭代反复迭代, ,得到得到 证明证明: : 解初值问题的梯形公式为解初值问题的梯形公式为证明用梯形公式证明用梯形公式2022-4-2543nnyhhy 221反复迭代反复迭代, ,得到得到 nnyhhy22110 ynnhhy 2201231222.2222yhhyhhyhhnnn 由于由于 ，有，有 nhx hxhnhhhy 22limlim00 xnhy e0limhxhhhhhhxhhhhh 222200221221limlim

31、x e证毕证毕 2022-4-25449.3 龙格龙格-库塔（库塔（Runge-Kutta）法）法)()(1iixyxy 只有对平均高度只有对平均高度 f(x,y) 提供一种算法，便可得到提供一种算法，便可得到一种计算格式。一种计算格式。hhxyhxfhhxyiii)(),()( )(,1 yhfyyii9.3.1 龙格龙格-库塔法的基本思想库塔法的基本思想上的平均斜率区间,)(,1 iixxyf),(yxfy 11),(iiiixxxxdxyxfdxy上的平均高度区间或,1 iixx2022-4-2545 111),(hkyyyxfkiiii )2(),(),(2111121kkhyyhky

32、xfkyxfkiiiiiiEuler 公式可改写成公式可改写成改进的改进的Euler公式又可改写成公式又可改写成用左端点的高度作为平均高度得到的计算格式用左端点的高度作为平均高度得到的计算格式用左、右端点高度的平均值作为平均高度得到的计用左、右端点高度的平均值作为平均高度得到的计算格式算格式2022-4-2546 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点: :都是用都是用f(x,y)在某些点上值在某些点上值( (高度高度) )的线性组合得出的线性组合得出 y(xi+1i+1) )的的近似值近似值 yi+1i+1, ,而且增加计算而且增加计算 f(x,y)的次数的次数,

33、,可提高截可提高截断误差的阶。如欧拉公式断误差的阶。如欧拉公式: :每步计算一次每步计算一次 f(x,y)的值的值, ,为一阶方法为一阶方法局部截断误差为局部截断误差为 。改进欧拉公式。改进欧拉公式需计算两次需计算两次 f(x,y) 的值，它是二阶方法。它的局部的值，它是二阶方法。它的局部截断误差为截断误差为 。)(3hO2( )Oh2022-4-2547 于是可考虑在于是可考虑在 内多内多预报几个点的函数预报几个点的函数值值( (高度高度)，然后将其加权平均作为平均，然后将其加权平均作为平均高度高度，构造时，构造时要求近似公式在要求近似公式在( (xi,yi) )处的处的Taylor展开式与

34、解展开式与解y(x)在在 xi i 处的处的Taylor展开式的前面几项重合，从而使近似公式展开式的前面几项重合，从而使近似公式达到所需要的阶数。则可构造出更高精度的计算格式，达到所需要的阶数。则可构造出更高精度的计算格式，这就是龙格这就是龙格库塔（库塔（Runge-Kutta）法的基本思想。）法的基本思想。 1, iixx2022-4-2548只有多计算几个点的函数值，平均高度就更精确只有多计算几个点的函数值，平均高度就更精确.龙格龙格-库塔法的基本思想库塔法的基本思想函数值的线性组合，可以表示为上的定积分区间 1),(,1iixxiidxyxfxx),(yxfy 1),(1iixxiidx

35、yxfyy或多预报几个点的斜率值，将其加权平均作为平或多预报几个点的斜率值，将其加权平均作为平均斜率，平均斜率就更精确。均斜率，平均斜率就更精确。2022-4-25492211kkk ),(12phkyphxfkii 率的近似值，即的加权平均作为平均斜和斜率值以两点处的和上取两点在211,kkphxxxxxipiiii 9.3.2 二阶龙格二阶龙格库塔法库塔法).,(),(),(,2211pipipiiiiiyxfkxkxyyxfkxk 处处切切线线的的斜斜率率点点为为处处切切线线的的斜斜率率点点为为式中式中:比照改进的欧拉法比照改进的欧拉法),(iiipiyxphfyy 2022-4-255

36、0)(,()( yfyK 式中式中 K K可看作是可看作是 y=y(x)在区间在区间 上的平均斜率。上的平均斜率。所以可得计算公式为：所以可得计算公式为： 1, iixxhKxyyii )(1)()(2211kkhxyi （9.14） )()()(11iiiixxyxyxy 也即也即 hKxyxyii )()(1（9.13）),()(12phkyphxfphxykiii 其中：其中： ),()(1iiiyxfxyk 由微分中值定理，存在点由微分中值定理，存在点 ，使得，使得),(1 iixx2022-4-2551二元函数的二元函数的TaylorTaylor展开式展开式 ),()(),(),(0

37、00000yxfykxhyxfkyhxf ),()(! 21002yxfykxh).,()()!1(1001kyhxfykxhnRnn )(),(),(),(),(200000000hOyxkfyxhfyxfkyhxfyx ),(yxfy yyxfyxfyyx ),(),(),(),(),(yxfyxfyxfyx 2022-4-2552将以上结果代入将以上结果代入)14. 9()()(22111kkhxyyii 2112()()()()()iiiiiyy xhy xy xph yxO h )16. 9()()()()()(32221hOxyphxyhxyiii 展展开开处处在在将将二二元元函函

38、数数TaylorPphkyphxfkiii),(12 将将 在在 x= =xi i 处进行二阶处进行二阶Taylor展开：展开： )()(! 2)()()(321hOxyhxyhxyxyiiii （9 9.15） )(1 ixy)(),(),(),(),(22hOyxfyxfyxfphyxfkiiyiiiixii )()()(2hOxyphxyii 2022-4-2553)16. 9()()()()()(322211hOxyphxyhxyyiiii 对式对式(9.15)(9.15)和和(9.16)(9.16)进行比较系数后可知进行比较系数后可知, ,只要只要 )17.9(211221 p成立成

39、立, ,格式格式(9.14)(9.14)的局部截断误差就等于的局部截断误差就等于)(3hO)15. 9()()(! 2)()()(321hOxyhxyhxyxyiiii 式式(9.17)(9.17)中具有三个未知量中具有三个未知量, ,但只有两个方程但只有两个方程, ,因而因而有无穷多解。有无穷多解。111iikyxye)(2022-4-2554式式(9.17)(9.17)中，若取中，若取 , ,则则P =1=1，这是无穷多，这是无穷多解中的一个解，将以上所解的值代入式解中的一个解，将以上所解的值代入式(9.14)(9.14)并改并改写可得写可得 2121 ),(),(21121211hkyx

40、fkyxfkkkhyyiiiiii 不难发现，上面的格式就是改进的欧拉格式。不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（凡满足条件式（9.179.17）有一簇形如上式的计算格式，）有一簇形如上式的计算格式，这些格式统称为二阶龙格这些格式统称为二阶龙格库塔格式。因此改进的库塔格式。因此改进的欧拉格式是众多的二阶龙格欧拉格式是众多的二阶龙格库塔法中的一种特殊库塔法中的一种特殊格式。格式。 )17.9(211221 p 2022-4-2555若取若取 , ,则则 ，此时二阶龙格，此时二阶龙格- -库塔库塔法的计算公式为法的计算公式为 0121, 12p 1,2 , 1 , 0 ni 此计算

41、公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间 的中点。的中点。 12ix 1, iixx)2,(),(12121khyxfkyxfkiiii 21hkyyii )17.9(211221 p 2022-4-25569.3.3 三阶龙格三阶龙格- -库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 pix )1( qpqhxxiqi并用三个点并用三个点 的斜率的斜率 k1 1, ,k2 2, ,k3 3 加权平均加权平均得出平均斜率得出平均斜率 k* *的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式:

42、 : qipiixxx ,（9.18） 为了预报点为了预报点 的斜率值的斜率值k3 3, ,在区间在区间 内有两内有两个斜率值个斜率值 k1 1和和 k2 2 可以用可以用, ,可将可将k1 1, ,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率, ,从而得到从而得到 的预报值的预报值 qix,qiixx )(qixyqiy 1122i qiyyqhkk,qiixx ),(),(213322111iiiiiiiphkyphxfkyxfkkkkhyy 2022-4-2557于是可得于是可得 ),(3qiqiyxfk 运用运用Taylor展开方法选择参数展开方法选择参数 , ,可以使格

43、式可以使格式( (9.18)的局部截断误差为的局部截断误差为 , ,即具有即具有三阶精度三阶精度. .)(4hO 613121112332223232121pqqpqp21321, qp满满足足：可可得得21321, qp。截断误差为截断误差为)(4hO2022-4-2558局部截断误差为局部截断误差为 , ,即具有三阶精度，这类格即具有三阶精度，这类格式统称为式统称为三阶龙格三阶龙格库塔方法库塔方法。)(4hO)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（9.19） 时时，当当2, 1,61,64,61, 1,212

44、1321 qp是其中的一种，称为库塔（是其中的一种，称为库塔（Kutta）公式。）公式。 2022-4-25599.3.4 四阶龙格四阶龙格库塔法库塔法 如果需要再提高精度，用类似上述的处理方法，如果需要再提高精度，用类似上述的处理方法，只需在区间只需在区间 上用四个点处的斜率加权平均作为上用四个点处的斜率加权平均作为平均斜率平均斜率 k*的近似值，构成一系列四阶龙格的近似值，构成一系列四阶龙格库塔公库塔公式。具有四阶精度，即局部截断误差是式。具有四阶精度，即局部截断误差是 。 由于推导复杂，这里从略，只介绍最常用的一种由于推导复杂，这里从略，只介绍最常用的一种四阶经典四阶经典RK公式公式。

45、,1 iixx)(5hO2022-4-2560 四阶经典四阶经典 RK公式公式 )22(6),()2,()2,(),(43211314221312121kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（9.20） 2022-4-25619.3.5 9.3.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1) 计算步骤计算步骤 输入输入 , ,h,Nh,N 使用龙格使用龙格库塔公式（库塔公式（9.209.20）计算出）计算出y y1 1 输出输出 ，并使，并使 转到转到 直至直至n n N N 结束。结束。 10, xx11, yx0101,yyxx2022

46、-4-2562(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 2022-4-2563例例9.4 9.4 取步长取步长 h= =0.2，用经典，用经典 R-C 格式求解初值问

47、题格式求解初值问题 1)0(2yxyy10 x解解: : 2 . 0, 1, 0,2),(00 hyxxyyxf0),(001 yxfk2 .0)1 , 1 .0()2,(10202 fkhyxfkh204.0)02.1 , 1 .0()2,(20203 fkhyxfkh41632. 0)0408. 1 , 2 . 0(),(3004 fhkyxfkh040811. 1)22(62 . 0432101 kkkkyy),()2,()2,(),(314221312121hkyxfkkhyxfkkhyxfkyxfkiiiiiiii )22(643211kkkkhyyii 由四阶龙格由四阶龙格-库塔公

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

49、适的算法。当针对问题的具体特点选择合适的算法。 2022-4-2565 以经典的四阶龙格以经典的四阶龙格- -库塔法库塔法( (9.20)为例。从节点为例。从节点x xi i出发，出发，先以先以h为步长求出一个近似值，记为为步长求出一个近似值，记为 ，由于局部截断误，由于局部截断误差为差为 ，故有，故有 )(1hiy )(5hO5)(11)(chyxyhii 当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。9.3.6 变步长的龙格变步长的龙格-库塔法库塔法 在微分方程的数值解中，选择适当的步长是非常重要在微分方程的数值解中，选择适当的步长是非常重要

50、的。单从每一步看，的。单从每一步看，步长越小，截断误差就越小；步长越小，截断误差就越小；但随着但随着步长的缩小，在一定的求解区间内所要完成的步长的缩小，在一定的求解区间内所要完成的步数步数就就增加增加了了。这样会引起计算量的增大，并且会这样会引起计算量的增大，并且会引起舍入误差引起舍入误差的大的大量量积累与传播。积累与传播。因此微分方程数值解法也有选择步长的问因此微分方程数值解法也有选择步长的问题。题。2022-4-2566然后将步长折半然后将步长折半, ,即以为即以为 步长步长, ,从节点从节点x xi i出发出发, ,跨跨两步到节点两步到节点x xi+1i+1, ,再求得一个近似值再求得一

51、个近似值 , ,每跨一步的每跨一步的截断误差是截断误差是 , ,因此有因此有2h)(21hiy52 hc5)2(1122)()( hcxyxyhii这样这样 161)()()(11)2(11 hiihiiyxyyxy)(151)()(1)2(1)2(11hihihiiyyyxy 由此可得由此可得 这表明以这表明以 作为作为 的近似值，其误差可用步的近似值，其误差可用步长折半前后两次计算结果的偏差长折半前后两次计算结果的偏差 )2(1hiy)(1ixy)(1)2(1hihiyy 来判断所选步长是否适当来判断所选步长是否适当2022-4-2567当要求的数值精度为当要求的数值精度为时：时： （1

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

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

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

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

56、式反复递推后，可得将上式反复递推后，可得 011)1(yhyii), 2 , 1()1 (00iyyhyiii或或h1式中式中2022-4-2571 要使要使 y yi i 有界，其充要条件为有界，其充要条件为 1 11 h即即 由于由于 ，故有，故有 0 20h（9.269.26） 可见，如欲保证算法的稳定，显式可见，如欲保证算法的稳定，显式 EulerEuler 格格式的步长式的步长 h h 的选取要受到式（的选取要受到式（9.269.26）的限制。）的限制。 的绝对值越大，则限制的的绝对值越大，则限制的 h h 值就越小。值就越小。 用隐式用隐式 EulerEuler 格式格式右矩形公式

57、，对模型方程的计，对模型方程的计算公式为算公式为 )(11 iiiyhyy iiyhy 111可化为可化为00)1(yyhyiii 2022-4-2572由于由于 , ,则恒有则恒有 ，故恒有，故恒有 0111 hiiyy 1 因此，隐式因此，隐式 Euler 格式是绝对稳定的（无条件稳格式是绝对稳定的（无条件稳定）的（对任何定）的（对任何 h h 0 0）。）。 2022-4-2573 常微分方程初值问题的求解，是将微分方程转常微分方程初值问题的求解，是将微分方程转化为差分方程来求解，并用计算值化为差分方程来求解，并用计算值y yi i来近似替代来近似替代y(xy(xi i) )，这种近似替

58、代是否合理，还须看分割区间，这种近似替代是否合理，还须看分割区间 的长度的长度 h h 越来越小时，即越来越小时，即 时，时， 是否成立。若成立，则称该方法是收敛的，否则称是否成立。若成立，则称该方法是收敛的，否则称为不收敛。为不收敛。 iixx,1 01 iixxh)(iixyy ),(1iiiiyxhfyy 9.4.2 收敛性收敛性 这里仍以这里仍以Euler方法为例，来分析其收敛性。方法为例，来分析其收敛性。Euler格式格式2022-4-2574设设 表示取表示取 时时, 按按Euler公式的计算结果公式的计算结果, 即即1iy)(iixyy )(,()(1iiiixyxhfxyyEu

59、ler方法局部截断误差为：方法局部截断误差为：)()(2)(1211 iiiixxyhyxy设有常数设有常数 ，则，则 )(max21xycbxa 211)(chyxyii（9.27） 总体截断误差总体截断误差 1111111)()(iiiiiiiyyyxyyxy)(,()(),(11iiiiiiiixyxhfxyyxhfyyy),()(,()(iiiiiiyxfxyxfhyxy（9.28） 又又 由于由于f(x,y)f(x,y)关于关于y y满足李普希茨条件满足李普希茨条件, ,即即 2022-4-2575iiiiiiyxyLyxfxyxf )(),()(代入代入(9.28)上式，有上式，有

60、 iiiiyxyhLyy )()1(11ihL )1(再利用式（再利用式（9.27），式（），式（9.28） 1111111)()( iiiiiiiyyyxyyxy2)1 (chhLi21)1 (chhLii即即 上式反复递推后，可得上式反复递推后，可得 1020)1 ()1 (ikkiihLchhL1)1 ()1 (0iihLLchhL（9.29） 2022-4-2576设设 （T T为常数）为常数） Tihxxi0hLehL 1TLihLieehL)1 (因为因为 所以所以 把上式代入式（把上式代入式（9.299.29），得），得 ) 1(0TLTLieLche若不计初值误差，即若不计初值