浙江大学：《数值分析》课程PPT教学课件（双语版）第九章 常微分方程数值解 Numerical Methods for Ordinary Differential Equations §4 线性多步法 Multistep Method

§4线性多步法My当时,为隐式公 用若干节点处的式;B则为显式公式 y(x;)。其通式可写为: yn=any2+ay1+…+yx+M(BJm1+月+月+…+Bf) >基于数值积分的构造法 将y=f(x,y)在[x;,xl上积分,得到 y(xi+1)-V(xi )=.f(x, y(x)dx 只要近似地算出右边的积分lk≈」f(x,(x),则可通 过y=y1+k近似y(x)。而选用不同近似式k,可得到不 同的计算公式

§4 线性多步法 /* Multistep Method */ 用若干节点处的 y 及 y’ 值的线性组合来近似 y(xi+1)。 ... ( ... ) i 1 0 i 1 i 1 k i k 1 i 1 0 i 1 i 1 k i k y y y y h f f f f + - - - + - - =a +a + +a + b + b + b + + b 其通式可写为： ( , ) j j j f = f x y 当 b-10 时，为隐式公 式; b-1=0 则为显式公式。 ➢ 基于数值积分的构造法 将 y = f (x, y) 在 [xi , xi+1 ] 上积分，得到  + + - = 1 ( ) ( ) ( , ( )) 1 i i x x y xi y xi f x y x dx 只要近似地算出右边的积分 ，则可通 过 近似y(xi+1) 。而选用不同近似式 Ik，可得到不 同的计算公式。  +  1 ( , ( )) i i x x I k f x y x dx i i k y = y + I +1

§4 Multistep method a亚当姆斯显式公式/ Adams explicit formulae 利用k+1个节点上的被积函数值f1,f1,…,f构造k阶牛顿 后插多项式N(x2+th),t∈[0,1,有 f(x,y(x)dx=N(x, +th)hdt+R,(x;+th)hdt P=男+小N4+2/显式计算公式eon 局部截断误差为:R=y(x1)-yn=hfR(x)+t塘余顶 例:k=1时有M(x2+th)=f+tV=f+(1-f1) h =y++(f-1)=y+(3f-fa) R d∫(x,y(5)1 2! th(t+1)ht=;h3y"(51)

§4 Multistep Method  亚当姆斯显式公式 /* Adams explicit formulae */ 利用k+1 个节点上的被积函数值 构造 k 阶牛顿 后插多项式 , 有 i i i k f f f - - , , ... , 1 N (x + th) , t [0, 1] k i    = + + + + 1 0 1 0 ( , ( )) ( ) ( ) 1 f x y x dx N x th hdt R x th hdt k i x x k i i i Newton 插值余项 y y h N x th dt i i k i ( ) 1 0 1 = + + +  /* 显式计算公式 */ 局部截断误差为：  = + - + = + 1 0 1 1 R y(x ) y h R (x th)dt i i i k i 例：k=1 时有 ( ) ( ) 1 i + = i +  i = i + i - i-1 N x th f t f f t f f + = +  + - - = + - - 1 0 1 1 1 (3 ) 2 [ ( )] i i i i i i i i f f h y y h f t f f dt y th t h dt dx d f y R h x x i ( 1) 2! 1 ( , ( )) 1 0 2 2 = +    ( ) 12 5 3 i = h y 

4 Multistep method 注:一般有R=B+y4+2(5),其中B与y1计算公式中 f,…,k各项的系数均可查表得到。 k f f f f B 1 Misprint 12 23 16 on p. 204 12 55 59 24 24 24 24 720 纔常用的是k=3的阶亚当姆斯显式公式 h y社1=y+(55f-591+372-9f3)

§4 Multistep Method 注：一般有 ，其中Bk 与yi+1 计算公式中 fi , …, fi-k 各项的系数均可查表得到 。 ( ) 2 ( 2) i k k i k R B h y  + + = 0 1 1 2 3 k 2 1 2 3 12 23 24 55 2 1 - 12 16 - 24 59 - 12 5 24 37 24 9 - 12 5 8 3 720 251 fi fi-1 fi-2 fi-3 … Bk … … … … … … … Misprint on p.204 常用的是 k = 3 的4阶亚当姆斯显式公式 (55 59 37 9 ) 24 i+1 = i + i - i-1 + i-2 - i-3 f f f f h y y

§4 Multistep Method a亚当姆斯隐式公式/ Adams implicit formulae 利用k+1个节点上的被积函数值f1,f,…,构造k阶 牛顿前插多项式。与显式多项式完全类似地可得到一系列 隐式公式,并有R2=B442y+2(mn),其中B与fFH,,…, k+1的系数亦可查表得到。 小于Bk 12 8 12 12 12 24 19 19 24 24 24 24 720 较同阶显 纔常用的是k=3的阶亚当姆斯隐式公式 式稳定 y社1=y+(9f++19-51+f2)

§4 Multistep Method  亚当姆斯隐式公式 /* Adams implicit formulae */ 利用k+1 个节点上的被积函数值 fi+1 , fi , …, fi-k+1构造 k 阶 牛顿前插多项式。与显式多项式完全类似地可得到一系列 隐式公式，并有 ，其中 与 fi+1 , fi , …, fi-k+1 的系数亦可查表得到。 ( ) 2 ( 2) i k k i k R B h y  + + = ~ Bk ~ 0 1 1 2 3 k 2 1 - 2 1 12 5 24 9 2 1 12 8 24 19 12 1 - 24 5 - 24 1 12 1 - 24 1 - 720 19 - fi+1 fi fi-1 fi-2 … Bk … … … … … … … ~ 常用的是 k = 3 的4阶亚当姆斯隐式公式 (9 19 5 ) 24 i+1 = i + i+1 + i - i-1 + i-2 f f f f h y y 小于Bk 较同阶显 式稳定

§4 Multistep Method 亚当姆斯预测校正系统 / Adams predictor-correp Predicted意:三步所用公 Spl:用Rnge-m法计、 value Pil AC的精度必须相同。 Se2:用dms显式计算预测值; 通常用经典 Runge- Kuta法配合4阶 Se3:用同阶 Adams隐式计算校正值 Adam公式。 251 4阶 adams显式公式的截断误差为y(x+)- Hey! look at the local truncation error of tRo(5r 4阶A隐成公武的截断糕为9=1的(m) 当h充分 时世峡以天发时: 251 Qy(+1)≈+1+(y 270 i+1 外推 fied final y(x+1)≈y+1 (v+1 extrapalategn+l 270 i+1

§4 Multistep Method  亚当姆斯预测-校正系统 /* Adams predictor-corrector system */ Step 1: 用Runge-Kutta 法计算前 k 个初值； Step 2: 用Adams 显式计算预测值； Step 3: 用同阶Adams 隐式计算校正值。 注意：三步所用公 式的精度必须相同。 通常用经典Runge￾Kutta 法配合4阶 Adams 公式。 Hey! Look at the local truncation error of the explicit and implicit Adams methods: and Don’t you think there’s something you can do? ( ) 720 251 5 (5) i h y  ( ) 720 19 5 (5) i 4阶Adams隐式公式的截断误差为 - h y  ( ) 720 19 ( ) 5 (5) i 1 i 1 i y x + - y + = - h y  4阶Adams显式公式的截断误差为 ( ) 720 251 ( ) 5 (5) i 1 i 1 i y x + - y + = h y  当 h 充分小时，可近似认为i  i，则： 19 251 ( ) ( ) 1 1 1 1  - - - + + + + i i i i y x y y x y ( ) 270 251 ( ) i+1  i+1 + i+1 - i+1 y x y y y ( ) 270 19 ( ) i+1  i+1 - i+1 - i+1 y x y y y Predicted value pi+1 Modified value mi+1 Corrected value ci+1 Modified final value yi+1 外推技术 /* extrapolation */

§4 Multistep method Adams 4th-Order predictor-corrector Algorithm To approximate thethe solution of the initial-value problem y=f(x,vacmtll-u At (N+l)equally s Input: endpoir 应为(cH1-p),但因cH1尚未 Output: approxm 算出,只好用(c-n;)取代之 Step I Set h=(b-a)7N, u Step 2 Fori=1, 2,3 Compute yi using classical Kutta method; Output(xi, yil? Step 3 Fori=4,.,n do steps 4-10 Step 5 pi1=y; +h(55f 1+37f1-2-9f3)/24;/ predict* Step 6 m 4 1=Pi+251(c, p, )/270; /*modify * Step 7 c=y, +h(9f(x u,m 41)+19f-5fi-+fi_)/24;/*correct*/ Step 8 Vii=Ci 1-19(Ci1-Pi1)/270;/*modify the final value*/ Step 9 output(x +1, v2+1); Step 10 forj=0, 1,2,3 Set xi=x+1; yi=V+1; /Prepare for next iteration Step 11 STOP

§4 Multistep Method Adams 4th-Order predictor-corrector Algorithm To approximate the the solution of the initial-value problem At (N+1) equally spaced numbers in the interval [a, b]. Input: endpoints a, b; integer N; initial value y0 . Output: approximation y at the (N+1) values of x. Step 1 Set h = (b - a) / N ; x0 = a; y0 = y0; Output ( x0 , y0 ); Step 2 For i = 1, 2, 3 Compute yi using classical Runge-Kutta method; Output ( xi , yi ); Step 3 For i = 4, …, N do steps 4-10 Step 5 ; /* predict */ Step 6 ; /* modify */ Step 7 ; /* correct */ Step 8 ; /* modify the final value */ Step 9 Output ( xi+1 , yi+1 ); Step 10 For j = 0, 1, 2, 3 Set xi = xi+1 ; yi = yi+1 ; /* Prepare for next iteration */ Step 11 STOP. 0 y'= f (x, y), a  x  b, y(a) = y pi+1 = yi + h(55 f i - 59 f i-1 + 37 f i-2 - 9 f i-3 )/ 24 mi+1 = pi+1 + 251(ci - pi )/ 270 (9 ( , ) 19 5 )/ 24 i+1 = i + i+1 i+1 + i - i-1 + i-2 c y h f x m f f f yi+1 = ci+1 -19(ci+1 - pi+1 )/ 270 应为( ci+1 - pi+1 ), 但因ci+1 尚未 算出，只好用( ci - pi )取代之

§4 Multistep method >基于泰勒展开的构造法 1=any+y1+…+yk+(B1fm1+Af+B、f+…+BfA) 将通式中的右端各项y1…,yk;f1,f1, ,fk分别在x;点作泰勒展开,与精确解 yx#1)在x;点的泰勒展开作比较。通过令同 类项系数相等,得到足以确定待定系数a ·· k B_1,B,…,B的等式,则可构造出线性多步法的 公式

§4 Multistep Method ➢ 基于泰勒展开的构造法 ... ( ... ) i 1 0 i 1 i 1 k i k 1 i 1 0 i 1 i 1 k i k y y y y h f f f f + - - - + - - =a +a + +a + b + b + b + + b 将通式中的右端各项 yi-1 , … , yi-k ; fi+1, fi-1 , … , fi-k 分别在 xi 点作泰勒展开，与精确解 y(xi+1) 在 xi 点的泰勒展开作比较。通过令同 类项系数相等，得到足以确定待定系数a0 , … , ak ; b -1 , b 0 , … , b k 的等式，则可构造出线性多步法的 公式

§4 Multistep Method 例:设pn=an1+aa1+a22+M(只y十y1+月y2+月y3 确定式中待定系数a,a,a2,A,B,月2,B,使得公式具有4阶 精度。 解:y-=y-hy+是h2y-是y+hy+O(h) Vi-2=y:-2hy'+2h2y-ih'y+2hty(+O(h) yi=y; -hy+2hym-6h'yi+o(h) yi_r=yi-2hy+2h y th'yi+o(h) ′3=y-3hy"+2h2ym-2h3y24)+O(h) y(xm1)=yi+hy+jhy+hv +1h'y()+o(h5)/*y(x )=yi* Ca+a1+a2=1 h(-a1-2a2+B+B1+B2)=h 7个未知数 h2(a1+2a2-B1-2B2-363)=lh h(-a1-a2+1月+2B2+263)=4h 5个方程 h(2a1+3a2-B1-月2-2B3)=是h

例：设 ( ) +1 0 1 -1 2 -2 0 1 -1 2 -2 3 -3 = + + + +  +  + i i i i i i i i y a y a y a y h b y b y b y b y 确定式中待定系数a0 , a1 , a2 , b0 , b1 , b2 , b3 , 使得公式具有4阶 精度。 §4 Multistep Method 解： ( ) 4 (4) 5 2 4 3 1 6 2 1 2 1 yi-1 = yi - hy i + h y i - h y i + h yi + O h 2 2 ( ) 4 (4) 5 3 3 2 3 2 4 yi-2 = yi - hy i + h y i - h y i + h yi + O h ( ) 3 (4) 4 6 2 1 2 1 y i-1 = y i - hy i + h y i - h yi + O h 2 2 ( ) 3 (4) 4 3 2 4 y i-2 = y i - hy i + h y i - h yi + O h 3 ( ) 3 (4) 4 2 2 9 2 9 y i-3 = y i - hy i + h y i - h yi + O h ( ) ( ) 4 (4) 5 2 4 3 1 6 2 1 2 1 y xi+1 = yi + hy i + h y i + h y i + h yi + O h /* y(xi ) = yi */ a0 +a1 +a2 = 1 h(-a1 - 2a2 + b 0 + b1 + b 2 ) = h 2 2 1 2 1 2 1 2 3 2 1 h ( a + 2a - b - 2b - 3b ) = h 3 6 1 2 3 9 2 1 2 1 3 2 4 6 1 3 1 h (- a - a + b + 2b + b ) = h 4 2 4 1 2 3 9 3 2 4 6 1 1 3 2 2 2 4 1 4 1 h ( a + a - b - b - b ) = h 个未知数 个方程 7 5

§4 Multistep Method 令a1=a2=0→Adm显式公式 以yh取代y1,并取a1=a2=0→ Adams隐式公式 以y3取代 则可导出-组4阶显式算法,其中 包含了著名的米尔尼pMn 取 得到辛甫生/ Simpson+公式 配使用 其”在区间x1,x1上积分,并用 i+1 Simpson数值积分公式来近似积 注:上分项,亦可得业mpn公式。区间 Ix13,xal上积分,侍到 PExi-3)+ f(x,y(x)dx 再过f,n,f2做f的插、试式即可。 辛甫生 Sim pson公式Pm=P1+3(m+4+y)

§4 Multistep Method  令 a1 = a2 = 0 Adams 显式公式  以 y i+1 取代 y i-1，并取a1 = a2 = 0 Adams 隐式公式  以 yi-3 取代 y i-3 ，则可导出另一组4 阶显式算法，其中 包含了著名的米尔尼 /* Milne */ 公式 (2 2 ) 3 4 +1 -3 -1 -2 = +  -  +  i i i i i y y y h y y 其局部截断误差为 ( ) , ( , ) 45 14 1 5 (5) Ri = h y  i  i  xi xi+ 注：上式也可通过数值积分导出，即将 在区间 上积分，得到 再过 做 f 的插值多项式即可。 y = f (x, y) [ , ] i-3 i+1 x x  + - + = - + 1 3 ( ) ( ) ( , ( )) , 1 3 i i x x y xi y xi f x y x dx 1 2 , , i i- i￾f f f 取 a1 = 1, a2 = 0 得到辛甫生 /* Simpson */ 公式 与Milne 公式匹配使用  辛甫生 /* Simpson */ 公式 ( 4 ) 3 +1 -1 +1 -1 = +  +  +  i i i i i y y y h y y 在区间[xi-1 , xi+1]上积分，并用 Simpson数值积分公式来近似积 分项，亦可得此Simpson公式

§4 Multistep method s Milne- Simpson系统的缺点是稳定性差,为改善稳定性, 考虑另一种隐式校正公式: VH1=any+ay21+a2y-2+(B1yn+A0y+A1y) 要求公式具有4阶料度。通过泰勒展开,可得到5个等式, 从中解出6个未知则有1个自由度 哈明/ Hamming用a不同鞫值进行试验,发现当a1=0 时,公式的稳定性较好m期公式 3 y#=(9y1-y=2)+(y+1+2y-y}1) 8 其局部截断误差为R 和y(9),5;∈(x1 注:哈明公式不能用数值积分方法推导出来

§4 Multistep Method  Milne-Simpson 系统的缺点是稳定性差，为改善稳定性， 考虑另一种隐式校正公式： ( ) +1 0 1 -1 2 -2 -1 +1 0 1 -1 = + + +  +  +  i i i i i i i y a y a y a y h b y b y b y 要求公式具有4 阶精度。通过泰勒展开，可得到 个等式， 从中解出 个未知数，则有 个自由度。 5 6 1 取 a1 = 1 得 Simpson 公式 哈明 /* Hamming */ 用a1 的不同数值进行试验，发现当a1 = 0 时，公式的稳定性较好，即： ( 2 ) 8 3 (9 ) 8 1 +1 -2 +1 -1 = - +  +  -  i i i i i i y y y h y y y 其局部截断误差为 ( ) , ( , ) 40 1 1 5 (5) i = - i i  xi xi+ R h y   注：哈明公式不能用数值积分方法推导出来