—常微分方程的数值解法.ppt

《—常微分方程的数值解法.ppt》由会员分享，可在线阅读，更多相关《—常微分方程的数值解法.ppt（63页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、本章内容本章内容n6.1 引言引言n6.2 欧拉方法欧拉方法n6.3 龙格龙格库塔方法库塔方法n6.4 边值问题的数值方法边值问题的数值方法n一一. 问题提出问题提出 有一个或多个导数及其函数的方程式称为微分有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。方程，在工程中常遇到求解微分方程的问题。6.1 引言引言)()(,),(dd00 xyyyxybaxyxfxy解函数初值问题如，一阶常微分方程的6.1 引言引言n二二. 两类定解问题两类定解问题常微分方程的定解问题有两种基本类型类常微分方程的定解问题有两种基本类型类:初值问题初值问题和和边值问题边值问题定解指已

2、知因变量和定解指已知因变量和/或其导数在某些点上是已知的或其导数在某些点上是已知的(约束条件约束条件)。1. 边值问题边值问题 约束条件为已知，在自变量的任一非初值上，已知约束条件为已知，在自变量的任一非初值上，已知函数值和函数值和/或其导数值，如或其导数值，如 例如例如,受连续分布横向荷载的变截面简支梁弯曲问题受连续分布横向荷载的变截面简支梁弯曲问题ybyaybaxyyxfy求解 )(,)(, ),(0)()0()()(dd22lwwxEIxMxw简支梁跨度弯曲刚度挠度，其中，lxEIxw,)()(q(x)xwOl2. 初值问题初值问题00)(),(yxyyxfy例如，单自由度系统的非线性受

3、迫振动例如，单自由度系统的非线性受迫振动 0000)(,)(),(yxyyxyyyxfyxm( )P t023200dd( )ddd(0)dx xxxmckxxP tttxxxxt，唯一。初值问题的解必存在且由常微分方程理论知：，使得即存在常数条件，满足李普希兹且关于适当光滑连续，只要函数yyLyxfyxfLLipschitzyyxf),(),()(),(实际问题中还存在初边值混合问题，如梁在横向激励下的弯曲振动。高阶常微分方程可以化成一阶的常微分方程组00)(),(yxyyxfy 很多微分方程的解不能用初等函数来表示，有时很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但

4、计算量太大而不实即使能够用解析式表示其解，但计算量太大而不实用用(表达式过于复杂表达式过于复杂)。 需要用数值方法来求解，一般只要求得到若干个需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式点上的近似值或者解的简单的近似表达式(精度要精度要求满足即可求满足即可)。计算）数值解（适合于计算机近似解析解，求近似解求精确解一般较为困难6.1 引言引言 ), 2 , 1 , 0(), 2 , 1 , 0()(, :1)(),(011000niihxxxxhniyxybxxxayxyyxfyiiiiin，有通常取等距节点，即。的近似值得到在数值解式6.1 引言引言yx( )

5、yf x 0y0 x0hxxxxx)()()()()()(1211211210nnnnnnnnxyyyyyyxyxyxyxyxxxxxxy的近似解上的值在一些离散点即解函数6.1 引言引言初值初值问题问题的的常见常见解法解法单步法：单步法： 利用前一个单步的信息利用前一个单步的信息(一个点一个点)，在，在y=f(x) 上找下一点上找下一点yi， 有欧拉法，龙格库格法。有欧拉法，龙格库格法。预测校正法：预测校正法： 多步法，利用一个以上的前点信息求多步法，利用一个以上的前点信息求f(x)上的下一个上的下一个yi， 常用迭代法，如改进欧拉法，阿当姆斯法。常用迭代法，如改进欧拉法，阿当姆斯法。6.1

6、 引言引言得到。，这可以通过递推公式，求出yy,y,知的y进，由已点：按节点顺序依次推初值问题的数值解法特1ii106.2 欧拉方法及其改进欧拉方法及其改进 Eulers Method n内容内容一一. 欧拉格式欧拉格式二二. Euler预估预估校正法校正法三三. 误差估计、收敛性和稳定性误差估计、收敛性和稳定性 6.2.1 欧拉公式：欧拉公式：/* Eulers Method */向前差商近似导数向前差商近似导数hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y记为记为)1,., 0(),(1 niyxfhyyiiii6.2 欧拉方法及其改进欧拉方

7、法及其改进1111111( )(,)(,)-(,),(,)-()nnnnnnnnnnnnnnnnnnnyy xPP xyPxyPPyyf xyyyhf xyxxEuler几何意义：折线逼近解曲线。设已做出折线的顶点 ，过依方向场的方向再推进，显然两个顶点 ，的坐标有关系即这就是欧拉公式。xP0P1P2P3P4PnyO( )yy x，计计算算结结果果如如下下表表：分分别别取取步步长长公公式式为为解解：此此时时的的。确确解解是是的的数数值值解解，此此问问题题的的精精，方方法法求求初初值值问问题题例例：利利用用05. 0 , 1 . 0 , 2 . 02 , 1 , 0,0)211()1/()(0)

8、0(202110221222 hiyyxhyyEulerxxxyyxyxyEuleriiii6.2 欧拉方法及其改进欧拉方法及其改进7.2 欧拉方法欧拉方法hxiyi真值真值y(xi)误差误差y(xi)-yih=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.00000

9、0.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227越大越好。阶方法则称该方法为。

10、误差为若某种方法的局部截断：定义的误差，记为步严格算出逐，由去掉假设）整体截断误差（这一步的误差，记为严格算出，按某方法由假设）局部截断误差（：定义计数值解法的截断误差估pphORyxyyyxyyyxyRyyxyyPiiiiiiiiiiiiii)(2)()(2)()(11 111111011116.2.26.2 欧拉方法及其改进欧拉方法及其改进计算方法09计11、617.2 欧拉方法欧拉方法，一阶方法相减得与式式显格式为：，则假设积分法）展开法（也可使用数值解：使用显格式的局部截断误差例：求)()(2)()()()()(,()(),()()()(2)(,()(2)()()()(21211112

11、21hOyhyxyRbabxyxhfxyyxhfyyEulerxyyayhxyxhfxyhyhxyxyxyTaylorEulerniiiiiiiiiiiiiiiiiiii 6.2 欧拉方法及其改进欧拉方法及其改进计算方法09计11、617.2 欧拉方法欧拉方法（具有普遍意义）低一阶。比局部截断误差，则若初始误差为（证略）为求解范围。，（则足估计式方法的整体截断误差满有界，则且条件有关于设：定理系整体”截断误差的关“局部iibxaababLiRbaxyMieLhMyxyeEulerxyyyLyxfyxfLipschitzyyxf0,)(max), 2 , 112)( )(),(),(),(102

12、)(200)(212211 6.2 欧拉方法及其改进欧拉方法及其改进6.2.3 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商来向后差商来近似导数近似导数hxyxyxykkk)()()(11)1,., 0(),(111 niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得同时出现在等式的两边，不能直接得到，故称为到，故称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再用隐式法（一般先用显式计算一个初值

13、，再用隐式法（迭代）迭代）求解。求解。)(,)(111kkkxyxfxy)211( 2 , 1 , 0,0)211(0)0(202112121)0(1)0(1102121122iiiiiiiiiiyxhyyEuleryyiyyxhyyEuleryxyxy公式得到可由显式，迭代初值迭代求解公式为此时的隐式，例：求初值问题6.2 欧拉方法及其改进欧拉方法及其改进1221111221111221111( )()( )()()()()2()()( )(, ()()( )2()(,)()( )2( )( )iiiiiiiiiiiiiiiiiiiiiiiiy xxTayloryy xy xy xxxxxy

14、y xy xf xy xhxxayyyf xyhxxbabEulerR将在作展开即又由式和得隐式法的截断误差为22111()()2iiiiyy xyh 6.2 欧拉方法及其改进欧拉方法及其改进显格式确定。，初值用、用迭代法求情况）、化为显格式（仅小数是隐格式。实用上（高一阶期望截断误差的阶数会若将两者做平均，可以正负号法的截断误差相差一个法与隐式显式EulerynyxfyxfhyyEulerEuleriiiiiii1111 21 ) 1,0,1,2,(i ),),2显、隐式两种算法的显、隐式两种算法的平均。需要迭代求解，能否不迭代？平均。需要迭代求解，能否不迭代？6.2 欧拉方法及其改进欧拉方

15、法及其改进6.2.4 梯形格式梯形格式 6.2 欧拉方法及其改进欧拉方法及其改进校正法预估Euler6.2.5/* predictor-corrector method */Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii 1,2 , 1 , 0,),(),()(20121211niyhKyhxfKyxf

16、KKKhyyEuleriiiiii 方法也可以写成：方法也可以写成：改进的改进的7.2 欧拉方法及其改进欧拉方法及其改进7.2 欧拉方法欧拉方法 9 ,2 , 1 , 0,11 . 0)1 . 0(21 . 0/2)(05. 0)2(9 ,2 , 1 , 0,1/2 . 01 . 1)1()21()(1 . 01)0(10201121211012/1iyKyxKyKyxyKKKyyEuleriyyxyyEulerxxyhyxyxyyiiiiiiiiiiii方方法法利利用用改改进进方方法法利利用用解解：）。（精精确确解解为为的的数数值值解解，取取步步长长，例例：求求初初值值问问题题7.2 欧拉方

17、法欧拉方法ixiEuler方法方法yi改进改进Euler法法yi精确解精确解y(xi)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.

18、6124521.6733201.732051计算结果如下计算结果如下:6.3 龙格龙格 库塔方法库塔方法n内容内容一一. 2阶龙格阶龙格 库塔格式库塔格式三三. 高阶龙格高阶龙格 库塔格式库塔格式单步法：单步法：即利用前一个节点的函数值即利用前一个节点的函数值yi，计算后，计算后一个节点的函数值一个节点的函数值yi+1。目的：目的：建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜率某一斜率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其点。欧拉法及其各种变形所能达到的最高

19、精度为各种变形所能达到的最高精度为2 阶。阶。6.3 龙格龙格 库塔方法库塔方法n二二. 2阶龙格阶龙格 库塔格式库塔格式斜率斜率一定取一定取K1 K2的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？ )(),(),(2121),(),(2),(00121211111111xyyhKyhxfKyxfKKKhyyyxfyxfhyyyxhfyyiiiiiikkkkkkkkkk可可将将其其改改写写为为：考考察察改改进进的的欧欧拉拉法法，6.3 龙格龙格 库塔方法库塔方法7.3 龙格龙格 库塔方法库塔方法首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有，使得到的算

20、法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 )(iixyy )()(311hOyxyRiii Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 3: 将将

21、 y( xi+1 ) 在在 xi 点的点的泰勒泰勒展开并与展开并与yi+1作比较作比较)()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，则必须有：，则必须有：)()(311hOyxyRiii 21,1221 p 这里有这里有3个未知个未知数，数，2个方程。个方程。存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙阶龙格格 - 库塔格式库塔格式。21, 121 p注意到，注意到， 就是就是改进的欧拉法改进的欧拉法。 Step 2: 将将 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyph

22、xyhyhOxyphxyxyhyyiiiiiiii 其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待定系数，确定这些系数的步骤均为待定系数，确定这些系数的步骤与前面相似。与前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 问题问题: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？6.3 龙格龙格 库塔方

23、法库塔方法n三三. 高阶龙格高阶龙格 库塔格式库塔格式3阶阶龙格龙格-库塔法库塔法)2,(),(),()4(2131222132161hKhKyhxfKKyxfKyxfKKKKyyiihihiiihii 6.3 龙格龙格 库塔方法库塔方法6.3 龙格龙格 库塔方法库塔方法最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法),(),(),(),()22(34222312221432161hKyhxfKKyxfKKyxfKyxfKKKKKyyiihihihihiiihii 注：注： 龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值

24、。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数)(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2 nhO8 n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。6.3 龙格龙格 库塔方法库塔方法计算方法09计11、617.3 龙格龙格 库塔方

25、法库塔方法例：使用高阶例：使用高阶R-K方法计算初值问题方法计算初值问题解解: （1）使用三阶）使用三阶R-K方法方法(公式公式) 1)0(5 . 002yxyy. 1 . 0 h取取1111. 1)4(61 . 02555. 1)2( 1 . 0(1025. 1)21 . 0(1132101212032102201 KKKyyKKyKKyKyKn时时6.3 龙格龙格 库塔方法库塔方法计算方法09计11、617.3 龙格龙格 库塔方法库塔方法其余结果如下其余结果如下: i xi K1 K2 K3 yi 1.0000 0.1000 1.0000 1.1025 1.2555 1.1111 2.00

26、00 0.2000 1.2345 1.3755 1.5945 1.2499 3.0000 0.3000 1.5624 1.7637 2.0922 1.4284 4.0000 0.4000 2.0404 2.3423 2.8658 1.6664 5.0000 0.5000 2.7768 3.2587 4.1634 1.9993 6.3 龙格龙格 库塔方法库塔方法计算方法09计11、617.3 龙格龙格 库塔方法库塔方法（2）如果使用四阶）如果使用四阶R-K方法方法(公式公式)1111. 1)22(61 . 02351. 1)1 . 0(1133. 1)21 . 0(1025. 1)21 . 0(

27、11432101230422032102201 KKKKyyKyKKyKKyKyKn时时6.3 龙格龙格 库塔方法库塔方法计算方法09计11、617.3 龙格龙格 库塔方法库塔方法其余结果如下其余结果如下: i xi K1 K2 K3 K4 yi 1.0000 0.1000 1.0000 1.1025 1.1133 1.2351 1.1111 2.0000 0.2000 1.2346 1.3756 1.3921 1.5633 1.2500 3.0000 0.3000 1.5625 1.7639 1.7908 2.0423 1.4286 4.0000 0.4000 2.0408 2.3428 2

28、.3892 2.7805 1.6667 5.0000 0.5000 2.7777 3.2600 3.3476 4.0057 2.00006.3 龙格龙格 库塔方法库塔方法 步长过大，达不到精度要求；步长过小，虽然局部截断误差小，加大了计算工作量，舍入误差的累积增大。解决途径之一引入变步长技术，常用的有Richardson外推法。从结点 xi出发，先以h为步长，通过一步计算出y(xi+1)的近似值6.4 步长的自动选择步长的自动选择( )1hiy( )1211()()(1 )ihppiiy xychO hchxh在一般情况下，系数既依赖于 ，又依赖于 ，但当 较小时，可以近似地看成常数。1( /

29、2)( /21)11112/ 2 ,()2()(,22)2ihiihpppiihy xycO hhxh然后，将步长折半，即以为步长，仍从出发经过两部计算求得y(x ）的另一个近似值y其中每一步的截断误差约为c故( /2)( )2( /2)( )211111121) ()2()(2(2)(1)2()()=(4)21)21)3)phhpiiipppphiiiphpyy xyyOyO hy xh以乘以式的两端，并与式相减，得到即 （( /2)( )11i/2)1(11+)12(24)()1)phhiihhpiiiyyyxyyy （式表明，若取 作为的近似值，精度比和5）（要高一阶。( /2( /2)

30、( /2)( )111)1( /211)(11)(5)()2)=1(1)hhhhihhiiiiiiipyyy xxyyyyy式整理成 若将作为的近似值，则其误差可以用前后两次步长的计算结果的差来表示，即可用来判断所取步长是（6）（否恰当。( /2)( )11( )=hhiiixy xyy 在当前结点 ，若，反复将步长折半进行计算，直至精度满足。若，反复将步长加倍，直至精度超出。上述方法，表面上看，虽然增加了计算量，但从总体考虑常常是合算的，尤其是方程的解变化剧烈的情况下。计算方法09计11、617.2 欧拉方法欧拉方法( )yy x 0yx0 nxx0y，则称该方法收敛。）时若有（同时当对某方

31、法，任意固定定义：收敛性0)(0.0iiiiyxyihihxxx1/* Convergency */6.5 收敛性与稳定性收敛性与稳定性均均收收敛敛）（当当校校法法：预预梯梯形形法法，）（当当显显式式方方法法： 00)(00)(2hhOEulerhhOEulernn 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。0)0()()(yyxyxy解：解：该问题的精确解为该问题的精确解为 xeyxy 0)( 欧拉公式为欧拉公式为iiiiyhyhyy)1 (1 0)1 (yhyii 对任意固定的对任意固定的 x = xi = i h ，有，有xhhxihyhyy)1(

32、)1 (/10/0ehhh /10)1(lim)(0ixxyeyi 6.5 收敛性与稳定性收敛性与稳定性计算方法09计11、617.2 欧拉方法欧拉方法例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。1)0()(30)(yxyxy0.00.10.20.30.40.5精确解改进欧拉法 欧拉隐式欧拉显式 节点 xixey30 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56

33、251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107原因原因 ?!稳定性2.取步长取步长 h = 0.16.5 收敛性与稳定性收敛性与稳定性计算方法09计11、617.2 欧拉方法欧拉方法稳定性比较复杂。无限扩大。不随绝对稳定则称该方法关于步长）（）。若（实际为步计算值有误差引起第），则（实际初值为现假设初值有微小误差，严格得出作计算，由初值用某方法固定步长ihiSyiyyyhiiiiii, 2 , 10000

34、0:稳定性的定义6.5 收敛性与稳定性收敛性与稳定性法是不稳定的法及其改进的时，当例如，初值问题EulerEulerhyxyxy1 . 01)0()(30)(一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程yy 常数，可常数，可以是复数以是复数hhhh若试验方程的数值解关于步长 绝对稳定，一般来说，满足绝对稳定的 有许多，即存在某一范围。则将在复平面上的相应范围称为绝对稳定区域。6.5 收敛性与稳定性收敛性与稳定性方法绝对稳定。时，即方法绝对稳定显然，与上式相减得，则设初值有小扰动法例：对EulerhhEulerhyhyyhyhyhyyEuleriiiiiiiiii11

35、1)1 ()()1 ()1 ()1 ()(0110011100116.5 收敛性与稳定性收敛性与稳定性RehImh( 1,0)h复平面Euler法的绝对稳定区域定。，隐式的梯形法绝对稳可以证明：对任意步长h法是稳定的时，当初值问题Eulerhyxyxy301,1)0()(30)(7.3 龙格龙格 库塔方法库塔方法例：例：隐式龙格隐式龙格-库塔法库塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为 )2,2(1111KhyhxfKhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域

36、为的绝对稳定区域为0ReImk=1k=2k=3k=4-1-2-3-123ReIm无条件稳定无条件稳定6.6 一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程 0000)(),()(),(zxzzyxgzyxyzyxfy 可以把单个方程 中的 f 和 y 看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。 ),(yxfy 6.6.1 一阶常微分方程组对于一阶常微分方程组的初值问题对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。程可转化为一阶常微分方程组来求解。 设

37、为节点上的近似解，则有改进的Euler格式为 iiizyiihxx,);, 3 , 2 , 1(0),(1iiiiizyxhfyy),(1iiiiizyxhgzz预报：),(),(21111iiiiiiiizyxfzyxfhyy),(),(21111iiiiiiiizyxgzyxghzz校正：校正： 又，相应的四阶龙格库塔格式（经典格式）为 )22(6)22(64321143211LLLLhzzKKKKhyyiiii),(),()2,2,()2,2,()2,2,()2,2,(),(),(331433142221322213112121121211hLzhKyxgLhLzhKyxfKLhzKhy

38、xgLLhzKhyxfKLhzKhyxgLLhzKhyxfKzyxgLzyxfKiiiiiiiiiiiiiiiiiiiiiiii式中式中 例 用改进的Euler法求解初值问题 2)0(1)0(zzyxzyzxyy2.00 x取步长h=0.1，保留六位小数。 解: 改进的Euler法公式为111111105. 0)()(05. 0iiiiiiiiiiiiiiiizyxzyxzzzyxzyxyyiiiiiiiiiizyxzzzyxyy1.0)(1.011由初值由初值 , ,计算得计算得 2)0(, 1)0(00zzyy050000.2800000.011zy046951. 2) 1 . 0(801

39、500. 0) 1 . 0(11zzyy090992.2604820.022zy088216. 2)2 . 0(604659. 0)2 . 0(22zzyy6.6.1 一阶常微分方程组 高阶微分方程(或方程组)的初值问题,原则上都可以归结为一阶方程组来求解。例如,有二阶微分方程的初值问题 0000)(,)(),(yxyyxyyyxfy在引入新的变量 后,即化为一阶方程组初值问题:0000)(,)(,),(yxzyxyzyzyxfz此为一个一阶方程组的初值问题，对此可用前面中介绍的方法来求解。例如应用四阶龙格-库塔公式得 zy)22(6)22(64321143211LLLLhzzKKKKhyyi

40、iii),()2,2,(2)2,2,(2),(3314342221323112121211hLzhKyxfLhLzKLhzKhyxfLLhzKLhzKhyxfLLhzKzyxfLzKiiiiiiiiiiiiiiii例例 求解下列二阶微分方程的初值问题求解下列二阶微分方程的初值问题 1)0(, 0)0(yyxyy10 x取步长取步长h=0.1 h=0.1 解解: :先作变换：令先作变换：令 ，代入上式，得一阶方程组，代入上式，得一阶方程组 yz1)0(, 0)0(,zyzyxzz用四阶龙格用四阶龙格- -库塔方法求解库塔方法求解, ,按式按式(9.37)(9.37)及及(9.38)(9.38)进

41、行计算：进行计算：取步长取步长 , , , , ,1 . 0h00 x00y10z0i时时 2105. 1)105. 11 . 01 () 1 . 00()()(1105. 1105. 11 . 01105. 1) 1 . 121 . 01 ()21 . 00()2()2(055. 11 . 121 . 0121 . 1) 121 . 01 ()21 . 00()2()2(05. 1121 . 012101),(130043042003203100210200000101hLzhxLhLzKLhzhxLLhzKLhzhxLLhzKxzzyxfLzK1104. 1)2105. 1105. 121 . 121 (61 . 01)22(61053. 0)1105. 1055. 1205. 121 (61 . 00)22(6432101432101LLLLhzzKKKKhyy然后计算然后计算 时的时的y y2 2和和z z2 2；依此类推，直到；依此类推，直到i=9i=9时的时的y y1010和和z z1010，即可得到，即可得到其数值解。其数值解。 1i44332211,LKLKLKLK6.7 边值问题的数值解法边值问题的数值解法 打靶法与有限差分法ybyaybaxyyxfy求解 )(,)(, ),(