常微分方程初值问题的数值方法基础

1、工程数学工程数学工程数学工程数学1第十章第十章 常微分方程数值解常微分方程数值解第一节第一节 求解初值问题数值方法的基本原理求解初值问题数值方法的基本原理第二节第二节 高精度的单步法高精度的单步法 第三节第三节 线性多步法线性多步法第四节第四节 一阶微分方程组的解法一阶微分方程组的解法第五节第五节 边值问题的打靶法和差分法边值问题的打靶法和差分法工程数学工程数学工程数学工程数学2考虑一阶常微分方程的初值问题考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b R1 上连续，且关于

2、上连续，且关于 y 满足满足 Lipschitz 条条件件，即存在与，即存在与 x, y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述都成立，则上述IVP存存在唯一解。在唯一解。| ),(),(|2121yyLyxfyxf 要计算出解函数要计算出解函数 y(x) 在一系列节点在一系列节点 a = x0 x10,使得使得),(hyxhQyynnnn 111()()pO hp (, )nnQ xyh( , , )( , , )Q x y hQ x y hL yy对一切对一切 成立成立,则该方法收敛则该方法收敛,且有且有 y

3、y和和)(pnhOe 工程数学工程数学工程数学工程数学14由该定理可知整体截断误差总比局部截断误差低一阶由该定理可知整体截断误差总比局部截断误差低一阶 对改进的对改进的Euler法法, ),(,(),(),(yxhfyhxfyxfhyxQ 21于是有于是有 1( , , )( , , )( , )( , )2(,( , )(,( , )Q x y hQ x y hf x yf x yf xh yhf x yf xh yhf x y设设L为为f关于关于y的的Lipschitz常数常数,则由上式可得则由上式可得( , , )( , , )(1/2)Q x y hQ x y hLhyy限定限定h即可

4、知即可知Q满足满足Lipschitz条件条件,故而改进的故而改进的Euler法收敛法收敛.工程数学工程数学工程数学工程数学15例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6

5、.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 73. 稳定性稳定性工程数学工程数学工程数学工程数学16定义定义若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的 /*absolutely stable */。一般分析时为简单起

6、见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */yy 常数，可以常数，可以是复数是复数当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定， 的全体构成的全体构成绝对稳定区域绝对稳定区域。我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。000yy h h h工程数学工程数学工程数学工程数学17例：例：考察显式欧拉

7、法考察显式欧拉法110(1)nnnnyyhyhy 000yy 110(1)nnyhy 11110(1)nnnnyyh 由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：hh 1|1| h0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法11nnnyyh y 111nnyyh 11011nnh 可见绝对稳定区域为：可见绝对稳定区域为：1|1| h210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。法的好。工程数学工程数学工程数学工程数学1812()11()()()()()

8、( )2!( )( ,),( )( ,)( ,)(), ,nppnnnnnnxypTaylory xhhyy xy xhyxyxyxPyxf x yyxfx yfx y f x y 若若用用 阶阶多多项项式式近近似似函函数数有有：其其中中。但但由由于于公公式式中中各各阶阶偏偏导导数数计计算算复复杂杂，不不实实用用。第二节第二节 高精度的单步法高精度的单步法在高精度的单步法中在高精度的单步法中, ,应用最广泛的是应用最广泛的是Runge-KuttaRunge-Kutta( (龙格龙格- -库塔库塔) )方法方法一一、Runge-Kutta法的基本思想（法的基本思想（1）工程数学工程数学工程数学工

9、程数学19(0)(0)(0)(1)(1)(1)(2)(2)(2)( )(1); ;2,3,jjjjyffffyffxyffyffxyffyffjxy 一般地有一般地有工程数学工程数学工程数学工程数学20111112121 (,)11()22(,)(, )nnnnnnnnnnE ulerE uleryyhKE ulerKfxyyyhKKE ulerKfxyKfxhyhK 如如 果果 将将公公 式式 与与 改改 进进公公 式式 写写 成成 下下 列列 形形 式式 ：公公 式式 改改 进进公公 式式Runge-Kutta法的基本思想（法的基本思想（2）工程数学工程数学工程数学工程数学2111( ,

10、)()( , )( , ) nnf x yy xyf x yf x y以上两组公式都使用函数在某些点上的以上两组公式都使用函数在某些点上的值的线性组合来计算的近似值。值的线性组合来计算的近似值。EulerEuler公式：每步计算一次的值，为一阶方法。公式：每步计算一次的值，为一阶方法。改进Euler公式：需计算两次的值，二阶方法。改进Euler公式：需计算两次的值，二阶方法。工程数学工程数学工程数学工程数学22Runge-Kutta法的基本思想（法的基本思想（3）( ,)(,)( ) nnnf x yxyTaylory xxTaylor于于是是可可考考虑虑用用函函数数在在若若干干点点上上的的函

11、函数数值值的的线线性性组组合合来来构构造造近近似似公公式式，构构造造是是要要求求近近似似公公式式在在处处的的展展开开式式与与解解在在处处的的展展开开式式的的前前面面几几项项重重合合，从从而而使使近近似似公公式式达达到到所所需需要要的的阶阶数数。即即避避免免求求偏偏导导，又又提提高高了了方方法法的的精精度度，此此为为R RK K方方法法的的基基本本思思想想。11111(,)(,)(2, 3,)pnniiinniininijjjyyhc KKfxyKfxa h yhb Kip 工程数学工程数学工程数学工程数学23二、二阶龙格库塔方法二、二阶龙格库塔方法11111RK(,)(,)(2,3,),(,)

12、()pnniiinniininijjjiijinnnyyhc KKfxyKfxa h yhb KipabcxyTaylory xxTaylor 一一般般地地，方方法法设设近近似似公公式式为为其其中中，都都是是参参数数，确确定定它它们们的的原原则则是是使使近近似似公公式式在在处处的的展展开开式式与与在在处处的的展展开开式式的的前前面面项项尽尽可可能能多多地地重重合合。工程数学工程数学工程数学工程数学2411122122211()2(,)(,)nnnnnnyyh c Kc KKf xyKf xa h yhb K 当p时，近似公式为 当p时，近似公式为 112221122321(,)(,)(,(,)

13、(,) (,)(,)(,) (,)()nnnnnnnnnnnnnnnxnnynnnnxyTayloryyh c f xyc f xa h yhb f xyyh c f xycf xya hfxyhb fxyf xyO h 上式在处的展开式为上式在处的展开式为工程数学工程数学工程数学工程数学2512232221()(,)(,)(,)(,)()nnnxnnynnnnyccf xyhc a fxybfxyf xyhO h 123123()()()()()()2(,)(,)(,)(,)()2nnnnnnnnnxnnynnnny xxTaylorhy xy xhyxyxO hyfxyhhfxyfxyfx

14、yO h 在在处处的的展展开开式式为为工程数学工程数学工程数学工程数学2612222211 1 / 2 1 / 2 ),ccc ac bO 3 3有有 无无 穷穷 多多 组组 解解 ， 每每 一一 组组 解解 得得 一一近近 似似 公公 式式 ， 局局 部部 截截 断断 误误 差差 均均 为为( (h h这这 些些 方方 法法 统统 称称 二二 阶阶 方方 法法 。122211121211,1,2() / 2(,)(,)nnnnnnccabE uleryyh KKKfxyKfxhyhK 取取此此 为为 改改 进进公公 式式 。近近 似似 公公 式式 为为 122211212110,1,2(,)

15、(2,2)nnnnnnccabyyhKKfxyKfxhyhK 取取此此 为为 常常 用用 的的 二二 阶阶 公公 式式 ，称称 为为 中中 点点 公公 式式 。 工程数学工程数学工程数学工程数学27三、三阶龙格库塔方法三、三阶龙格库塔方法11231213123 (4)6(,)( ,)22(,2)nnnnnnnnpRKRKhyyKKKKf xyhhKf xyKKf xh yhKhK 类类似似地地，对对，即即三三个个点点，通通过过更更复复杂杂的的计计算算，可可导导出出三三阶阶公公式式。常常用用的的三三阶阶公公式式为为：工程数学工程数学工程数学工程数学28四、四阶龙格库塔方法四、四阶龙格库塔方法11

16、23412132434 (22)6(,)(,)2 2(,)22(,)nnnnnnnnnnpRKRKhyyKKKKKf xyhhKf xyKhhKf xyKKf xh yhK 对对，即即四四个个点点，可可导导出出四四阶阶公公式式。常常用用的的四四 阶阶公公式式为为： 工程数学工程数学工程数学工程数学29 h=0.2,x=0 x=12(01);(0)1.xyyxyy 设设取取步步长长从从直直到到用用四四阶阶龙龙格格库库塔塔方方法法求求解解初初值值问问题题例例：工程数学工程数学工程数学工程数学30112341211322433(22);62;2;222;222().nnnnnnnnnnnnnnhyy

17、KKKKxKyyxhhKyKhyKxhhKyKhyKxhKyh Kyh K 由由 经经 典典 的的 四四 阶阶 龙龙 格格 库库 塔塔解解公公 式式 得得：工程数学工程数学工程数学工程数学312RK RKRK RK）方方法法的的导导出出基基于于T Ta ay yl lo or r展展开开，故故要要求求所所求求问问题题的的解解具具有有较较高高的的光光滑滑度度。当当解解充充分分光光滑滑时时，四四阶阶方方法法确确实实优优于于改改进进E Eu ul le er r法法。对对一一般般实实际际问问题题，四四阶阶方方法法一一般般可可达达到到精精度度要要求求。如如果果解解的的光光滑滑性性差差，则则用用四四阶阶

18、方方法法解解的的效效果果不不如如改改进进E Eu ul le er r法法。两点说明两点说明:1RKRK46RK5）当当p=1,2,3,4p=1,2,3,4时时，公公式式的的最最高高阶阶数数恰恰好好是是p,p,当当p4p4时时，公公式式的的最最高高阶阶数数不不是是p p，如如p=5p=5时时仍仍为为 ，p=p= 时时公公式式的的最最高高阶阶数数为为 。工程数学工程数学工程数学工程数学32五、变步长的龙格五、变步长的龙格库塔方法库塔方法()1()51115(2)15(2)11(2)11()11,(),2,2()2,2()1.()16hnhnnnnhnhnnhnnhnnhyy xychhxxhyc

19、hy xycy xyy xy n n以以经经典典四四阶阶龙龙格格库库塔塔公公式式为为例例。从从节节点点x x 出出发发，以以 为为步步长长求求一一近近似似值值将将步步长长折折半半，即即取取为为步步长长从从跨跨两两步步到到，求求一一近近似似值值每每跨跨一一步步的的截截断断误误差差是是因因此此有有由由上上两两式式 (2)(2)()11111().15hhhnnnny xyyy 工程数学工程数学工程数学工程数学33R-K方法的绝对稳定区域/2121233223443 11 ,22111 )22411 )24 ( , ) () K()()()()()nnnnnnnRKhhhhf x yyhhyyyyh

20、KKyyhhKKhyyhhhKK 代入公式：代入公式：将将234112341 2)6111 12624(2()()()nnnhyyKKKKyhhh 工程数学工程数学工程数学工程数学34234111112624()()()nnhhhh 则则234111 112624()()()hhhh 绝绝对对稳稳定定区区域域： 2 1 -3 -2 -1 0 -1 -2 工程数学工程数学工程数学工程数学3504年：设年：设f(x)四阶连续可导，四阶连续可导，试建立如下数值微分公式试建立如下数值微分公式 并推导该公式的截断误差。并推导该公式的截断误差。0,0,1,2.ixxih i01212( ) 2 ( )( )( )f xf xf xf xh234(4)011111234(4)21111222(4)012128.(10) ()( )( )( )( )( )23!4!()( )( )( )( )()23!4!()2 ( )()( )( )12hhhf xf xhfxfxfxfhhhf xf xhfxfxfxfhf xf xf xhfxfh分两式相加除 ，得解：解：工程数学工程数学工程数学工程数学3605年：设年：设f(x)四阶连