计算机数学基础(下)

1、计算机数学基础(下)第5编 数值分析第第1414章章 常微分方程的数值解法常微分方程的数值解法本章主要内容本章主要内容：1. 欧拉法欧拉法 2. 欧拉欧拉公式的截断误差公式的截断误差 3. 改进欧拉改进欧拉法法4. 龙格龙格库塔法库塔法 重点：重点：欧拉法欧拉法、龙格、龙格库塔法库塔法 难点：局部截断误差难点：局部截断误差 ，龙格，龙格库塔法库塔法14.1 欧拉法 14.1.1 欧拉公式 本章的学习目的是要通过数值解法，求解一阶微分方程的初值问题。 设一阶微分方程的初值问题为 怎样求 呢？ 我们过点 以 为斜率作切线切线方程为以 代入，得再以x1作为x0，用上面的方法作切线求y2，如此下去得到

2、 ，从而求得 的近似值。这种解法称为欧拉法。00)(),(yxyyxfdxdyy)(xy),(00yx)(,(0000 xxyxfyy),()( 000yxfxy1xx )(,()(0100011xxyxfyyxy)(,()(111nnnnnnnxxyxfyyxy)(),(),(121nxyxyxy如果 的取法是等步长的。记 上面用欧拉法求解的公式可表示成： 这一公式称为欧拉公式。欧拉法的几何意义： nnxxxxxxh11201X1210,nxxxxhkxhxxyxfhyyxyknnnnnn) 1(),()(0111Y)(xyy 0 x1x2x3x4x5x例1 用欧拉法求初值问题在 处的近似值

3、。解： 列表计算，)2(1 . 0),(1kkkkkkkkyxyyyxfhyy0 . 1 , 2 . 0 , 1 . 0 x1 . 0, 1, 0,2),(00hyxyxyyxfy)9 , 2 , 1 , 0(k1)0(2yyxydxdy264911. 1)3 . 0(,183216. 1)2 . 0(,095445. 1) 1 . 0(yyy得784771. 1)0 . 1 (,y局部截断误差 精确解 与近似解 之间的误差称为局部截断误差。 由泰勒展开式可知： 因此，欧拉法具有一阶精度。 它的局部截断误差是关于步长 h 的二阶无穷小量。 即： 21! 2)()( )(hxyhxyyxykkk

4、k)()(211hOyxykk)(1kxy1ky2! 2)(),(hxyhyxfykkkk)(21hOyk2002年1月试卷选择题5 解常微分方程初值问题的欧拉法的局部截断误差是( )。)()()45hOBhOA)()()23hODhOCD14.1.2 改进欧拉法 求微分方程的解，最简单的想法就是直接积分这在高等数学中可能是不可行的，但在数值分析中却是可行的。用数值积分的梯形公式代入，有：这是求 的近似值的梯形公式。由于等式的两边都有 ，不好直接计算，称为隐式形式。得替代用),(),(,11kkkkxyxyyydxxyxfxyxykkxxkk)(,()()(11)(,()(,(2)()(111

5、kkkkkkkxyxfxyxfhxyxy),(),(2)(1111kkkkkkkkyxfyxfhyyxy)(1nxy1ny但由于是近似计算，我们可以先用欧拉公式求出的一个近似值，然后把这个近似值代入等号的右边，计算等号左边的 ，前者称为预报值，后者称为校正值。这种方法称为改进欧拉法。公式可表示为：公式也可表示为：为常数取等距节点上式中的hxxhkkk,1),(),(2),(1111kkkkkkkkkkkkyxfyxfhyyyxfhyy校正值预报值1ny1ny),(),(2),(1111kkkkkkkkkkyxfyxfhyyyxhfyy校正值预报值写成一个式子。就是：也可表示为平均值的形式：改进

6、欧拉法具有二阶精度。它的局部截断误差是关于步长h 的三阶无穷小量。即： 2/ )(),(),(11cpkpkkckkkpyyyyxhfyyyxhfyy),(,(),(211kkkkkkkkyxhfyxfyxfhyy)()(311hOyxykk例2 用改进欧拉法求例1初值问题的近似解。解：本题用改进欧拉法求解的公式为： 列表计算，)(21)2(1 . 0)2(1 . 011cpkpkpkckkkkpyyyyxyyyyxyyy266201. 1)3 . 0(,184097. 1)2 . 0(,095909. 1) 1 . 0(yyy得737867. 1)0 . 1 (,y2001年7月试卷计算题1

7、4 用改进的欧拉法预报校正公式，取步长h=0.2求解初值得问题解：本题改进欧拉法的预报校正公式为：当 k=0 时，1) 0() 4 . 00(2yxxyy)22( 1 . 04 . 01111kkkkkkkkkkyxyxyyyxyy校正值预报值2 . 0, 1, 0100 xyx96. 02 . 021 . 01)22( 1 . 014 . 01100010001yxyxyyyxyy当 k=1 时，4 . 0,96. 0, 2 . 0211xyx8509. 0)8832. 04 . 096. 02 . 0(2 . 096. 0)22( 1 . 08832. 096. 02 . 04 . 096

8、. 04 . 02211121112yxyxyyyxyy8509. 0)4 . 0(y14.2 龙格库塔法 14.2.1 龙格库塔法的基本思想 设一阶微分方程的初值问题的解为 。 如果是等距节点，记步长 根据微分中值定理上式中 而根据题设故有：00)(),(yxyyxfdxdyy)(xyy hyyxxxyxyhxyykkkkkkk111)()()( )( kkxxh1)(,()( ),(hxyhxfhxyyxfykkk10),(1则kkxx)(,(1hxyhxhfyykkkk由于 是 中满足微分中值定理的点的导数因此， 称为平均斜率则：如取 处的斜率 作为平均斜率 的近似值，则得欧拉公式，它的

9、近似等级为 。如取 处的斜率 的平均数作为平均斜率的近似值，则得改进欧拉公式，它的近似等级为 。公式为：)( y),(1kkyxf)(,(hxyhxfkk),(1kkxxhyykk1)(,(hxyhxfkkkx),(kkyxf)(2hO1,kkxx21,)(3hO),(,(),(1112kkkkkkyxhfyxfyxf2/ )(211hyykk由此想到，我们是否可以在 中多取几个点，再以它们的某种平均值作为 ，把精度进一步提高呢？假设在 中取n个点其中 ，斜率为如果取第一个点为 ，则 取平均斜率的计算方法为则：其中：恰当的选取公式中的常数，就可以使精度尽量地高。这就是龙格库塔法的一般公式。 ,

10、1kkxx1,021n),(11mjjmjkmkmhyhxfnmmmkkhyy11,1kkxxhxhxhxnkkk,21n,21kxnn2211),(, 011kkyxf14.2.2 二阶龙格库塔法 假设在 中再取一个点如果取第一个点为 ，则 取点 为第二个点，则：公式为：恰当的选取 ，可以使精度达到 。这就是二阶龙格库塔公式。 ),(1kkyxf),(),()(12122111phyphxfyxfhyykkkkkk,1kkxxphxxkpkkxn2211pkx),(),(12phyphxfyxfkkpkpkp,21)(3hO常见的二阶龙格库塔公式称为中点公式。它是二阶龙格库塔公式中的结果。2

11、1, 1, 021p)2,2(),(12121hyhxfyxfhyykkkkkk14.2.3 三、四阶龙格库塔法 三阶龙格库塔法的一般公式为其中：该公式的局部截断误差为 。)(4hO)2(,()2,2(),(213121hyhxfhyhxfyxfkkkkkk3211616461kkyy四阶龙格库塔法的一般公式为其中：该公式的局部截断误差可达 。)(5hO),()2,2()2,2(),(3423121hyhxfhyhxfhyhxfyxfkkkkkkkk)22(643211hyykk四阶龙格库塔法的优点是； 它是一步法，即已知yk就可以求出yk+1，不需要知道其它的数据。 精确度高。 缺点是计算量大。在一个步长的计算中，要四次计算f(x,y)的值。例如，用四阶龙格库塔法再解在 处的近似值。解：1)0(2yyxydxdy0 . 1 , 2 . 0 , 1 . 0 x)22(643211hyykk4 , 3 , 2 , 1 , 0, 2 . 0, 1, 000khyx我们仍然用 Excel 列表计算，为了清楚，先用公式计算出 ，再计算出 ，然后按步长增加，计算下一步长。计算结果为：341667. 1)4 . 0(,183229. 1)2