数值分析第八微分方程

1、二元 Taylor 级数：设 f Cn(a, b × c, d)，(x, y), (x + x, y + y) a, b × c, d。令 F (t) = f (x + tx, y + ty)，则 F (0) = F (x, y)，F (1) = F (x + x, y + y)。f= fx ，f = fy 。当 n = 1 时，记xyF (0) = fx (x, y)x + fy (x, y)y；1/95常微分方程数值分析2fx22f2 fy2= f ，= f 。再记= f，xxxyyyxy1 f+ x, y + y)(x)2R1(1) =2+fxy(x + x, y +

2、 y)(x)(y) +1 fyy(x + x, y + y)(y)2，2其中 0 < < 1。f (x + x, y + y) = F (1)= F (0) + F (0) + R1(1)。误差主项：1F (0)22/95常微分方程数值分析当 n = 2 时，F (0) = f, y)(x)2+2fxy(x, y)(x)(y) + fyy(x, y)(y)2；1 f(x + x, y + y)(x)3R2(1) =6+1fxx y(x + x, y + y)(x)2(y) +21 fxyy(x + x, y + y)(x)(y)22+1fyyy(x + x, y + y)(y)3，

3、6其中 0 < < 1。3/95常微分方程数值分析f (x + x, y + y) = F (1)= F (0) + F (0) + 1F (0) + R2(1)。2误差主项：1F (0)。6误差估计可差高阶无穷小。更一般地，f (x + x, y + y) = F (1)= F (0) + F (0) + F (0) + · · · +12 1F (n1)(0) + Rn(1)。(n1)!4/95常微分方程数值分析第八章 常微分方程的数值解法§1引言科学技术中有很多问题都可用常微分方程的定解问题来描述，主要有初值问题与边值问题两大类。本章只

4、考虑初值问题。dydx= f (x, y), a x b y(x0) = y05/95常微分方程数值分析当 f (x, y) = f (x) 仅与 x 有关时，xf (x)dx。y(x) =x0一阶常微分方程组解在一点 (x0, y0) 附近的性质可通过 f (x, y) 的泰勒展开加以分析。f (x, y) =f (x0, y0) + c(x x0) + (y y0) + · · · ，| 大小反映解对扰动的敏感性。6/95常微分方程数值分析当 为正时，y(x) 随时间 t 增大而增大，是指数增长的；当 为负时，y(x) 随时间 t 增大而减少， 是指数衰减的，

5、是稳定的。7/95常微分方程数值分析虽然求解常微分方程有各种各样的方法，但方法只能用来求解一些特殊类型的方程，实际问题归纳出来的微分方程主要靠数值解法。8/95常微分方程数值分析虽然求解常微分方程有各种各样的方法，但方法只能用来求解一些特殊类型的方程，实际问题归纳出来的微分方程主要靠数值解法。所谓数值解法，就是寻求解 y(x) 在一系列离散节点< · · · 上的近x1 < x2 < · · · < xn< xn+1似值 y1, y2, · · · , yn, yn+1, &

6、#183; · · 。8/95常微分方程数值分析相邻两个节点的间距 hn为步长。今后如不特别 xn 称= xn+1，总假定hi = h(i = 0, 1, · · · ) 为常数。这时xn = x0 + nh, n = 0, 1, 2, · · · 。9/95常微分方程数值分析相邻两个节点的间距 hn为步长。今后如不特别 xn 称= xn+1，总假定hi = h(i = 0, 1, · · · ) 为常数。这时xn = x0 + nh, n = 0, 1, 2, · 

7、3; · 。本章首先对常微分方程离散化，建立求数值解的递推公式。有两类方法，一类是单步法，另一类是 k 步法。9/95常微分方程数值分析§2简单的数值方法一、Euler 公式dydx= f (x, y), a x b y(x0) = y010/95常微分方程数值分析§2简单的数值方法一、Euler 公式dydx= f (x, y), a x by(x ) = y00n = 0, 1, · · ·y= yn + hf (xn, yn),n+1 y(x0) = y010/95常微分方程数值分析若 y0 已知，则可逐次求出y1y2.= y0

8、 + hf (x0, y0)= y1 + hf (x1, y1)11/95常微分方程数值分析若 y0 已知，则可逐次求出y1 = y0 + hf (x0, y0) y2 = y1 + hf (x1, y1).几何意义由 (x0, y0) 出发取曲线 y = y(x) 的切线dy（存在！），则斜率= f (x0, y0) dx(x0,y0)11/95常微分方程数值分析由于 (x0, y0) 及 f (x0, y0) 已知，必有切线方程。由点斜式写出切线方程：dyy = y0 + (x x0)= dx(x0,y0)y0 + (x x0)f (x0, y0)等步长为 h，则 x1 x0算出 y1：y

9、1 = y0 + hf (x0, y0)= h，可由切线12/95常微分方程数值分析逐步计算出 y = y(x) 在 xn+1 点值：yn+1= yn + hf (xn, yn), n = 0, 1, 2, · · ·13/95常微分方程数值分析逐步计算出 y = y(x) 在 xn+1 点值：yn+1= yn + hf (xn, yn), n = 0, 1, 2, · · ·用分段的折线逼近函数, 此为“折线法” 而非“切线法”，除第一段是曲线上的切线，其它都不是。13/95常微分方程数值分析例：求解初值问题2x= y , 0 &l

10、t; x < 1,yy y(0) = 114/95常微分方程数值分析例：求解初值问题2x= y , 0 < x < 1,yyy(0) = 1解：欧拉方法，取步长 h = 0.1，()2x= yn + hy yn+1nnyn精确解 y = 1 + 2x, y(xn) = 1 + 2xn。列表14/95常微分方程数值分析15/95常微分方程数值分析xnyny(xn)0.11.10001.09540.21.19181.1832.0.91.71781.67331.01.78481.7321为了分析计算公式的精度，通常可用泰勒展开将 y(xn+1) 在 xn 处展开，则有y(xn+1)

11、 = y(xn + h)h2 y( ),= y(x ) + y(x )h +nnn2n (xn, xn+1)16/95常微分方程数值分析在 yn = y(xn) 的前提下，f (xn, yn) = f (xn, yn(xn) = y(xn)。于是h2y(x )。欧拉法的误差 y(x) yn+1n+1n2称为此方法的局部截断误差。17/95常微分方程数值分析二、后退 Euler 方法向后差商：y(xn+1)y(xn)y(xn+1) hy= yn + hf (xn+1, yn+1)n+1y(x ) = y00后退 Euler 方法为隐式算法，也称隐式Euler 方法。18/95常微分方程数值分析显

12、式和隐式各有特点，显式易于计算，而隐式算法比较稳定。为避免解非线性方程，与 Euler 法结合：(0)y= yn + hf (xn, yn)n+1y(k+1)= yn + hf (xn+1, y(k) ),n+1n+1 k = 0, 1, 2, · · ·19/95常微分方程数值分析由于 f (x, y) 对 y 满足条件，且yn+1= yn + hf (xn+1, yn+1)(k+1)n+1y y n+1(k)= h f (x, y) f (xn+1, yn+1) n+1n+1(k)n+1 hL y y n+1只要 hL < 1，则上述迭代收敛。20/95

13、常微分方程数值分析考虑 y = yy(0) = 1用向后 Euler 方法产生数值解：y0 = 1, yn+1 = yn + hyn+1因此，在第 n 步，近似解为yn = (1 h)n。对负的 ，要求 |1 h|1切正的步长 h 都成立。< 1。这对一21/95常微分方程数值分析三、梯形公式微分方程 y = f (x, y) 等价于方程xn+1y(x) = y(x ) +f (x, y)dxn+1nxn用梯形公式，且令：= y(xn+1), yn = y(xn)yn+1则得：= yn + h(f (xn, yn) + f (xn+1, yn+1)。yn+12这是二阶隐式算法。22/95

14、常微分方程数值分析同样与 Euler 法结合，形成迭代算法，对n = 0, 1, 2, · · · ，(0)y= yn + hf (xn, yn)n+1()y(k+1)= yn +f (xn, yn) + f (xn+1, y(k) )hn+1n+12 k = 0, 1, 2, · · ·23/95常微分方程数值分析同样与 Euler 法结合，形成迭代算法，对n = 0, 1, 2, · · · ，(0)y= yn + hf (xn, yn)n+1()y(k+1)= yn +f (xn, yn) + f

15、(xn+1, y(k) )hn+1n+12k = 0, 1, 2, · · ·梯形公式的收敛性：(k+1)n+1(k)n+1hL2y yy yn+1 n+1收敛条件 0 < hL/2 < 1。23/95常微分方程数值分析梯形公式比 Euler 法精度要高，但每次迭代均多算一次函数值提高精度的计算代价。24/95常微分方程数值分析四、改进的欧拉公式梯形公式虽然提高了精度，但使算法复杂。而在实际计算中只迭代一次，这样建立的公式。y¯n+1校正校正系统称作改进的欧拉= yn + hf (xn, yn)；= yn + hf (xn, yn) + f

16、(xn+1, y¯n+1)。yn+1225/95常微分方程数值分析这一计算公式也可表为yn+1=yn + hf (xn, yn)+f (xn +h, yn +hf (xn, yn)，2或表为下列平均化形式y = yn + hf (xn, yn);pyc = yn + hf (xn+1, yp); yn+1= (yp + yc)/2.26/95常微分方程数值分析例：用欧拉公式和改进的欧拉公式解初值问题2x= y , (0 < x < 1)yy y(0) = 1.27/95常微分方程数值分析例：用欧拉公式和改进的欧拉公式解初值问题2x, (0 < x < 1)=

17、y yyy(0) = 1.解：取步长 h = 0.1，2xn )。欧拉公式为：yn+1 = yn + h(yn yn27/95常微分方程数值分析改进的欧拉公式为：2xn );y = yn + h(yn pyn2xn+1 );yc= yn + h(yp yp= 1(yp + yc).yn+12列表28/95常微分方程数值分析改进的 Euler 法明显了精度。29/95常微分方程数值分析xnyny(xn)0.11.09591.09540.21.18411.1832.0.91.67821.67331.01.73791.7321五、单步法的局部截断误差与阶dydx= f (x, y), a x b初值

18、问题:y(x ) = y00的单步法可用一般形式表示yn+1= yn + h(xn, yn, yn+1, h)其中多元函数 与 f (x, y) 有关。当 含有 yn+1 时，方法是隐式的；不含则是显式。30/95常微分方程数值分析显式单步法yn+1= yn + h(xn, yn, h)Euler 法：(x, y, h) = f (x, y)。31/95常微分方程数值分析显式单步法yn+1= yn + h(xn, yn, h)Euler 法：(x, y, h) = f (x, y)。显式单步法的局部截断误差定义如下 定义：设 y(x) 是初值问题的准确解，称= y(xn+1) y(xn) h(

19、xn, y(xn), h)Tn+1为局部截断误差。31/95常微分方程数值分析局部截断误差是假设在 xn 之前没有误差，当 yn = y(xn) 时，计算一步的误差。显式 Euler 法的局部截断误差= y(xn+1) y(xn) hf (xn, y(xn)Tn+1= y(xn + h) y(xn) hy(xn)= h2 y(x ) + O(h3)n232/95常微分方程数值分析定义：设 y(x) 是初值问题的准确解，若存在最大整数 p 使显式单步法的局部截断误差满足Tn+1= y(x + h) y(x) h(x, y(x), h) =O(hp+1)则称该方法为 p 阶方法。33/95常微分方

20、程数值分析写：Tn+1 = (xn, y(xn)hp+1 + O(hp+2)，局部截断误差主项：(xn, y(xn)hp+1。以上定义对隐式单步法也适用。34/95常微分方程数值分析例：后退 Euler 法Tn+1= y(xn+1) yn+1= y(xn+1) y(xn) hf (xn+1, y(xn+1)+hf (xn+1, y(xn+1) hf (xn+1, yn+1)注意到|hf (xn+1, y(xn+1) hf (xn+1, yn+1)| hL|y(xn+1) yn+1| = hLTn+1是 Tn+1的高阶无穷小。35/95常微分方程数值分析Tn+1= y(xn+1) y(xn) h

21、f (xn+1, y(xn+1)+O(hp+2)= hy(x ) + h2 y(x ) + O(h3)nn2hy(xn) + hy(xn) + O(h2) + O(hp+2)2h= y (x ) + O(h3) + O(hp+2)n22h故 p = 1，局部截断误差主项 y (x )。n236/95常微分方程数值分析梯形法Tn+1h= y(xn+1) y(xn) y (xn) + y (xn+1)2= hy(x ) + h2 y(x ) +h3 y(x)nnn23!h2 y(x y (x ) + y(x ) + hy(x ) +h)nnnn22+O(h4)3h= y (x ) + O(h4)n

22、123h二阶方法，主项 y (x )。n1237/95常微分方程数值分析改进欧拉公式K1 = f (xn, yn)，K2 = f (xn + h, yn + hK1)，= yn + 1h(K1 + K2)，xn+1= xn + h。yn+12二阶显式算法。38/95常微分方程数值分析K1 = f (xn, yn)，K2 = f (xn + h, yn + hK1)= f (xn, yn) + fx (xn, yn)h+fy (xn, yn)hK1 + O(h2)= f (xn, yn) + fx (xn, yn)h+fy (xn, yn)f (xn, yn)h + O(h2)39/95常微分方

23、程数值分析又y(xn) = f (xn, yn)，y(xn) =f (xn, y(xn)ddx= fx (xn, yn) + fy (xn, yn)f (xn, yn) K2 K1 主项为 y(xn)h。= h(K2 K1) 作为欧拉法的局部以 en+12截断误差。K2 是下一步的 K1。40/95常微分方程数值分析对一个简单问题用 Euler 方法说明数值解法的难点。考虑 y = yy(0) = 1用 Euler 方法产生数值解：y0 = 1, yn+1 = yn + hyn = (1 + h)yn因此，在第 n 步，近似解为yn = (1 + h)n。41/95常微分方程数值分析另一方面，

24、方程实际解为 y = et。若 < 0，则实际解是指数衰减的。当 t 趋于无穷时它趋于 0。实际解是一个稳定状态。而近似解趋于 0 当且仅当 |1 + h| < 1。必须选择步长：0 < h < 2/。42/95常微分方程数值分析例如，若 = 20，必须选择步长h < 0.1。而实际解为y(x) = e20x 2.1 × 109, x 1 时。43/95常微分方程数值分析例如，若 = 20，必须选择步长h < 0.1。而实际解为y(x) = e20x 2.1 × 109, x 1 时。因此，数值解必须小步长进行，而解的性质指出可以采用大步

25、长。Euler 方法需要改进。43/95常微分方程数值分析§3龙格 -库塔方法一、Runge-Kutta 法的基本思想将 Euler 公式与改进 Euler 公式写成下列形式：Euler 公式y= yn + hK1n+1 K1= f (xn, yn)44/95常微分方程数值分析改进 Euler 公式()11y= yn + hK + Kn+11222K1= f (xn, yn)K= f (xn + h, yn + hK1)2以上两组公式都使用函数 f (x, y) 在某些点上的值的线性组合来计算 y(xn+1) 的近似值 yn+1。45/95常微分方程数值分析Euler 公式：每步计算

26、一次 f (x, y) 的值， 为一阶方法。改进 Euler 公式：需计算两次 f (x, y) 的值，二阶方法。46/95常微分方程数值分析一般形式，yn+1= yn + h(xn, yn, h)r(x , y , h) =c Knniii=1K= f (xn, yn)1i1K= f (x + h, y+ hµ K )ininijjj=1(i = 2, 3 · · · , r)这里 ci, i, µij 均为常数。称上述方法为r 级显式龙格 -库塔法，简称 R-K 法。47/95常微分方程数值分析Euler 法：r = 1, (xn, yn,

27、 h) = f (xn, yn)， 阶 p = 1。改进 Euler 法：r = 2，p = 2。要具有更高的阶 p，就要增加点数 r。龙格 -库塔法阶的估计需二元 Taylor 定理：48/95常微分方程数值分析二、二阶显式龙格库塔方法推导 r = 2 的情形y= yn + h(c1K1 + c2K2)n+1K1= f (xn, yn)= f (xn + 2h, yn + µ21hK1)K2局部截断误差= y(xn+1) y(xn) hc1f (xn, yn)Tn+1+c2f (xn + 2h, yn + µ21hK1)49/95常微分方程数值分析为得到阶 p，要将上式各

28、项在 (xn, yn) 展开。记：yn = y(xn), fn = f (xn, yn)。23) = y + hy +hh4y(xy +y+ O(h )n+1nnnn23!其中50/95常微分方程数值分析y = f (xn, yn) = fnny = f (x , y) + f (x , y )fnnnnnnxyyn = fn, yn) + 2fnfxy(xn, yn)2+f f(x , y )nnn xy+f (x , y )f (x , y) + f (x , y )fnnnnnnnyxyf (xn + 2h, yn + µ21hfn)= fn + fx (xn, yn)2h +

29、 fy (xn, yn)µ21hfn+O(h2)代入得51/95常微分方程数值分析= (1 c1 c2)fnTn+1+ c22) fx (xn, yn)h(12()1h + O(h3)2+ c µf (x , y )f221nnny2c + c2 = 11c22 = 1/2c µ= 1/2221有无穷多组解，每一组解得一近似公式，局部截断误差均为 O(h3)。52/95常微分方程数值分析上述方法统称二阶方法。取 c1 = c2 = 1, a2 = b21 = 1，此为改进2Euler 公式。近似公式为y= yn + h(K1 + K2)/2n+1K1= f (xn

30、, yn)= f (xn + h, yn + hK1) K253/95常微分方程数值分析取 c1 = 0, c2 = 1, a2 = b21 = 1，此为常用2的二阶公式，称为中点公式。y= yn + hK2n+1K1= f (xn, yn)= f (xn + h/2, yn + hK1/2) K254/95常微分方程数值分析三、三阶龙格库塔方法类似地，对 r = 3，即三个点，通过更复杂的计算，可导出三阶 RK 公式。常用的三阶 RK 公式为：h K1y= yn + 4K2 + K3)n+16K= f (xn, yn)1K = f (xn + h, yn + hK1)222K= f (xn

31、+ h, yn hK1 + 2hK2)355/95常微分方程数值分析四、四阶龙格库塔方法对 r = 4，即四个点，可导出四阶 RK 公式。常用的四阶 RK 公式为：= yn + h(K1 + 2K2 + 2K3 + K4)yn+16K= f (xn, yn)1K2 = f (xn + h, yn + hK1)22= f (xn + h, yn + hK2)K322K= f (xn + h, yn + hK3)456/95常微分方程数值分析例：设取步长 h = 0.2，从 x = 0 直到x = 1 用四阶龙格库塔方法求解初值问题2x= y , (0 < x < 1);yy y(0)

32、 = 1.57/95常微分方程数值分析例：设取步长 h = 0.2，从 x = 0 直到x = 1 用四阶龙格库塔方法求解初值问题2x= y , (0 < x < 1);yyy(0) = 1.解：由经典的四阶龙格库塔公式得57/95常微分方程数值分析 yh= yn +(K1 + 2K2 + 2K3 + K4);n+162xn= yn = yn +K;1ynh2x + hh2nK2K1 +K1;2ynh2xn + hhK= yn +K2 2+K2;23yn2(xn + h)= yn + hK3 yK.4+ hKn358/95常微分方程数值分析两点说明：1） 当 r = 1, 2, 3

33、, 4 时，RK 公式的最高阶数恰好是 r，当 r > 4 时，RK 公式的最高阶数不是 r，如 r = 5 时仍为 4，r = 6 时 RK 公式的最高阶数为 5。2） RK 方法的导出基于 Taylor 展开，故要求所求问题的解具有较高的光滑度。59/95常微分方程数值分析当解充分光滑时，四阶 RK 方法确实优于改进 Euler 法。对一般实际问题，四阶RK 方法一般可达到精度要求。如果解的光滑性差，则用四阶 RK 方法解的效果不如改进 Euler 法。Runge-Kutta 方法也分显式和隐式。60/95常微分方程数值分析r 级显式如下：ry= yn + hciKin+1i=1K1

34、= f (xn, yn)i1K= f (x + a h, y+ hb K )ininijji=1(i = 2, 3 · · · , r)61/95常微分方程数值分析r 级隐式：ry= y + hc Kn+1niii=1rKi= f (xn + aih, yn + hbijKj)j=1(i = 2, 3 · · · , r)系数矩阵 bij 为 r 阶方阵，显式时为严格上三角。r 级显式 RK 方法阶数 r。通常使用不超过 4 阶方法。62/95常微分方程数值分析隐式 RK 方法：梯形方法：隐式 RK 方法（二阶）()+KK12y= y

35、n + h+n+122K1= f (xn, yn)= f (x, y)hK1+hK2 K2n+1n2263/95常微分方程数值分析一级二阶隐式中点方法 y= yn + hK1n+1()hhK1K= fxn + 2 , yn +12或写成= y + hf (xn)hyn + yn+1y+,2n+1n2四阶隐式 RK 方法。此外，还有隐式 RK 方法有较好的稳定性，是目前求解刚性方程组的重要方法。64/95常微分方程数值分析龙格 -库塔法：()i1，Ki = ftn + aih, yn + hbijKjj=1i = 1, 2, · · · , kk= yn + hci

36、Ki。yn+1误差en+1i=1k= hdiKi。i=165/95常微分方程数值分析的龙格 -库塔公式：ode45，ode23。ode23 同时包括一个二阶和三阶龙格 -库 塔公式。其差 en+1 作为选取下一步步长的依据。的多步法：ode113，· · ·66/95常微分方程数值分析f=(t,x)2*t*x; t,x=ode23(f,0,1,1);plot(t,x)显示解：x = et2 。注：当 f略。与 t 无关时，变量 t 也不能省67/95常微分方程数值分析68/95常微分方程数值分析单步法用于实验方程 y = y。= E(h)yn，yn+1E(h) 依

37、赖于所选的方法。欧拉法：E(h) = 1 + h； 经典四阶龙格 -库塔法：(h)2(h)3(h)4E(h) = 1 + h +2!3!4!69/95常微分方程数值分析当 |E(h)| < 1 时，方法绝对稳定。复变量 h，绝对稳定区域。实变量 h，绝对稳定区间。欧拉法绝对稳定区间 (2, 0)。经典四阶龙格 -库塔法绝对稳定区间(2.785, 0)。70/95常微分方程数值分析有各种常微分方程的解算器，例：ode23，ode45。用法：第一个参数：函数句柄或第二个参数：t0, tf ，tf 可以是 Inf；第三个参数：初始值 y0。函数；区间；相对误差，绝对误差有默认值。71/95常微

38、分方程数值分析误差是 ode23，ode45 的一个属性，可用odeset 修改。opts = odeset(reltol,1.e-5, abstol,1.e-8, . outputfcn,odeplot)修改了相对误差，绝对误差，添加了每步后执行的函数。此时，语句 ode23(F, t0, tf , y0, opts)，产生图形。72/95常微分方程数值分析Events 是 ode23，ode45 的另一个重要属性。实际问题，往往终止时间是需要求出的。例：y¨ = 1 + y2，y(0) = 1，y(0) = 0。 有空气阻力的自由落体模型。当 y(t) = 0时，过程终止。73/

39、95常微分方程数值分析例：地球绕太阳轨道。旋转一周即结束。Events 属性需调用一个函数，ode23 或ode45 每执行一步，有三个输出量。74/95常微分方程数值分析自由落体模型 Events：function gs,isterm,dir gs = y(1);isterm = 1;=g(t,y)dir=;75/95常微分方程数值分析第一个输出量是一个函数，终止时是零；第二个输出量为 1 时，表示方程在函数gs 为零时，应该结束，为 0 时，不应该结束；第三个输出量，为 0 或 时，表示可以从任意方向接近零；为 1 时，从上方（正的方向）趋于零；为 1 时，从下方（负的方向）趋于零。76/

40、95常微分方程数值分析§4线性多步法在逐步推进的求解过程中，计算 yn+1前事实上已经求出一系列近似值之y0, y1, . . . , yn，如果充分利用前面多步的信息来yn+1，则可期望获得较高的精度，这就是多步法的基本思想。77/95常微分方程数值分析一、线性多步法的一般公式如果计算 yn+k 时，要用到前 k 步的信息yn, yn+1, . . . , yn+k1，则称此方法为线性 k 步法。一般公式为：k1kyn+k=iyn+i+ hifn+ii=0i=0其中 i, i 均为常数，= f (xn+i, yn+i)，(i = 0, 1, . . . , k)。fn+i若 k =

41、 0，显式；k = 0，隐式。78/95常微分方程数值分析定义：设 y(x) 是初值问题的准确解，线性多步法在 xn+k 的局部截断误差为Tn+k=k1kiy(xn+i) hiy (xn+i)。y(xn+k) i=0i=0若 Tn+k的。= O(hp+1)，则称该方法是 p 阶79/95常微分方程数值分析Taylor 展开：y(xn+i) = y(xn+ih) = y(xn)+ihy(xn)+· · ·y(xn+i) = y(xn + ih) =y(xn) + ihy(xn) + · · ·代入到 Tn+k，得Tn+k=c0y(xn)

42、 + c1hy(xn) +· · · + cphpy(p)(xn) +· · ·80/95常微分方程数值分析其中c0 = 1 (0 + · · · + k1)，c1 = k 1 · 1 + · · · + (k 1)k1(0 + · · · + k)，q!cq = kq 1q1 + · · · + (k 1)qk1q1q1 · 1 + · · · + kq1k，q =

43、 2, 3, . . .81/95常微分方程数值分析若在公式中选取系数 i 及 i，使它满足c0 = c1 = · · · = cp = 0，cp+1 = 0。局部截断误差= cp+1hp+1y(p+1)(xn) + O(hp+2)Tn+k称右端第一项为局部截断误差主项，cp+1为误差系数。82/95常微分方程数值分析例：k = 1 时，若 1 = 0，则 0 = 1，0 = 1。c0 = c1 = 0，得欧拉法：yn+1 = yn + hfn。此时，c2 = 1/2，局部截断误差为12 3=h y (xn) + O(h )。Tn+12一阶精度。83/95常微分方程数值分析例：k = 1 时，若 1 = 0，隐式