第十三讲：微分方程数值解2

1、5 线性多步法线性多步法 /* multistep method */用用若干若干节点处的节点处的 y 及及 y 值的值的线性组合线性组合来近似来近似y(xi+1)。).(.110111101kikiiikikiiiffffhyyyy- - -+ +- - - -+ + + + + + + + + += =b bb bb bb ba aa aa a其通式可写为：其通式可写为：),(jjjyxff = =当当 b b- -1 0 时，为时，为隐式公隐式公式式; b b- -1=0 则为则为显式公式显式公式。 基于数值积分的构造法基于数值积分的构造法将将 在在 上积分，得到上积分，得到),(yxf

2、y =,1+iixx + += =- -+ +1)(,()()(1iixxiidxxyxfxyxy只要只要近似地算出右边的积分近似地算出右边的积分 ，则可通，则可通过过 近似近似y(xi+1) 。而。而选用不同近似式选用不同近似式 ik，可得到不，可得到不同的计算公式同的计算公式。 + + 1)(,(iixxkdxxyxfikiiiyy+ += =+ +15 multistep method 亚当姆斯显式公式亚当姆斯显式公式 /* adams explicit formulae */利用利用k+1 个节点上的被积函数值个节点上的被积函数值 构造构造 k 阶牛顿阶牛顿后插后插多项式多项式 , 有

3、有kiiifff- - -,.,11, 0, )( + +thtxnik + + + += =+ +1010)()()(,(1dthhtxrdthhtxndxxyxfikxxikiinewton插值余项插值余项dthtxnhyyikii)(101+ + += = + +/* 显式计算公式显式计算公式 */局部截断误差为：局部截断误差为： + += =- -= =+ + +1011)()(dthtxrhyxyrikiii例：例：k=1 时有时有)()(11- - -+ += = + += =+ +iiiiiifftfftfhtxn - - -+ +- -+ += =- -+ + += =1011

4、1)3(2)(iiiiiiiiffhydtfftfhyydththtdxyfdhrxxi)1(!21)(,(1022+ += = )(1253iyh = =5 multistep method注：注：一般有一般有 ，其中，其中bk 与与yi+1 计算公式计算公式中中 fi , , fi- -k 各项的各项的系数系数均可查表得到均可查表得到 。 )()2(2ikkkiyhbr + + += =10123k21231223245521-1216-2459-1252437249-12583720251fifi- -1fi- -2fi- -3bk常用的是常用的是 k = 3 的的4阶阿达姆斯显式公式阶

5、阿达姆斯显式公式)9375955(243211- - - -+ +- -+ +- -+ += =iiiiiiffffhyy5 multistep method阿达姆斯隐式公式阿达姆斯隐式公式 /* adams implicit formulae */利用利用k+1 个节点上的被积函数值个节点上的被积函数值 fi+1 , fi , , fi- -k+1 构造构造 k 阶阶牛顿牛顿前插前插多项式。与显式多项式完全类似地可得到一系列多项式。与显式多项式完全类似地可得到一系列隐式公式隐式公式，并有，并有 ，其中，其中 与与 fi+1 , fi , , fi- -k+1 的系数亦可查表得到。的系数亦可查

6、表得到。)()2(2ikkkiyhbr + + += =kb10123k21- -21125249211282419121- -245- -241121- -241- -72019- -fi+1fifi- -1fi- -2bk常用的是常用的是 k = 3 的的4阶阿达姆斯隐式公式阶阿达姆斯隐式公式)5199(242111- - -+ + + +- -+ + += =iiiiiiffffhyy小于小于bk较同阶显较同阶显式式稳定稳定5 multistep method阿达姆斯预估阿达姆斯预估-校正方法校正方法 /* adams predictor-corrector system */step

7、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 dont you think theres something

8、 you can do?)(720251)5(5iyh )(72019)5(5iyh - -4阶阶adams隐式公式的截断误差为隐式公式的截断误差为)(72019)()5(511iiiyhyxy - -= =- -+ + +4阶阶adams显式公式的截断误差为显式公式的截断误差为)(720251)()5(511iiiyhyxy = =- -+ + +当当 h 充分小充分小时，可近似认为时，可近似认为 i i ，则：，则：19251)()(1111- - - - -+ + + + +iiiiyxyyxy)(270251)(1111+ + + + +- -+ + iiiiyyyxy)(27019)

9、(1111+ + + + +- - - iiiiyyyxypredicted value pi+1modified value mi+1corrected value ci+1modified final value yi+1外推技术外推技术 /* extrapolation */5 multistep method adams 4th-order predictor-corrector algorithmto approximate the the solution of the initial-value problemat (n+1) equally spaced numbers in

10、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 s

11、teps 4-10step 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)(,),(yaybxayxfy= = = =24/ )9375955(3211- - - -+ +- -+

12、+- -+ += =iiiiiiffffhyp270/ )(25111iiiipcpm- -+ += =+ + +24/ )519),(9(21111- - -+ + + + +- -+ + += =iiiiiiifffmxfhyc270/ )(191111+ + + + +- - -= =iiiipccy应为应为( ci+1 - - pi+1 ), 但因但因ci+1 尚未尚未算出，只好用算出，只好用( ci - - pi )取代之。取代之。5 multistep method 基于泰勒展开的构造法基于泰勒展开的构造法).(.110111101kikiiikikiiiffffhyyyy- -

13、-+ +- - - -+ + + + + + + + + += =b bb bb bb ba aa aa a 将通式中的右端各项将通式中的右端各项 yi- -1, , yi- -k ; fi+1, fi- -1, , fi- -k 分别在分别在 xi 点作点作泰勒展开泰勒展开，与精确解，与精确解 y(xi+1) 在在 xi 点的泰勒展开作点的泰勒展开作比较比较。通过令。通过令同同类项系数相等类项系数相等，得到足以确定待定系数，得到足以确定待定系数a a0, , a ak ; b b- -1, b b0, , b bk 的等式，则可构造出线性多步法的的等式，则可构造出线性多步法的公式。公式。例：

14、例：设设)(3322110221101- - - - - -+ + + + + + + += =iiiiiiiiyyyyhyyyyb bb bb bb ba aa aa a确定式中待定系数确定式中待定系数a a0, a a1, a a2, b b0, b b1, b b2, b b3, 使得公式具有使得公式具有4阶阶精度。精度。5 multistep method解：解：)(5)4(42413612211hoyhyhyhyhyyiiiiii+ + + - - + + - -= =- -)(225)4(43233422hoyhyhyhyhyyiiiiii+ + + - - + + - -= =-

15、 -)(4)4(3612211hoyhyhyhyyiiiii+ +- - + + - - = = - -)(224)4(33422hoyhyhyhyyiiiii+ +- - + + - - = = - -)(34)4(3292293hoyhyhyhyyiiiii+ +- - + + - - = = - -)()(5)4(42413612211hoyhyhyhyhyxyiiiiii+ + + + + + + + += =+ +/* y(xi) = yi */1210= =+ + +a aa aa ahh= =+ + + +- - -)2(21021b bb bb ba aa a221321212

16、12)322(hh= =- - - -+ +b bb bb ba aa a36132921212341613)2(hh= =+ + + +- - -b bb bb ba aa a424132923416123212414)(hh= =- - - -+ +b bb bb ba aa a个未知数个未知数个方程个方程755 multistep method 令令 a a1 = a a2 = 0adams 显式显式公式公式 以以 y i+1 取代取代 y i- -1，并取，并取 a a1 = a a2 = 0adams 隐式隐式公式公式 以以 yi- -3 取代取代 y i- -3 ，则可导出另一组，

17、则可导出另一组4 阶显式算法，其中阶显式算法，其中包含了著名的包含了著名的米尔恩米尔恩 /* milne */ 公式公式)22(342131- - - -+ + + + - - + += =iiiiiyyyhyy其局部截断误差为其局部截断误差为),(,)(45141)5(5+ + = =iiiiixxyhr 注：注：上式也可通过上式也可通过数值积分数值积分导出，即将导出，即将 在区间在区间 上积分，得到上积分，得到 再过再过 做做 f 的插值多项式即可。的插值多项式即可。),(yxfy =,13+-iixx + +- -+ += =- -+ +13,)(,()()(31iixxiidxxyxf

18、xyxy21,- - -iiifff取取 a a1 = 1, a a2 = 0得到得到辛甫生辛甫生 /* simpson */ 公式公式与与milne 公式匹配使用公式匹配使用 辛甫生辛甫生 /* simpson */ 公式公式)4(31111- -+ +- -+ + + + + + + += =iiiiiyyyhyy在区间在区间xi- -1, xi+1上积分，并用上积分，并用simpson数值积分数值积分公式来近似积公式来近似积分项，亦可得此分项，亦可得此simpson公式。公式。4 multistep method milne-simpson 方法的缺点是方法的缺点是稳定性差稳定性差，为改善稳定性，为改善稳定性，考虑另一种隐式校正公式：考虑另一种隐式校正公式：)(11011221101-+-+=iiiiiiiyyyhyyyyb bb bb ba aa aa a要求公式具有要求公式具有4 阶精度。通过泰勒展开，可得到阶精度。通过泰勒展开，可得到 个等式，个等式，从中解出从中解出