龙格库塔算法

常微分方程的数值解法本章主要内容：1.欧拉法2.欧拉公式的截断误差3.改进欧拉法4.龙格—库塔法重点：欧拉法、龙格—库塔法难点：局部截断误差，龙格—库塔法14.1欧拉法14.1.1欧拉公式本章的学习目的是要通过数值解法，求解一阶微分方程的初值问题。设一阶微分方程的初值问题为怎样求呢？我们过点以为斜率作切线切线方程为以代入，得再以x1作为x0，用上面的方法作切线求y2，如此下去得到，从而求得的近似值。这种解法称为欧拉法。00)(),('yxyyxfdxdyy)(xy),(00yx))(,(0000xxyxfyy),()('000yxfxy1xx))(,()(0100011xxyxfyyxy))(,()(111nnnnnnnxxyxfyyxy)(),(),(121nxyxyxy如果的取法是等步长的。记上面用欧拉法求解的公式可表示成：这一公式称为欧拉公式。欧拉法的几何意义：nnxxxxxxh11201X1210,,nxxxxhkxhxxyxfhyyxyknnnnnn)1(),()(0111Y)(xyy0x1x2x3x4x5x例1用欧拉法求初值问题在处的近似值。解：列表计算，)2(1.0),(1kkkkkkkkyxyyyxfhyy0.1,,2.0,1.0x1.0,1,0,2),('00hyxyxyyxfy)9,,2,1,0(k1)0(2yyxydxdy264911.1)3.0(,183216.1)2.0(,095445.1)1.0(yyy得784771.1)0.1(,y局部截断误差精确解与近似解之间的误差称为局部截断误差。由泰勒展开式可知：因此，欧拉法具有一阶精度。它的局部截断误差是关于步长h的二阶无穷小量。即：21!2)()(')(hxyhxyyxykkkk)()(211hOyxykk)(1kxy1ky2!2)(),(hxyhyxfykkkk)(21hOyk14.1.2改进欧拉法求微分方程的解，最简单的想法就是直接积分这在高等数学中可能是不可行的，但在数值分析中却是可行的。用数值积分的梯形公式代入，有：这是求的近似值的梯形公式。由于等式的两边都有，不好直接计算，称为隐式形式。得替代用),(),(,11kkkkxyxyyydxxyxfxyxykkxxkk))(,()()(11))](,())(,([2)()(111kkkkkkkxyxfxyxfhxyxy)],(),([2)(1111kkkkkkkkyxfyxfhyyxy)(1nxy1ny但由于是近似计算，我们可以先用欧拉公式求出的一个近似值，然后把这个近似值代入等号的右边，计算等号左边的，前者称为预报值，后者称为校正值。这种方法称为改进欧拉法。公式可表示为：公式也可表示为：为常数取等距节点上式中的hxxhkkk,,1)],(),([2),(1111kkkkkkkkkkkkyxfyxfhyyyxfhyy校正值预报值1ny1ny)],(),([2),(1111kkkkkkkkkkyxfyxfhyyyxhfyy校正值预报值写成一个式子。就是：也可表示为平均值的形式：改进欧拉法具有二阶精度。它的局部截断误差是关于步长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(,y[2001年7月试卷计算题14]用改进的欧拉法预报—校正公式，取步长h=0.2求解初值得问题解：本题改进欧拉法的预报—校正公式为：当k=0时，1)0()4.00(2'yxxyy)22(1.04.01111kkkkkkkkkkyxyxyyyxyy校正值预报值2.0,1,0100xyx96.02.021.01)22(1.014.01100010001yxyxyyyxyy当k=1时，4.0,96.0,2.0211xyx8509.0)8832.04.096.02.0(2.096.0)22(1.08832.096.02.04.096.04.02211121112yxyxyyyxyy8509.0)4.0(y14.2龙格—库塔法14.2.1龙格—库塔法的基本思想设一阶微分方程的初值问题的解为。如果是等距节点，记步长根据微分中值定理上式中而根据题设故有：00)(),('yxyyxfdxdyy)(xyyhyyxxxyxyhxyykkkkkkk111)()()(')('kkxxh1))(,()('),,('hxyhxfhxyyxfykkk10),,(1则kkxx))(,(1hxyhxhfyykkkk由于是中满足微分中值定理的点的导数因此，称为平均斜率则：如取处的斜率作为平均斜率的近似值，则得欧拉公式，它的近似等级为。如取处的斜率的平均数作为平均斜率的近似值，则得改进欧拉公式，它的近似等级为。公式为：)('y),(1kkyxf))(,(hxyhxfkk),(1kkxxhyykk1))(,(hxyhxfkkkx),(kkyxf)(2hO1,kkxx21,)(3hO)),(,(),(1112kkkkkkyxhfyxfyxf2/)(211hyykk由此想到，我们是否可以在中多取几个点，再以它们的某种平均值作为，把精度进一步提高呢？假设在中取n个点其中，斜率为如果取第一个点为，则取平均斜率的计算方法为则：其中：恰当的选取公式中的常数，就可以使精度尽量地高。这就是龙格—库塔法的一般公式。],[1kkxx1,,,021n),(11mjjmjkmkmhyhxfnmmmkkhyy11],[1kkxxhxhxhxnkkk,,,21n,,,21kxnn2211),(,011kkyxf14.2.2二阶龙格—库塔法假设在中再取一个点如果取第一个点为，则取点为第二个点，则：公式为：恰当的选取，可以使精度达到。这就是二阶龙格—库塔公式。),(1kkyxf),(),()(12122111phyphxfyxfhyykkkkkk],[1kkxxphxxkpkkxn2211pkx),(),(12phyphxfyxfkkpkpkp,,21)(3hO常见的二阶龙格—库塔公式称为中点公式。它是二阶龙格—库塔公式中的结果。21,1,021p)2,2(),(12121hyhxfyxfhyykkkkkk14.2.3三、四阶龙格—库塔法三阶龙格—库塔法的一般公式为其中：该公式的局部截断误差为。)(4hO))2(,()2,2(),(213121hyhxfhyhxfyxfkkkkkk3211616461kkyy四阶龙格—库塔法的一般公式为其中：该公式的局部截断误差可达。)(5hO),()2,2()2,2(),(3423121hyhxfhyhxfhyhxfyxfkkkkkkkk)22(643211hyykk四阶龙格—库塔法的优点是；⑴它是一步法，即已知yk就可以求出yk+1，不需要知道其它的数据。⑵精确度高。缺点是计算量大。在一个步长的计算中，要四次计算f(x,y)的值。例如，用四阶龙格—库塔法再解在处的近似值。解：1)0(2yyxydxdy0.1,,2.0,1.0x)22(643211hyykk4,3,2,1,0,2.0,1,000khyx我们仍然用Excel列表计算，为了清楚，先用公式计算出，再计算出，然后按步长增加，计算下一步长。计算结果为：341667.1)4.0(,183229.1)2.0(yy4321,,,1ky),()2,2()2,2(),(3423121hyhxfhyhxfhyhxfyxfkkkkkkkk612514.1)8.0(,483281.1)6.0(yy732142.1)0.1(y解常微分方程初值问题的三阶龙格——库塔法的局部截断误差是。作业：P.185P.203带※的练习题)(4hO