第七章(微分方程-3)

1、第七章常微分方程的数值解法§ 1引言1、 一阶常微分方程初值问题(微分方程加初值条件)y f(x,y),a x by(x。)yo特征：f(x, y)为的二元函数，只有一个变元x，一阶导数y' dy，要求满足微分方dx程且过点(xoy。)的曲线，一个未知函数y y(x)(微分方程的解函数)。2、解的存在性假设f (x, y)连续且对y满足Lipschitz条件，即存在常数 L 0，使对 力皿 R，有 f (x, y1) f (x, y2) L y1 y2，那么初值问题 的解存在唯一，本章均假设 f (x, y)对y满足 Lipschitz 条件。3、数值解法问题：求y y(x)

2、精确解(解析表达式)极其困难，实际应用中只要求数值解。数值解：在a,b取X。 X1 X2Xn Xn 1 系列等距离散节点，步长h Xn 1 Xn , n 0,1,，求 y(Xn)的近似值 yn。XoXixohx2 xo 2hXn Xo nhXn 1Xo (n 1)hy(xo)ygyX)y(xn)y(xn 1)=yoy2ynyn 1方法：建立yn yn 1的递推公式，y(xo) yo，从而yoy1yn y 1按节点排列顺序(步进式)。§ 2 Euler 方法1、显示、隐式Euler法与梯形法Xn i对(1)的方程两端由Xn到Xn i定积分得：y(xn J y(xn)f (x, y(x)

3、dxxn从而有：y(Xn 1)y(Xn)hf ( ,y( ),XnXn 1 ( *)Xn 1(积分中值定理 x f (x,y(x)dx (Xn 1 Xn)f( ,y( ) h f( , y( ),XnXn 1 )x()假设取 Xn，得 f(x,y)dx h f (xn,y(xn)(左矩形公式)xn由(* )得:y(Xn 1) y(Xn)h f(Xn, y(Xn)Yny(Xn),yn1y(Xn1),那么有：yn1 ynhf (Xn,yn), n0,1,称式为解初值问题的Euler法(显格式)YoY1ynyn 1几何意义：初值问题的解曲线y y(x)过点Po(xo,Yo)，从Xo出发，以f(Xo,

4、y。)为斜率作一 段直线，与直线 X X1交点于R(X1,Y1)，显然有Y1Yo hf (Xo, Yo)，再从R出发，以f(X1,yJ为斜率作直线推进到X X2上一点P2，其余类推，这样得到解曲线的一条近似曲线，它就是折线P0RP2(某点导数值的几何意义为该点切线的斜率)。(2)假设取Xn 1，得 n 1 f (X, y)dx h f(Xn 1, y(Xn 1)(右矩形公式)Cn类似可得：Yn 1 Yn hf(Xn1,Yn1), n 0,1,称式为隐式Euler法(向后Euler法，隐格式)Xn 1h(3)取 f (X, y)dx - f (Xn, Y(Xn)f (Xn 1, y(Xn 1 )

5、(梯形公式)Cn2那么得:Yn 1hYn f(Xn, Yn)f ( Xn 1, Yn 1 ), n 0,1,2称式为 梯形法 .(隐格式)0.1 ，例1用Euler法、隐式Euler法、梯形法及改进Euler法求解y' y x 1,y(0)1，取h计算到 x 0.5，并与精确解比较 .解：由于 f (x,y) y x 1,h 0.1 x00,y01Euler 法：yn 1yn h( ynXn 1)(1 h)yn hXn h 0.9yn O.IXn 0.1隐式 Euler 法: yn 1Ynh(Yn 1Xn 11)1解出Yn 1(Yn1 hXn1h)(Yn1.10.1xn 0.11)1当

6、n 0时，Y1'(Y00.1X00.11)1.00901n 0时，y10.9y0 O.1xO 0.1 1.00000.1.1h梯形法：yn 1yn -f(Xn,yn) f (Xn 1, yn 1)2hYn 1 Yn ( Yn Xn 1)(% 1 Xn 11)21 1解出yn1(2h)ynh(XnXnh)2h(1.9yn0.2Xn0.21)2 h2.11当 n 0时，y1(1.9 0.21) 1.0047622.1改进Euler法：yn 1yn f (Xn,yn) f X 1nhf X，n )h即 yn 1yn( yn Xn 1)n h(n X. 1) X. 1 121h(2h(1 h)

7、2Xn2(Xn2h(2 h)20.905 yn 0.095xn 0.1当n 0时,y10.905 y00.095x00.11.005000精确解：yX x e X,梯形法效果最好,!)X0Euler 法 yn隐式 Euler法yn梯形法Yn改进Euler法yn精确解yXn0.1 11.000 0001.009 0911.004 7621.005 0001.004 837:0.21.010 0001.026 4461.018 5941.019 0251.019 7310.31.029 0001.051 3151.040 6331.041 2181.040 8180.41.056 1001.083

8、 0141.070 0971.070 8021.070 3200.51.090 4901.120 9221.106 2781.107 0761.106 531改进Euler较好：具有相同的误差数量级,其它不好2、隐式法yn 1的计算Euler法及梯形法方法1:显示化f对y线性时，以yn 1为未知量的一元线性方程，见上例 方法2:迭代法f对y非线性时，可看作一个关于 yn 1方程，利用迭代法求解Yn 1)Euler 法：yn 1 yn hf (Xn 1,yn 1),n0,1,(*1)以yn 1为未知量的一兀非线性方程(k 1)y f (Xn(0) yn 1 yn1, y(k),k0,1,h(*3

9、 )(*4 )梯形法：yn1 yn 畀X5 fXni，yni，n 0,1,ynvyn £f(Xn,yn)彳风小补叫)k 01 yi yn hf(xn,yn)计算步骤：(1) 初值 y10)yo hf(Xo, yo)；迭代yi(0)yiyi(k)y?1yi(当 |y1k1)yj)(2) y20)y1 hf(X1，y1)；迭代 y2°) yjy)k)y；k ° 纸只要步长h足够小，就可保证迭代收敛1a当h *2收敛L证明：b迭代序列ynk1收敛收敛性:Jk 1)-hyn 1yn 12f (Xn 1,ynk)1)f (XnO2(b)当 h -L,有Yn)-L21, yn

10、(*4)收敛.(k)1yn 1yn(k)1 1yn 10 ;由*3 、*4式知:满足李氏条件、，k、，(hL)k、，° 、，y n 1yn 1y n 1yn 12yn0)1yn 10，从而收敛得证。(即1hL 1时有3、 改进Euler法梯形法隐式迭代式仅限迭代一次，为防止迭代，可先用 Euler法计算出yn1的近似yn 1yn 1yn hf Xn, ynh显式预报，隐式校正yn 1ynfXn,yn f Xn 1,yn 12h或yn 1 y fXn,yn f Xn 1, yn hf Xn , yn 62称以上两式为改进Euler.法校正就是迭代一次 右端已不含yn 1 ,显式方法。例

11、2:用梯形法的迭代格式求y' y ，y0 1的数值解，h 0.1，计算到x3 0.3解：梯形公式yn 1hyn -(yn2Xn)y(yn 12xn 1n 1)0.05yn1 0.1Xn 1Xn1.05yn 0.12ynyn 1yn 1yn(k 1) yn 10.05ynk1 0.1X n 1(k)1.05yXnn0.1Y11.095656迭代式yn 1y n，解为：Y21.183594(0) yn 1ynhf(Xn,yn)1.1ynXn0.2ynY31.4150593、 Euler法的局部截断误差定义1在yn y(Xn)的前提下，称y(Xn 1)y“i为在Xn 1的局部截断误差定义2假

12、设一种数值方法的局部截断误差：y(Xn i) yn 1 o(hp1)，那么称这种方法是精 度为p阶的。(含hp 1的项称为局部截断误差主项.)注：按某种方法由yn算出yn 1这一步的误差。(1)Euler 法y(Xn1) y(Xn h)ynynf(x, y(x)1ynh f(Xn, yn)y(xn), y (x)由(*)(* )式有:y(Xn 1)' h2y(Xn) h y (Xn) y (Xn)(*)21y (Xn ) f (Xn, y(Xn)h f(Xn, y(Xn)y(Xn) h y'(Xn) ( * )h2“y (Xn)2y(Xn)yn 1O(h2)根据定义2, Eul

13、er法中的p 1故此方法为一阶方法隐式Euler法(分析时隐式公式右边的yn 1 y(xn 1)y(Xn), y'(x)f(X, y(x)y'(Xn)f(Xn,y(Xn)yn h f(Xn 1, yn 1)yX) h f(X.y(Xn) h.y'(Xn h)'2y(Xn) h y (Xn) hynyny(Xn) hy"(Xn)由(*)及(*)式:y(Xn 1)y'(Xn) hy"(* ) h2 2yyn故局部截断误差主项为与 y'd),21,(3)梯形法y(Xn 1)y(Xnh)y(Xn) hy (Xn)yn 1 ynh一f(

14、Xn, yn)2f (Xn1, y(xn类似有y(x1) yn 1h3，“12 y(Xn)(4)改进Euler法y'(x)那么 y"(Xn)fx(Xn，yn)f(Xn 1,ynf(X, y)hf (Xn，yn)(x)fx (x, y)fy(Xn, Yn) Yf(Xn h,yn1,y(xn(Xn)(Xn)1) y(Xn) h y (Xn 1)也是一阶方法.O(h2)u2. 3h- h'"y (Xn) y (Xn)23!1)O(h3)，方法是二阶的.fy(X, y) y (X)(Xn)hy'(Xn)f(Xn°n) fx(Xn,yn) h fy(

15、Xn，Yn) 1'区皿) y'(Xn) hy''(Xn)y (Xn) h y(Xn)hyn 1 yn f(Xn,yn)2h '= y(Xn) -y (Xn)2f (Xn 1, yn hf (Xn, yn)y(Xn)h y (Xn)h2''Ty(Xn)y(Xn 1) y(Xnh)y(Xn) h y (Xn)h2Ty(Xn)£3!y (Xn)y(Xn 1) yn 1O(h3),方法是二阶的.§ 3 Runge-Kutta 方法显式Runge-Kutta法的一般形式y iyn h(Gki C2k2ki f (Xn,yn)般形

16、式:k2f(Xna2h, ynb2iki)krf(Xn arh,yn hbnk1hbr ,r 1kr 1)2级显式R-K方法(r=2)Yn 1 yn h(Gk1 qk2)设想构造R K公式：k1f (Xn, yn)(*)k2f (Xn a2h, yn b2*2)试确定G,C2,a2,b21，使上式为二阶格式，即y(Xn 1) yn 1 O(h3):假设y(Xn)yn，对(7.3.6)式在(Xn, yn)处按Taylor公式展开，由于y'(Xn)y (Xn)f(Xn，y(Xn)df(x,y(x)f (XnWn)；11f (7"(7 WfW/X, W f(7WdXx Xn収八n，

17、八八n1 y 1，y & 丨 '八n，y '八nfx(Xnn)fy(Xn,yn)f (Xn,yn)'h2"h '"y(Xn 1)y(xn h)y(Xn) hy (Xn)y (xn) =y ()23!h；211h3，ynhf (Xn,yn)-2-fx(Xn,yn) fy(Xn,yn)f (Xn, yn) y ()k2f(Xnh,ynkJf (Xn,yn)a2hfx(Xn,yn)b21k1fy(Xn,yn)O(a；h2b：k2)O(h3)将上述结果代入(*)得：y(Xn 1)yn1(1cC2)hf(Xn,yn)(-。2玄2)fxX,n)(

18、-妙21)f y(Xn, y.)f X,n)2 2要使公式为二阶，必须：c1c211ga21(三个方程四个求知数，无穷多解)21C2b212yn 1 yn(K k2)1 2(1)取 C1 -,a2 b21 1，有k1 f(Xn,yn)(改进欧拉法)2k2 f (Xn 1,yn hK)1(2)取 c10, a2b21，有2k2kif(Xn,yn)(中点法)f (Xn - , yn22取其他数时，也可得到其他公式，但系数较复杂，一般不再给出 四阶R-K方法及步长的自动选择类似可得到，经典的四阶 R-K方法是：、yn 1yn(k162k2 2k3 k4)k1f (Xn, yn)hhhhk2f (Xn

19、2,yn2k1)k3f (xn, ynk2)2 2k4f (Xnh,ynhk3)它的局部截断误差0(h5)，故p=4,这是最常用的四阶R-K方法，数学库中都有用此方法 求解初值问题的软件.这种方法的优点是精度较高，缺点是每步要算4个右端函数值，计算量较大例7.3 用经典四阶R-K方法解例7.1的初值问题y'y x 1, y(0) 1，取h 0.1 ,并与改进Euler法、梯形法在处比较其误差大小，于是解 用四阶 R-K 方法公式 (7.3.12)，此处 当 n=0 时于是，按公式 (7.3.12)可算出此方法误差：改进 Euler法误差:梯形法误差： 可见四阶R-K方法的精度比二阶方法

20、高得多.y(Xn) yn的估计式 h?通常是通过不同步用四阶R-K方法求解初值问题(1)精度较高，但要从理论上给出误差 那么比较困难.那么应如何判断计算结果的精度以及如何选择适宜的步长长在计算机上的 计算结果近似估计.设 在处的值当 时， 的近似为y(xn Jynh)iCnh5假设以h为步长'计算两步到Xn 1，那么有(h)y(Xn 1) yn21于是得:即或y(Xn i)ynh)iy(Xn l)y(Xn l)(h)yn 1h(2)yn 1珈216ynh)1)y(Xn i)yn1(h) yn 1(7313)1 一 它给出了误差的近似估计.如果-yn2115hh果y：2；满足精度要求，假

21、设 丄y；2； ynh115曲站给定精度，那么认为以号为步长的计算结，那么还可放大步长因此提供了自动选择步长的方法.，于是由四阶R-K方法有§ 4线性多步法线性多步法的一般公式前面给出了求解初值问题(1)的单步法，其特点是计算yn1时只用到yn的值(要计算多点 处的函数值，计算量较大)，此时y。，, , yn 1, yn的值均已算出.如果在计算y* 1时除用yn的值 ,yn r 1的值，这就是多步法.假设记Xk Xo kh , h为步长， 01,n)，那么线性多步法可表示为(a r yn r 1 h( 1:r 10，称(7.5.1)为线性r步法.计算时用到前面已算出的r个0时，为显式

22、方法，当10那么称为隐式多步法.隐式方法与梯形外,yk其中值.y还用到yn 1 , yn 2, y(Xk), fkf (Xk,yk)(ky n 1 a0 yn a1 yn 1i, i为常数，假设2 r n, yn 1, yn r 1 当 1Xo般形式)：r f n r 1 )方法一样，计算时要用迭代法求yn1多步法(7.5.1)的局部截断误差定义也与单步法类似. 线性多步法在xn 1处的局部截断误差定义为：假定 yi y(x),i n,n 1, ,n r 1时 假设y(Xn 1) yn 1 O(hp 1)，称其为是p阶的。二、Adams显式与隐式方法y n 1 y n h( 1 f n 10

23、f n 1 f n 1r f n r 1 )(即 11， i 0, i 2, ,n r 1)的r步法称为Adams方法，当 1 0时为Adams显式方法，当1 0时，称为Adams隐式方法.对初值问题(1)的方程两端从xn到xn 1积分得Xn 1y(Xn 1) y(Xn) x f(x,y(x)dxxn选取不同的节点，对上式中被积函数进行不同的多项式插值，从而得到不同的的线性多步法。1、Adams四步显式公式(外插公式)选(xn, f(xn,y(xn)，( Xn1, f (Xn1, y(Xn1)，( Xn2, f (Xn 2, y(Xn 2 )，(Xn 3, f (Xn 3, y(Xn 3)四节

24、点作三次多项式插值：f (X, y(x)卩3&) R3 (x)从而 y(Xn 1) y(Xn)冷 1 R(x) R3(X)dXXn取Xn 3,Xn 2,Xn 1与Xn为插值节点做f(X, y(X)的插值多项式：P3(X)(X Xn 1)(x Xn 2)(X Xn 3)(Xn XnJ(XnXn 2)( XnXn 3)f (Xn, y(Xn)(X Xn)(X Xn2)(X Xn 3)(Xn 1Xn)(Xn 1Xn 2)(Xn 1Xn 3)f (Xn 1, y(Xn 1)(X Xn)(XXn 1)(XXn 3 )(Xn 2Xn)(Xn 2 XnJ(Xn2 Xn 3)f (Xn 2, y(Xn

25、 2)(X Xn)(X XnJ(X Xn 2)(Xn 3Xn)(Xn 3Xn 1)(Xn 3Xn 2)f (Xn 3, y(Xn 3)用 P3(x)代替 f (x, y(x)y(Xni) y(Xn)Xn 1XnP3 (x)dx,Xn th，并有 Xn Xn 1Xn 1 Xn 2Xn 2 Xn 3h11y(Xnl) y(Xn)6【f(Xn，y(Xn) 0(t 1)(t2)(t 3)dt f 区 1, 丫区 J) °t(t 2)(t1 1f(Xn 2,y(Xn 2) 0t(t 1)(t 3)dt f (Xn 3(Xn 3) °t(t 1)(t2)dt3)dty(Xn 1)y(X

26、n)刃55f(Xn,y(Xn) 59f(Xn1，y(Xn1)37 f (Xn 2, y(Xn 2)9f (Xn 3, y(Xn 3)令yi代替y(Xi),i n 1,n,n 1,n 2,n 3，从而得到四步Admas (外插)公式那么有:hyn 1 yn 一 55f(Xn,yn) 59彳区1皿1)37彳区2皿2)9心33)24yn55fn59 fn137 fn 29 fn3(注：记fif(Xj,yJ,i n, n 1, n 2,n3)24r步的Adams方法计算时必须先用其他方法求出前面r个初值y. 1, yn 2, , yn r 1才能按 给定公式算出后面各点的值，它每步只需计算一个新的f值

27、，计算量少，但改变步长时前面的yn 1, yn 2, yn r 1也要跟着重算，不如单步法简便局部截断误差：当 y(Xi) y,iYn 1 Ynn,n 1,n2,n 3时，由上知:Xn 1F3 (x)dxXny(Xn 1)从而：y(xnR3(X)Xn 1x P3(x) R3(x)dxxnXn 11)yn 1R3(X)dXn(4)1 d()f(x, y(x)y(Xn)4! dx41 d(4)(y'(x)(X Xn)(X XnJ(X Xn2)(XXn 3)Xn 1X (X4!dx44/)()(x 皿Xn)(XXm)(XXn2)(XXnXn 1)(X Xn 2)(X Xn 3)3)R3 (x

28、)dxXnXXn,Xn令 X Xn th，Xn 11(5)xy ( )(X Xn)(X Xn J(X Xn 2)(XXn 4!,(X Xn)(X Xn 1)(X Xn2)(X Xn 3)0 保号Xn 2 Xn 2 Xn 3，从而1并有Xny(Xn 1) yn 1XnXn因此是四阶方法R3(x)dxXn 1 Xnh5(5)(4!Xn 3)dX119MW 2)(t 3)dt 莎h5y(5)()5O(h )2、Adams三步隐式公式(内插公式)选(Xn1, f (Xn 1, y(Xn1),(Xn, f (Xn, y(Xn),(Xn1, f (Xn 1, y(Xn 1),hyn 1 yn 55f(Xn

29、,yn)24局部截断误差：y(Xn 1) yn 159f (Xn 1, yn 1)37f (Xn 2 , Yn 2)9f (Xn 3, yn 3)Xn 1xR3(x)dxx n251颓h5y () O(h5)，三步四阶隐式方法(xn 2, f (xn 2, y(Xn 2),四节点作三次多项式插值：f(X,y(X) P3(X)只3( X) 类似可得：y(Xn 1 )h必)2455f(Xn,y(Xn) 59f(Xn 小区 J) 37心小风 J) 9f(Xn 3,y(xn 3)用yn代替y(Xn 1)那么有:、Adams预测-校正方法Adams显式方法计算简单，但精度比隐式方法差，而隐式方法由于每步

30、要做迭代，计算不方便为了防止迭代，通常可将同阶的显式Adams方法与隐式Adams方法结合，组成预测-校正方法.以四阶方法为例，可用显式方法()计算初始近似yn°)i，这个步骤称为预测 (Predictor)，以P表示，接着计算f值(Evaluatio n),曾f(Xni,yno)i)，这个步骤用E表示，然后用隐式公式(7.5.10)计算yn 1，称为校正(Corrector)，以C表示，最后再计算，为下一步计算做准备.整个算法如下:h预测:yn 1 yn (55fn 59 fn 1 37 fn 2 9 fn 2)24_求值:fn 1 f ( Xn 1, yn 1)校正:yn1yn

31、(9fn119 fn5fn_1fn2)24 求值:yn 1 f(Xn, yn 1)上式称为四阶Adams预测-校正方法(PECE).例7.6用四阶显式Adams方法及四阶隐式Adams方法解初值问题，步长h=0.1计算得到.用到的初始值由精确解解 此题直接由公式 (7.5.8)及 (7.5.10)计算得到 .对于显式方法，将直接代入式 (7.5.8)得到其中对于隐式方法，由式 (7.5.10)可得到直接求出，而不用迭代，得到计算结果如表所示§ 5 单步法的收敛性与绝对稳定性是单单步法的收敛性定义4.1设y(x)是初值问题 的精确解,步 法 (7.3.2) 在处 产 生 的 近 似 解

32、 , 假设收敛于那么 称 方 法 (7.3.2)产 生 的 数 值 解实际上，定义中是一固定点，当 hO时nfg, n不是固定的.因显然方法收敛那么在固定点处的整体误差当p>1时面定理给出方法 收敛的条件 .定理4.1 设初值问题的单步法732是p阶方法p > 1,且函数对y满足Lipschitz条件，即存在常数L > 0,使对，均有定理证明略 .可见 3那么方法 (7.3.2)收敛，且二、 绝对稳定性，由于原始数据及计用单步法 (7.3.2)求数值解算过程舍入误差影响,实际得到的不是而是，其中误差，再计算下一步得到以 Euler 法为例，假设令,那么(7.4.1)如果那么从

33、冷口误，用 Euler 法求计算到差不增长，它是稳定的 .但如果条件不满足就不稳定例 7.4 y' =00y,y(0)=1,精确解为 解得假设 取 h=0.025 ， 那么，而 ，显 然计算是不稳定的 .如果用后退Euler法(725)解此例，仍取h=0.025，那么显然当，计算是稳定的有关，在此由此看到稳定性与方法有关，也与例中.在研究方法的稳定性时，通常不必对一般的f(x,y) 进行讨论，而只针对模型方程这里(7.4.2)可能为复数.规是因为 时 微分方程 (7.4.2)本身是不稳定的，而讨论数值方法 (7.3.2)的稳定性，必须在微分方程本身稳定 的 前 提 下 进 行 . 另 一 方 面 ， 对 初 值 问 题 (1) ， 假设 将 f(x,y) 在处线性展开，可得于是方程 (1)可近似表示为它说明用模型方程(742)是合理的，至于模型方程(742)中所以用复数 入是因为初值问题 如 果 是 方 程 组 ， 即， 那么是(mXm)阶矩阵，其特征值可能是复数.当然对单个方程，