第3章_常微分方程的差分方法

1、数值计算方法数值计算方法 （数值分析）（数值分析） 第三章：第三章： 常微分方程的差分方法常微分方程的差分方法 包含自变量、未知函数及未知函数的导数或微 分的方程称为微分方程。在微分方程中, 自变量的 个数只有一个, 称为常微分方程.。自变量的个数 为两个或两个以上的微分方程叫偏微分方程。微分 方程中出现的未知函数最高阶导数的阶数称为微分 方程的阶数。如果未知函数y及其各阶导数 都是一次的,则称它是线性的,否则称为非线性的。 )( , n yyy 引言 在高等数学中，对于常微分方程的求解，给出 了一些典型方程求解析解的基本方法，如可分离变 量法、常系数齐次线性方程的解法、常系数非齐次 线性方程

2、的解法等。但能求解的常微分方程仍然是 有限的，大多数的常微分方程是不可能给出解析解 。 譬如 22 yxy 这个一阶微分方程就不能用初等函数及其积 分来表达它的解。 再如，方程 1)0(y yy 的解 ,虽然有表可查,但对于表上 没有给出 的值,仍需插值方法来计算 x ey x e 从实际问题当中归纳出来的微分方程，通常主要依 靠数值解法来解决。本章主要讨论一阶常微分方程 初值问题 00 )( ),( yxy yxfy 在区间a x b上的数值解法。 可以证明,如果函数在带形区域 R=axb, -y内连续，且关于y满足李普希兹(Lipschitz) 条件，即存在常数L(它与x,y无关)使 21

3、21 ),(),(yyLyxfyxf 对R内任意两个 都成立,则方程的解 在a, b上存在且唯一。 21, y y )(xyy 数值方法的基本思想数值方法的基本思想 对常微分方程初值问题数值解法，就是要算出 精确解y(x)在区间a,b上的一系列离散节点处的函 数值 的近似值 。相邻两个节点的间距 称为步长，步 长可以相等，也可以不等。本章总是假定h为定数， 称为定步长，这时节点可表示为 数值解法需要把连续性的问题加以离散化，从而求 出离散节点的数值解。 bxxxxa nn 110 )(,),(),( 10n xyxyxy n yyy, 10 ii xxh 1 niihxxi, 2 , 1 ,

4、0 对常微分方程数值解法的基本出发点就是离散 化。其数值解法有两个基本特点，它们都采用“步 进式”，即求解过程顺着节点排列的次序一步一步 地向前推进，描述这类算法，要求给出用已知信息 计算 的递推公式。建立 这类递推公式的基本方法是在这些节点上用数值积 分、数值微分、泰勒展开等离散化方法，对初值问 题 中的导数 进行不同的离散化处理。 021 ,yyyy iii 1i y y 00 )( ),( yxy yxfy 对于初值问题 的数值解法，首先要解决的问题就是如何对微分方程 进行离散化，建立求数值解的递推公式。递推公式通 常有两类，一类是计算yi+1时只用到xi+1, xi 和yi，即 前一步

5、的值，因此有了初值以后就可以逐步往下计算 ，此类方法称为单步法；其代表是龙格库塔法。另 一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用 到 ，即前面k步的值，此类方法称为多步法；其代表 是亚当姆斯法。 00 )( ),( yxy yxfy ), 2 , 1( ,kpyx pipi EulerEuler公式公式 欧拉（Euler）方法是解初值问题的最简单的 数值方法。初值问题 的解y=y(x)代表通过点 的一条称之为微分方 程的积分曲线。积分曲线上每一点 的切线的斜率 等于函数 在这点的 值。 00 )( ),( yxy yxfy ),( 00 yx ),(yx )(x y ),(

6、yxf 简单数值方法 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 EulerEuler法的求解过程是法的求解过程是: :从初从初 始点始点P P0 0( (即点即点(x(x0 0,y,y0 0)出发出发, , 作积分曲线作积分曲线y=y(x)y=y(x)在在P P0 0点上点上 切线切线 ( (其斜率为其斜率为 ),),与与x=xx=x1 1直线直线 10P P ),()( 000 yxfxy 相交于相交于P P1 1点点( (即点即点(x(x1 1,y,y1 1),),得到得到y y1 1作为作为y(xy(x1 1) )的近似

7、值的近似值, , 如上图所示。过点如上图所示。过点(x(x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0) )为斜率的切线为斜率的切线 方程为方程为 当当 时时, ,得得 )(,( 0000 xxyxfyy 1 xx )(,( 010001 xxyxfyy 这样就获得了这样就获得了P P1 1点的坐标。点的坐标。 Pi+1 Pn y=y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 同样同样, , 过点过点P P1 1(x(x1 1,y,y1 1),),作积分曲线作积分曲线y=y(x)y=y(x)的切线的切线 交直线交直线x=xx=x2

8、 2于于P P2 2点点, ,切线切线 的斜率的斜率 = = 直线方程为直线方程为 21P P )( 1 x y ),( 11 yxf )(,( 1111 xxyxfyy )(,( 121112 xxyxfyy 当当 时时, ,得得 2 xx 当 时,得 Pi+1 Pn y= y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此获得了P2的坐标。重复以上过程,就可获得一系 列的点:P1,P1,Pn。对已求得点 以 = 为斜率作直线 ),( nnn yxP )( n x y ),( nn yxf )(,( nnnn xxyxfyy 1 n xx)(,(

9、11nxnnnn xxyxfyy nn yxy)(取 从图形上看从图形上看, ,就获得了一条近似于曲线就获得了一条近似于曲线 y=y(x)y=y(x) 的折线的折线 。 Pi+ 1 P n y= y(x) P 1 Pi Pn Pi+ 1 P 0 x0 x1 xi xi+ 1 xn P i P 1 这样这样, ,从从x x0 0逐个算出逐个算出 对应的数值解对应的数值解 n xxx, 21 n yyy, 21 n PPPP 321 通常取通常取 ( (常数常数),),则则EulerEuler法的计算格法的计算格 式式 hhxx iii 1 )( ),( 00 1 xyy yxhfyy iiii

10、i i=0,1,=0,1,n n ( ( 3 3.2 ) .2 ) 还可用还可用数值微分数值微分、数值积分法数值积分法和和泰勒展开法泰勒展开法推导推导 EulerEuler格式。以数值积分为例进行推导。格式。以数值积分为例进行推导。 将方程将方程 的两端在区间的两端在区间 上积分得，上积分得， ),(yxfy 1 , ii xx 11 ),( i i i i x x x x dxyxfdxy 11 )(,)(),()()( 1 i i i i x x i x x ii dxxyxfxydxyxfxyxy 选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项 , ,就会得到不同的

11、计算公式。就会得到不同的计算公式。 1 )(, i i x x dxxyxf (3.3) 用左矩形方法计算积分项用左矩形方法计算积分项 )(,)()(, 1 1 iiii x x xyxfxxdxxyxf i i 代入代入( (3 3.3).3)式式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到即可得到 向前欧拉（向前欧拉（EulerEuler）公式）公式 ),( 1iiii yxhfyy 由于数值积分的矩形方法精度很低，所以由于数值积分的矩形方法精度很低，所以 欧拉（欧拉（EulerEuler）公式当然很粗糙。）公式当然很粗糙。 例例1 用欧拉法解初值问题

12、用欧拉法解初值问题 1)0( )6 . 00( 2 y xxyyy 取步长取步长h=0.2 ,计算过程保留计算过程保留4位小数位小数 解解: h=0.2, 欧拉迭代格式欧拉迭代格式 2 ),(xyyyxf 2 1 ),( iiiiiiii yhxhyyyxhfyy )2 , 1 , 0()4(2 . 0iyxy iii 当当 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

13、=2, x3 =0.6时，已知时，已知x2 =0.4， y2 =0.6144，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 0.1, 2 , 01) (0)1, hEuler x yyx y y 取利用公式求解 （ 例例 . . 2 . 0 1.1) 2 ( 1 计算及结果如下解： n n n n n nnn y x y y x yhyy clear; y=1, x=0, %初始化 for n=1:10 y=1.1*y-0.2*x/y, x=x+0.1, end y = 1 x = 0 y = 1.1000 x = 0.1000 y = 1.1918

14、x = 0.2000 y = 1.2774 x = 0.3000 y = 1.3582 x = 0.4000 y = 1.4351 x = 0.5000 y = 1.5090 x = 0.6000 y = 1.5803 x = 0.7000 y = 1.6498 x = 0.8000 y = 1.7178 x = 0.9000 y = 1.7848 x = 1.0000 梯形公式梯形公式 为了提高精度,对方程 的两端在区间上 积分得， 改用梯形方法计算其积分项，即 ),(yxfy 1 , ii xx 1 )(,)()( 1 i i x x ii dxxyxfxyxy )(,()(,( 2 )(

15、, 11 1 1 iiii ii x x xyxfxyxf xx dxxyxf i i (3.4 ) 代入(3.4)式,并用近似代替式中即可得到梯形公式 ),(),( 2 111 iiiiii yxfyxf h yy ( 3.5 ) 梯形公式（3.5）比欧拉公式( 3.2 )的精度高 ),(),( 2 111 iiiiii yxfyxf h yy ( 3.5 ) ( (3 3.5).5)式的右端含有未知的式的右端含有未知的y yi+1 i+1, ,它是一个关于 它是一个关于 y yi+1 i+1的函数方程 的函数方程, ,这类数值方法称为这类数值方法称为隐式方法隐式方法。相。相 反地反地, ,

16、欧拉法是关于欧拉法是关于y yi+1 i+1的一个直接的计算公式， 的一个直接的计算公式， 这类数值方法称为这类数值方法称为显式方法显式方法。 两步欧拉公式两步欧拉公式 对方程 的两端在区间上 积分得 ),(yxfy 11,ii xx 1 1 )(,)()( 11 i i x x ii dxxyxfxyxy ( 3.6 ) 改用中矩形公式计算其积分项，即 )(,)(, 11 1 1 iiii x x xyxfxxdxxyxf i i 代入上式,并用yi近似代替式中y(xi)即可得到两步 欧拉公式 ),(2 11iiii yxhfyy ( 3.7 ) 前面介绍过的数值方法前面介绍过的数值方法,

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

18、下, 即 时, 用数值 方法计算yi+1的误差 , 称为该数值方法 计算时yi+1的局部截断误差。 对于欧拉公式，假定 ，则有 )( ii xyy 11) ( iii yxyR )( ii xyy )()()(,()( 1iiiiii xyhxyxyxfhxyy 而将真解y(x)在xi处按二阶泰勒展开 ),()( ! 2 )()()( 1 2 1 iiiii xxy h xyhxyxy )( !2 )( 2 11 y h yxy ii 因此有 定义 2 数值方法的局部截断误差为 ,则称这 种数值方法的阶数是P。步长(h N n N 结束。结束。 10 ,xx )( 2 1 ),( ),( 1

19、1 cpi piic iiip yyy yxhfyy yxhfyy 11, y x 0101 ,yyxx ( (3)3) 程序实现程序实现( (改进欧拉法计算常微改进欧拉法计算常微 分方程初值问题分方程初值问题 ) 例例9.2 9.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 1)0( 2 y y x yy 区间为区间为 0,10,1 , ,取步长取步长h=0.1h=0.1 解解: : 改进欧拉法的具体形式改进欧拉法的具体形式 )( 2 1 ) 2 (1.0 ) 2 (1.0 1 1 cpi p i pic i i iip yyy y x yyy y x yyy 本题的精确解为本题的精确解

20、为 , , xxy21)( clear x=0,yn=1 %初始化 for n=1:10 yp=yn+0.1*(yn-2*x/yn); %预测 x=x+0.1; yc=yn+0.1*(yp-2*x/yp) ; yn=(yp+yc)/2 %校正 end 例例 3 对初值问题对初值问题 1)0( 0 y yy 证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为 n n h h y nhx 并证明当步长并证明当步长h0时时,yn收敛于精确解收敛于精确解 证明证明: 解初值问题的梯形公式为解初值问题的梯形公式为 x e ),(),( nnnnnn yxfyxf h yy yyxf),( 2 11

21、 nnnn yy h yy 整理成显式整理成显式 nn y h h y 反复迭代反复迭代,得到得到 y h h y h h y h h y h h y n nnnn . 1 0 y n n h h y 2 2 由于由于 ，有，有 nhx x x x x h x h h h x h n h h h h h y e e e 2 1 2 1 lim 2 2 limlim 2 2 2 2 2 2 000 n n h h y 2 2 x n h y elim 0 证毕证毕 作业题！ 1、列出下列初值问题的欧拉格式：、列出下列初值问题的欧拉格式： （1） 0=x=0.4, y(0)=1, 取取h=0.2

22、（2） 1=x=1.2 y(1)=1, 取取h=0.1 2、用梯形方法做、用梯形方法做 1=x N N 结束。结束。 10 ,xx 11 , yx 0101 ,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 k3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 输 出 x1, y1 n + 1 n n =

23、N ? x1 x0 y1 y0 结 束 n y 程序实现程序实现( (四阶龙格四阶龙格- -库塔法计库塔法计 算常微分方程初值问题算常微分方程初值问题) ) 例例4 4 取步长取步长h=0.2h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 1)0( 2 y xyy 10 x 解解: : 由四阶龙格由四阶龙格- -库塔公式可得库塔公式可得 2 . 0, 1, 0,2),( 00 hyxxyyxf 0),( 001 yxfk 2 . 0) 1 , 1 . 0() 2 ,( 10 2 0 2 fk h yxfk h 204. 0)02. 1 , 1 . 0() 2 ,( 20 2 0 3

24、 fk h yxfk h 41632. 0)0408. 1 , 2 . 0(),( 3004 fhkyxfk h 040811. 1)22( 6 2 . 0 432101 kkkkyy 可同样进行其余可同样进行其余y yi i的计算。本例方程的解为的计算。本例方程的解为 ，从表中看到所求的数值解具有，从表中看到所求的数值解具有4位有效数字。位有效数字。 2 x ey 龙格龙格库塔方法的库塔方法的推导基于推导基于Taylor展开方法展开方法 ，因而它要求所求的解具有较好的光滑性。如果解，因而它要求所求的解具有较好的光滑性。如果解 的光滑性差，那么，使用四阶龙格的光滑性差，那么，使用四阶龙格库塔方

25、法求得库塔方法求得 的数值解，其精度可能反而不如改进的欧拉方法的数值解，其精度可能反而不如改进的欧拉方法。 在实际计算时，应当针对问题的具体特点选择合适在实际计算时，应当针对问题的具体特点选择合适 的算法。的算法。 变步长的龙格变步长的龙格-库塔法库塔法 在微分方程的数值解中，选择适当的步长是非常在微分方程的数值解中，选择适当的步长是非常 重要的。单从每一步看，重要的。单从每一步看，步长越小，截断误差就越步长越小，截断误差就越 小；但随着步长的缩小，在一定的求解区间内所要小；但随着步长的缩小，在一定的求解区间内所要 完成的步数就增加了。这样会引起计算量的增大，完成的步数就增加了。这样会引起计算

26、量的增大， 并且会引起舍入误差的大量积累与传播。因此微分并且会引起舍入误差的大量积累与传播。因此微分 方程数值解法也有选择步长的问题。方程数值解法也有选择步长的问题。 以经典的四阶龙格以经典的四阶龙格- -库塔法库塔法( (3 3.20).20)为例。从节点为例。从节点 x xi i出发，先以出发，先以h h为步长求出一个近似值，记为为步长求出一个近似值，记为 ， 由于局部截断误差为由于局部截断误差为 ，故有，故有 )( 1 h i y )( 5 hO 5)( 11 )(chyxy h ii 当当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。 然后

27、将步长折半然后将步长折半, ,即以为即以为 步长步长, ,从节点从节点x xi i出发出发, ,跨跨 两步到节点两步到节点x xi+1 i+1, ,再求得一个近似值 再求得一个近似值 , ,每跨一步的每跨一步的 截断误差是截断误差是 , ,因此有因此有 2 h ) 2 ( 1 h i y 5 2 h c 5 ) 2 ( 11 2 2)()( h cxyxy h ii 这样这样 16 1 )( )( )( 11 ) 2 ( 11 h ii h ii yxy yxy )( 15 1 )( )( 1 ) 2 ( 1 ) 2 ( 11 h i h i h ii yyyxy 由此可得由此可得 这表明以这

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

29、方法称为这种通过步长加倍或折半来处理步长的方法称为 变步长法。变步长法。表面上看，为了选择步长，每一步都要表面上看，为了选择步长，每一步都要 反复判断反复判断，增加了计算工作量，但在方程的解，增加了计算工作量，但在方程的解 y(x)y(x)变化剧烈的情况下，总的计算工作量得到减少变化剧烈的情况下，总的计算工作量得到减少 ，结果还是合算的。，结果还是合算的。 1i y 1i y 作业题！ 3、取、取h=0.1用欧拉方法，用欧拉方法，取取h=0.2用改进欧拉方法，取用改进欧拉方法，取 h=0.4用四阶龙格用四阶龙格-库塔方法求解初值问题：库塔方法求解初值问题： 1) 0( ) 20(1 y xxy

30、y 4 4、求解微分方程初值问题：求解微分方程初值问题： （1 1）取）取h=0.1h=0.1，用欧拉方法求解。，用欧拉方法求解。( (5 5分分) ) （2 2）取）取h=0.4h=0.4，用经典四阶龙格，用经典四阶龙格- -库塔方法求解。库塔方法求解。( (6 6分分) ) （3 3）精确解）精确解y y（0.40.4）=1.=1.5127751277，分别确定上述两种方法计，分别确定上述两种方法计 算结果有几位有效数字。算结果有几位有效数字。( (4 4分分) ) 5 . 0)0( 4 . 00122 y xxyy 亚当姆斯格式亚当姆斯格式 龙格龙格-库塔库塔方法是一类重要算法，但这类算

31、法在方法是一类重要算法，但这类算法在 每一步都需要先预报几个点上的斜率值每一步都需要先预报几个点上的斜率值，计算量比，计算量比 较大。较大。 考虑到计算考虑到计算yi+1之前已得出一系列节点上之前已得出一系列节点上 的斜率值，能否利用这些已知值的斜率值，能否利用这些已知值 来减少计算量呢？来减少计算量呢？ 这就是亚当姆斯（这就是亚当姆斯（Adams）方法的设计思想。）方法的设计思想。 11 ,xxx ii 亚当姆斯方法 设用设用x xi i,x,xi i- -1 1两点的斜率值加权平均作为区间两点的斜率值加权平均作为区间 上的平均斜率，有计算格式上的平均斜率，有计算格式 1 , ii xx )

32、,( ),( )1 ( 111 11 iii iii iiii yxfy yxfy yyhyy （3 3.21.21） 选取参数选取参数，使格式（，使格式（3 3.21.21）具有二阶精度）具有二阶精度 。 亚当姆斯方法 将将 在在x xi i处处TaylorTaylor展开展开 1 i y )()( ! 2 1 )( 32 1 hOhyhyyy iii 代入计算格式代入计算格式( (3 3.21).21)化简化简, ,并假设并假设, , 因此有因此有 )(),( 11 iiii xyyxyy )()()()( 32 1 hOxyhxyhxyy iiii 与与y(xy(xi+1 i+1) )在

33、 在x xi i处的处的TaylorTaylor展开式展开式 )()( 2 1 )()()( 32 1 hOxyhxyhxyxy iiii 相比较相比较, , 需取需取 2 1 才使格式才使格式( (3 3.21).21)具有二阶精度。这样具有二阶精度。这样导出的计算格式导出的计算格式 )3( 2 11 iiii yy h yy 称之为二阶亚当姆斯格式。类似地可以导出三阶亚称之为二阶亚当姆斯格式。类似地可以导出三阶亚 当姆斯格式当姆斯格式。 )51623( 12 211 iiiii yyy h yy 和四阶亚当姆斯格式。和四阶亚当姆斯格式。 )9375955( 24 3211 iiiiii y

34、yyy h yy (3.22) 这里和下面均记这里和下面均记 ),( kikiki yxfy 上述几种亚当姆斯格式都是上述几种亚当姆斯格式都是显式的显式的，算法比较，算法比较 简单，但用节点简单，但用节点 的斜率值来预报的斜率值来预报 区间区间 上的平均斜率是个上的平均斜率是个外推过程，外推过程，效果不够效果不够 理想。为了进一步改善精度，变外推为理想。为了进一步改善精度，变外推为内插内插，即即 增加节点增加节点xi+1的斜率值来得出的斜率值来得出 上的平均斜率上的平均斜率 。譬如考察形如。譬如考察形如 11 ,xxx ii 1 , ii xx 1 , ii xx ),( ),( )1 ( 1

35、11 11 iii iii iiii yxfy yxfy yyhyy （3.23） 的隐式格式的隐式格式,设设(3.23)右端的右端的 Taylor展开有展开有 )(),( 11 iiii xyyxyy )()()1 ()()( 32 1 hOxyhxyhxyy iiii 可见要使格式可见要使格式(3.23)具有具有二阶精度二阶精度,需令需令 , 这样就可构造二阶隐式亚当姆斯格式这样就可构造二阶隐式亚当姆斯格式 2 1 )( 2 11iiii yy h yy 其实是梯形格式其实是梯形格式 。 类似可导出三阶隐式亚当姆斯格式类似可导出三阶隐式亚当姆斯格式 )85( 12 111 iiiii yy

36、y h yy 和四阶隐式亚当姆斯格式和四阶隐式亚当姆斯格式 )5199( 24 2111 iiiiii yyyy h yy （3.24） 参照改进的欧拉格式的构造方法，以四阶亚当参照改进的欧拉格式的构造方法，以四阶亚当 姆斯为例，将显式（姆斯为例，将显式（3 3.22.22）和隐式（）和隐式（3 3.24.24）相结合）相结合 ，用，用显式公式做预报显式公式做预报，再用隐式公式做校正再用隐式公式做校正，可构，可构 成亚当姆斯预报成亚当姆斯预报- -校正格式校正格式 )9375955( 24 3211 iiiiii yyyy h yy ),( 111 iii yxfy )5199 ( 24 21

37、11 iiiiii yyyy h yy ),( 111 iii yxfy （3 3.25.25） 预报预报 ： 校正校正 ： 亚当姆斯预报校正亚当姆斯预报校正-格式格式 这种预报这种预报- -校正格式是四步法，它在计校正格式是四步法，它在计 算算y yi+1 i+1时不但用到前一步的信息 时不但用到前一步的信息 ，而且，而且 要用到再前面三步的信息要用到再前面三步的信息 ，因，因 此它此它不能自行启动不能自行启动。在实际计算时，可借助。在实际计算时，可借助 于某种单步法，譬如四阶龙格于某种单步法，譬如四阶龙格库塔法提供库塔法提供 开始值开始值 。 ii yy, 321 , iii yyy 32

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

39、结果的精度。 再使用亚当姆斯预报再使用亚当姆斯预报- -校正公式校正公式( (3 3.25),.25), i y 12xey x 稳定性在微分方程的数值解法中是一个非常稳定性在微分方程的数值解法中是一个非常 重要的问题。因为微分方程初值问题的数值方法是重要的问题。因为微分方程初值问题的数值方法是 用差分格式进行计算的，而在差分方程的求解过程用差分格式进行计算的，而在差分方程的求解过程 中，存在着各种计算误差，这些计算误差如中，存在着各种计算误差，这些计算误差如舍入误舍入误 差差等引起的扰动，在传播过程中，可能会等引起的扰动，在传播过程中，可能会大量积累大量积累 ，对计算结果的准确性将产生影响。

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

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

42、递推后，可得将上式反复递推后，可得 0 1 1 )1(yhy i i ), 2 , 1()1 ( 00 iyyhy ii i 或或 h1式中式中 要使要使y yi i有界，其充要条件为有界，其充要条件为 1 11h即即 由于由于 ，故有，故有 0 2 0 h（3 3.26.26） 可见，如欲保证算法的稳定，可见，如欲保证算法的稳定，显式显式EulerEuler格式格式 的步长的步长h h的选取要受到式（的选取要受到式（3 3.26.26）的限制）的限制。 的绝的绝 对值越大，则限制的对值越大，则限制的h h值就越小值就越小。 用隐式用隐式EulerEuler格式，对模型方程格式，对模型方程 的

43、计算公式为，可化为的计算公式为，可化为 )( 11 iii yhyy ii y h y 1 1 1 由于由于 , ,则恒有则恒有 , ,故恒有故恒有 01 1 1 h ii yy 1 因此，隐式因此，隐式EulerEuler格式是绝对稳定的（无条件稳格式是绝对稳定的（无条件稳 定）的（对任何定）的（对任何h0h0）。）。 常微分方程初值问题的求解，是将微分方程转化常微分方程初值问题的求解，是将微分方程转化 为差分方程来求解，并用计算值为差分方程来求解，并用计算值y yi i来近似替代来近似替代y(xy(xi i) ) ， 这 种 近 似 替 代 是 否 合 理 ， 还 须 看 分 割 区 间，

44、 这 种 近 似 替 代 是 否 合 理 ， 还 须 看 分 割 区 间 的 长 度的 长 度 h h 越 来 越 小 时 ， 即越 来 越 小 时 ， 即 时 ，时 ， 是否成立。若成立，则称该方法是收敛的，否则称是否成立。若成立，则称该方法是收敛的，否则称 为不收敛。为不收敛。 这里仍以这里仍以EulerEuler方法为例，来分析其收敛性。方法为例，来分析其收敛性。 EulerEuler格式格式 ii xx, 1 0 1 ii xxh)( ii xyy ),( 1iiii yxhfyy 算法收敛性 设设 表示取表示取 时时, 按按Euler公式的计算结果公式的计算结果, 即即1i y)(

45、ii xyy )(,()( 1iiii xyxhfxyy Euler方法局部截断误差为：方法局部截断误差为： )()( 2 )( 1 2 11 iiii xxy h yxy 设有常数设有常数 ，则，则 )(max 2 1 xyc bxa 2 11 )(chyxy ii （3.27） 总体截断误差总体截断误差 1111111 )()( iiiiiii yyyxyyxy )(,()(),( 11iiiiiiii xyxhfxyyxhfyyy ),()(,()( iiiiii yxfxyxfhyxy （3.28） 又又 由于由于f(x,y)关于关于y满足李普希茨条件满足李普希茨条件,即即 iiiii

46、i yxyLyxfxyxf)(),()( 代入代入(3.28)上式，有上式，有 iiii yxyhLyy )()1( 11 i hL)1( 再利用式（再利用式（3.27），式（），式（3.28） 1111111 )()( iiiiiii yyyxyyxy 2 )1 (chhL i 2 1 )1 (chhL ii 即即 上式反复递推后，可得上式反复递推后，可得 1 0 2 0 )1 ()1 ( i k ki i hLchhL1)1 ()1 ( 0 ii hL L ch hL （3.29） 设设 （T T为常数）为常数） Tihxxi 0 hL ehL 1 TLihLi eehL)1 ( 因为因为

47、 所以所以 把上式代入式（把上式代入式（3 3.29.29），得），得 ) 1( 0 TLTL i e L ch e 若不计初值误差，即若不计初值误差，即 ，则有，则有 0 0 he L c TL i )1( （3 3.30.30） 式式( (3 3.30).30)说明说明, ,当时当时 , , , ,即即 , ,所以所以 EulerEuler方法是收敛的，且其收敛速度为方法是收敛的，且其收敛速度为 ，即具，即具 有一阶收敛速度。有一阶收敛速度。同时还说明同时还说明EulerEuler方法的整体截断方法的整体截断 误差为误差为 ，因此算法的精度为一阶。，因此算法的精度为一阶。 0h0 i )(

48、 ii xyy )(hO )(hO 一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程 我们已介绍了一阶常微分方程初值问题的我们已介绍了一阶常微分方程初值问题的 各种数值解法，对于一阶常微分方程组，可类似得各种数值解法，对于一阶常微分方程组，可类似得 到各种解法，而到各种解法，而高阶常微分方程可转化为一阶常微高阶常微分方程可转化为一阶常微 分方程组来求解。分方程组来求解。 一阶常微分方程组一阶常微分方程组 对于一阶常微分方程组的初值问题对于一阶常微分方程组的初值问题 00 00 )(),( )(),( zxzzyxgz yxyzyxfy （3 3.31.31） 可以把单个方程可以把单个方程

49、中中 的的f f 和和y y看作向量来处理，这样就可把前面介绍的各看作向量来处理，这样就可把前面介绍的各 种差分算法推广到求一阶方程组初值问题中来。种差分算法推广到求一阶方程组初值问题中来。 ),(yxfy ),( ),( ) 2 , 2 ,( ) 2 , 2 ,( ) 2 , 2 ,( ) 2 , 2 ,( ),( ),( 3314 3314 22 2 13 22 2 13 11 2 12 11 2 12 1 1 hLzhKyxgL hLzhKyxfK L h zK h yxgL L h zK h yxfK L h zK h yxgL L h zK h yxfK zyxgL zyxfK ii

50、i iii ii i ii i ii i ii i iii iii 式中式中 （3.34） 设设 为节点上的近似解，为节点上的近似解， 则有改进的则有改进的EulerEuler格式为格式为 iii zyiihxx,);, 3 , 2 , 1( 0 ),( 1iiiii zyxhfyy ),(1 iiii izyxhgzz 预报：预报： ),(),( 2 1 111 i iiiiiii zyxfzyxf h yy ),(),( 2 1 111 i iiiiiii zyxgzyxg h zz 校正：校正： （3 3.32.32） 又，相应的四阶龙格又，相应的四阶龙格库塔格式（经典格式）为库塔格式（

51、经典格式）为 )22( 6 )22( 6 43211 43211 LLLL h zz KKKK h yy ii ii （3 3.33.33） 把节点把节点xi上的上的yi和和zi值代入式值代入式(7.34), 依次算出依次算出 , 然后把它们代入式然后把它们代入式(7.33), 算出节点算出节点xi+1上的上的yi+1 和和zi+1值。值。 对于具有三个或三个以上方程的方程组的初值问题对于具有三个或三个以上方程的方程组的初值问题, 也可用类似方法处理也可用类似方法处理,只是更复杂一些而已。此外只是更复杂一些而已。此外,多步多步 方法也同样可以应用于求解方程组初值问题。方法也同样可以应用于求解方

52、程组初值问题。 44332211 ,LKLKLKLK 例例7.6 用改进的用改进的Euler法求解初值问题法求解初值问题 2)0( 1)0( z z yx z yzxyy 2.00 x 取步长取步长h=0.1，保留六位小数，保留六位小数 。 解解: 改进的改进的Euler法公式为法公式为 ),( 1iiiii zyxhfyy ),(1 iiii izyxhgzz 预报预报 ： ),(),( 2 1 111 i iiiiiii zyxfzyxf h yy ),(),( 2 1 111 i iiiiiii zyxgzyxg h zz 校正：校正： 将将 及及h=0.1代入上式代入上式,得得 z y

53、x zyxgzxyzyxf iiiiii ),(,),( i ii i i iiiii z yx zz zyxyy 1.0 )(1.0 1 1 1 1 1 1 111 05. 0 )()(05. 0 i ii i ii ii i iiiiiii z yx z yx zz zyxzyxyy 由初值由初值 , ,计算得计算得 2)0(, 1)0( 00 zzyy 050000.2 800000.0 1 1 z y 046951. 2) 1 . 0( 801500. 0) 1 . 0( 1 1 zz yy 090992.2 604820.0 2 2 z y 088216. 2)2 . 0( 6046

54、59. 0)2 . 0( 2 2 zz yy 高阶方程组高阶方程组 高阶微分方程高阶微分方程( (或方程组或方程组) )的初值问题的初值问题, ,原则上都原则上都 可以归结为一阶方程组来求解可以归结为一阶方程组来求解。例如。例如, ,有二阶微分方有二阶微分方 程的初值问题程的初值问题 0000 )(,)( ),( yxyyxy yyxfy （3 3.35.35） 在在引入新的变量引入新的变量 后后, ,即化为一阶方程组初值问题即化为一阶方程组初值问题: : yz 0000 )(,)(, ),( yxzyxyzy zyxfz （3 3.36.36） 式（式（3 3.36.36）为一个一阶方程组的

55、初值问题，对此可）为一个一阶方程组的初值问题，对此可 用前面中介绍的方法来求解。例如应用四阶龙格用前面中介绍的方法来求解。例如应用四阶龙格- - 库塔公式得库塔公式得 )22( 6 )22( 6 43211 43211 LLLL h zz KKKK h yy ii ii ),( ) 2 , 2 ,( 2 ) 2 , 2 ,( 2 ),( 3314 34 22 2 13 23 11 2 12 12 1 1 hLzhKyxfL hLzK L h zK h yxfL L h zK L h zK h yxfL L h zK zyxfL zK iii i ii i i ii i i iii i （3 3.37.37） （3 3.38.38） 消去消去 ，上式简化为：，上式简化为： )4 , 3 , 2 , 1( iKi )22( 6 )( 6 43211 321 2 1 LLLL h zz LLL h hzyy ii iii （3 3.39.39） （3 3.40.40） ), 2 ,( ) 2 , 42 ,( ) 2 , 2 ,( ),( 32 2 14 21 2 2 13 1 2 12 1 hLzL h hzyxfL L h zL h z h yxfL L h zz h yxfL zyxfL iiii iii i iii