NA005b微分方程求解

华中科技大学 李红 §4 线性多步法 /* Multistep Method */ 用若干节点处的 y及 y’值的线性组合来近似 y(xn+1)。 )...(... 110111101 knknnnknknnn ffffhyyyy −−+−−−+ ++++++++= ββββααα 其通式可写为： ¾基于数值积分的构造法 将 在 上积分，得到),( yxfy =′ ],[ 1+nn xx ∫ +=−+ 1 ))(,()()( 1 n n x xnn dxxyxfxyxy 只要近似地算出右边的积分 ，则可通 过 近似y(xn+1)。而选用不同近似式 In,可得到不 同的计算。 ∫ +≈ 1 ))(,(n n x xn dxxyxfI nnn Iyy +=+1 华中科技大学 李红 §4 Multistep Method # 亚当姆斯显式公式 /* Adams explicit formulae */ 利用k+1个节点上的被积函数值 构造 k阶牛顿 后插多项式 , 有 knnn fff −− ,...,, 1 ]1,0[,)( ∈+ thtxN nk ∫∫ ∫ +++=+ 1010 )()())(,(1 dthhtxRdthhtxNdxxyxf nkxx nknn dthtxNhyy nknn )( 1 01 ++= ∫+ /* 显式计算公式 */ 局部截断误差为： ∫ +=−= +++ 10111 )()( dthtxRhyxyT nknnn 例：k=1时有 )()( 11 −−+=∇+=+ nnnnnn fftfftfhtxN ∫ −−+ −+=−++= 10 111 )3(2)]([ nnnnnnnn ffhydtfftfhyy dththt dx yfdhT xxn )1(!2 1))(,(1 0 2 2 1 += ∫+ ξξ )(125 3 nyh ξ′′′= 华中科技大学 李红 §4 Multistep Method 注：一般有 ，其中Bk与yn+1 计算公式 中 fn , …, fn−k各项的系数均可查表得到。 )()2(21 n kk kn yhBT ξ+++ = 10 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 fn fn−1 fn−2 fn−3 … Bk … … … … … … … 常用的是 k = 3的4阶亚当姆斯显式公式 )9375955( 24 3211 −−−+ −+−+= nnnnnn ffffhyy 华中科技大学 李红 §4 Multistep Method # 亚当姆斯隐式公式 /* Adams implicit formulae */ 利用k+1个节点上的被积函数值 fn+1 , fn , …, fn−k+1构造 k阶 牛顿前插多项式。与显式多项式完全类似地可得到一系列 隐式公式，并有 ，其中 与 fn+1 , fn , …, fn−k+1的系数亦可查表得到。 )()2(21 n kk kn yhBT η+++ = ~ kB~ 10 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− fn+1 fn fn−1 fn−2 … Bk … … … … … … … ~ 常用的是 k = 3的4阶亚当姆斯隐式公式 )5199( 24 2111 −−++ +−++= nnnnnn ffffhyy 小于Bk 华中科技大学 李红 §4 Multistep Method ¾基于泰勒展开的构造法 )...(... 110111101 knknnnknknnn ffffhyyyy −−+−−−+ ++++++++= ββββααα 将通式中的右端各项 yn−1, … , yn−k ; fn+1, fn−1, … , fn−k分别在 xn点作泰勒展开，与精确解 y(xn+1) 在 xn点的泰勒展开作比较。通过令 同类项系数相等，得到足以确定待定系数α0, … , αk ; β−1, β0, … , βk的等式，则可构造出线性多步法 的公式。 华中科技大学 李红 例：设 )( 3322110221101 −−−−−+ ′+′+′+′+++= nnnnnnnn yyyyhyyyy ββββααα 确定式中待定系数α0, α1, α2, β0, β1, β2, β3,使得公式具有4阶 精度。 §4 Multistep Method 解： )( 5)4(42413612211 hOyhyhyhyhyy nnnnnn ++′′′−′′+′−=− )(22 5)4(432 3 3 42 2 hOyhyhyhyhyy nnnnnn ++′′′−′′+′−=− )( 4)4(361 2 2 1 1 hOyhyhyhyy nnnnn +−′′′+′′−′=′ − )(22 4)4(334 2 2 hOyhyhyhyy nnnnn +−′′′+′′−′=′ − )(3 4)4(329 2 2 9 3 hOyhyhyhyy nnnnn +−′′′+′′−′=′ − )()( 5)4(4241 3 6 12 2 1 1 hOyhyhyhyhyxy nnnnnn ++′′′+′′+′+=+ /* y(xn) = yn */ 1210 =++ ααα hh =++++−− )2( 321021 ββββαα 2 2 1 321212 12 )322( hh =−−−+ βββαα 3 6 1 32 9 212 1 23 4 16 13 )2( hh =+++−− βββαα 4 24 1 32 9 23 4 16 1 23 2 124 14 )( hh =−−−+ βββαα 个未知数 个方程 7 5 华中科技大学 李红 §4 Multistep Method ) 令 α1 = α2 = 0 Adams显式公式 ) 以 y′n+1取代 y′n−3，并取 α1 = α2 = 0 Adams隐式公式 ) 以 yn−3 取代 y′n−3，则可导出另一组4阶显式算法，其中 包含了著名的米尔尼 /* Milne */公式 )22( 3 4 2131 −−−+ ′+′−′+= nnnnn yyyhyy 其局部截断误差为 ),(,)( 45 14 1 )5(5 1 ++ ∈= nnnnn xxyhT ξξ 注：上式也可通过数值积分导出，即将 在区间 上积分，得到 再过 做 f的插值多项式即可。 ),( yxfy =′ ],[ 13 +− nn xx ∫ +−+= −+ 13 ,))(,()()( 31 nnxxnn dxxyxfxyxy 21 ,, −− nnn fff ) 辛甫生 /* Simpson */公式 )4(3 1111 −+−+ ′+′+′+= nnnnn yyy hyy 华中科技大学 李红 §4 Multistep Method ) Milne-Simpson 系统的缺点是稳定性差，为改善稳定性， 考虑另一种隐式校正公式： )( 11011221101 −+−−−+ ′+′+′+++= nnnnnnn yyyhyyyy βββααα 要求公式具有4 阶精度。通过泰勒展开，可得到 个等 式，从中解出 个未知数，则有 个自由度。 5 6 1 哈明 /* Hamming */用α1 的不同数值进行试验，发现当α1 = 0 时，公式的稳定性较好，即： )2( 8 3)9( 8 1 1121 −+−+ ′−′+′+−= nnnnnn yyyhyyy 其局部截断误差为 ),(,)( 40 1 1 )5(5 1 ++ ∈−= nnnnn xxyhT ξξ 注：哈明公式不能用数值积分方法推导出来。 华中科技大学 李红 §4 Multistep Method # 亚当姆斯预测-校正系统 /* Adams predictor-corrector system */ Step 1:用Runge-Kutta法计算前 k个初值； Step 2:用Adams 显式计算预测值； Step 3:用同阶Adams 隐式计算校正值。 注意：三步所用公 式的精度必须相 同。通常用经典 Runge-Kutta 法配 合4阶Adams 公 式。 4阶Adams隐式公式的截断误差为 )( 720 19)( )5(511 nnn yhyxy η−=− ++ 4阶Adams显式公式的截断误差为 )(720 251)( )5(511 nnn yhyxy ξ=− ++ 当 h充分小时，可近似认为ξn ≈ ηn，则： 19 251 )( )( 11 11 −≈− − ++ ++ nn nn yxy yxy )( 270 251)( 1111 ++++ −+≈ nnnn yyyxy )( 270 19)( 1111 ++++ −−≈ nnnn yyyxy 华中科技大学 李红 §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 n = 1, 2, 3 Compute yn using classical Runge-Kutta method; Output ( xn , yn ); Step 3 For n = 4, …, N do steps 4-10 Step 5 ; /* predict */ Step 6 ; /* modify */ Step 7 ; /* correct */ Step 8 ; /* modify the final value */ Step 9 Output ( xn+1 , yn+1 ); Step 10 For j = 0, 1, 2, 3 Set xn = xn+1 ; yn = yn+1 ; /* Prepare for next iteration */ Step 11 STOP. 0)(,),,(' yaybxayxfy =≤≤= 24/)9375955( 3211 −−−+ −+−+= nnnnnn ffffhyp 270/)(25111 nnnn pcpm −+= ++ 24/)519),(9( 21111 −−+++ +−++= nnnnnnn fffmxfhyc 270/)(19 1111 ++++ −−= nnnn pccy HW: p. 148-149 #9, 10