RK求解微分方程

Runge-Kutta积分方法所以得到：是精确的，中的平均速度。设是动点在其中为：，一般的解法可以表示对)(!3)(2)()()()(),(),().,(),(32111nnnnnnnnnnnnnnntYhtYhtYhtYhtYtYYttYtDYtDhYYYtFdtdY)(!)(!3)(2)()()(),()(121nrrnnnnnnntYrhtYhtYhtYhtYtYYtD由此得到高阶的单步法。但是，往往右函数的高阶导数或者无法直接得到、或者计算太过复杂。所以实际的做法是：用[tn,tn+1]区间中解曲线邻域的一些已知点函数值的线性组合来代替F(t,Y)的导数，从而得到高阶的单步法公式。),(),(2121)()),(,(),(2)()()),(,(),(2EulerEuler1212113111KYhthFKYthFKKKYYhOYtFhYhtFYtFhtYtYdYtFhYhtFYtFhYYnnnnnnnnnnnnnnnnnnnnnnn：上还可以写成如下形式另一方面，该公式实际所以它是二阶的。它的局部截断误差为校正方法：型的预估法，即改进的例所以它是一阶的。：，它的局部截断误差为法：).(),()()(),(Euler2111hOYthFtYtYdYtFhYYnnnnnnnnn此处即通过计算已知点的函数值(K1,K2)的线性组合代替高阶导数，得到了较高的精度。Runge-Kutta方法的推导Runge-Kutta方法的一般形式：1111111,,3,2,),(),(ijijiijjijnininnriiinnriKhYhtFKYtFKKchYY确定了阶数之后，再通过Taylor展开、比较两边系数的方法，确定各待定系数：,,iiijc二阶显式Runge-Kutta方法.),(,),()()(),(),()(212211112122122111nnnnnnnnnnnnnnnYthFYhtFcYtFchtYtYdhKYhtFKYtFKKcKchYY：此处的局部截断误差为)()(!3)(!2)()()()(4321hOtYhtYhtYhtYhtYtYnnnnnn)(),(),(),(2212212hOYtFhFYtFhFhFYhtFnnYnnntnnnn展开各项如下：nnnYnntnnYnnnYYnnntYnnttnnnnYnntnnnYnntnnnnnnn)F,Y(tF),Y(tF),Y(tF)F,Y(tF)F,Y(tF),Y(tFY)F,Y(tF),Y(tFdt)dY(t),Y(tF),Y(tFdt)),Y(tdF(tYF),YF(tY22其中二阶显式Runge-Kutta方法).(),(21),(21)1(32212222211hOFYtFhcYtFhchFccdnnnYnntnn将展开式代入,得到局部截断误差：要使得方法是二阶的，则局部截断误差应该为三阶小量，即：021;021;012122221cccc三个方程,四个未知数,所以其解不唯一.2122110,1;.2cacaa可令则有：法。即为改进的积分公式：时Euler),(),(2)(.1,21,2112121121221hKYhtFKYtFKKKhYYccannnnnn即为中点公式。积分公式：时)2,2(),(.21,1,0,11212121221hKYhtFKYtFKhKYYccannnnnn例.1)0(,:ODE2yydxdy求解初值问题.11xy易知其精确解为：21.01.0.1,21Euler,21.1222121221nnnnnyyyyycca积分公式：法：改进的321.0321.0.23,31,32,31.2222121221nnnnnyyyyycca积分公式：分别用以下两种系数：步长都取为1.0h结果及比较三阶显式Runge-Kutta方法nnnYnntnnYnnnYYnnntYnnttn)F,Y(tF),Y(tF),Y(tF)F,Y(tF)F,Y(tF),Y(tFY22在推导二阶显式方法的过程中，注意到局部截断误差表达式中h3项包含了以下表达式：因此若要在局部截断误差中消去h3项，必须增加包含了以上各项的多个方程，同时我们注意到r=2时，只有等四个待定系数，少于方程的数目，所以这样的系数不存在。故：r=2时Runge-Kutta方法只能是二阶的。要得到三阶的方法，则必须有r=3。),(),(),()(232131331212213322111hKhKYhtFKhKYhtFKYtFKKcKcKchYYnnnnnnnn其局部截断误差为：).()()(33221111KcKcKchtYtYdnnn1,2121,,c三阶显式Runge-Kutta方法程如下：可得待定系数满足的方展开，使得作以及将),(Taylor)(,41132hOdtYKKnn61312113223233222332232313212321cccccccc)2,()2,2(),()4(62131213211hKhKYhtFKKhYhtFKYtFKKKKhYYnnnnnnnn三阶公式，如下：常见的公式称为公式。特别地，一个三阶，它们统称为因此可以得到众多公式个未知数，解不唯一。个方程要决定KuttaKuttaRunge86四阶显式Runge-Kutta方法个：。下面列出最常见的一的局部截断误差满足公式，它们导出各种四阶的类似前面的推导，可以)(KuttaRunge51hOdn),()2,21()2,21(),(226342312143211hKYhtFKKhYhtFKKhYhtFKYtFKKKKKhYYnnnnnnnnnn四阶显式Runge-Kutta方法xnxn+h/2xn+hf1f2f3f443212261ffffff.1)0(,:ODE2yxydxdy求解初值问题xexxy222易知其精确解为：方法求解：分别用二阶、四阶步长都取为KR1.0hx四阶二阶真解四阶误差二阶误差0.01.0000001.0000001.0000000.00000.0000000.11.1048291.1024501.1048291.60E-72.38E-30.21.2185971.2115071.2185973.40E-77.09E-30.31.3401411.3257661.3401415.48E-71.44E-20.41.4681751.4436711.4681757.69E-72.45E-20.51.6012781.5635061.6012799.95E-73.78E-20.61.7378801.6833741.7378811.20E-65.45E-20.71.8762461.8011791.8762471.42E-67.51E-20.82.0144571.9146032.0144591.68E-69.99E-20.92.1503952.0210862.1503971.96E-61.29E-11.02.2817162.1178002.2817182.32E-61.64E-1例结果及比较结果及比较.1045.1505.0).10(1060.1101040.2510.01514时误差为而二阶公式，相对误差仅为仍然是相当精确的结果时误差为时误差为对四阶公式，xhxxh关于Runge-Kutta方法RungeKuttaRungeKutta类似前面的推导,可以导出更高阶的公式.关于方法,有以下几点需要特别指出：。解曲线比较光滑的情形别适用于展开的方法，因此它特方法的推导基于TaylorKuttaRunge.1次右函数。、、分别须计算阶数相同，即它们每步数的次数和方法，每一步计算右函二阶、三阶、四阶的432Kutta-Runge.27)9(,6)8(,6)7(,5)6(,4)5()(KuttaRunge)4(.3NNNNNvNvN阶数，则有：次右函数可达到的最高表示计算若用比阶数大。次的次数方法每步须计算右函数阶的的波动。，局部误差会有比较大如果采用固定步长计算的步长等等因素相关。微分方程的性质、采用方法具体的系数、待解阶数、比较复杂，它和方法的法的局部截断误差估计KuttaRunge.4提高Runge-Kutta方法的精度的方法提高积分方法的精度,我们最熟悉的(不一定是最好的)措施是1()212()2212Euler()()211()()24nnnhhyyhyyxyxchchhyxyxchch我们用一个例子予以说明如下法的近似解：将步长减半为时，有外推法Richardson0h将上二式作适当线性组合,可使的一次项为：)()(2)()()2(xyxyxyhh。计算量同时增加了一倍但我们注意到右函数的得精度提高一阶。利用此式计算时，可使提高精度最简单的方法是缩短步长，但要以牺牲计算速度和积累舍入误差为代价。变步长的Runge-Kutta方法作为妥协，如果能在计算过程中实时控制步长的大小，就可以在获得较高的计算速度的同时，保证较高的精度。()()()111111()()()()2222111pRungeKutta()2()222hhhpnnnphhhhpnnnpdYtYchhhhdYtYcc一般地,设有阶的公式,其局部截断误差为将步长减半为时,有()()2()()1211()()()()()()2221111111(1)211;.1221hhhhpnnphhhhhhnnnnnnppcccYYchdYYdYY假定,记为,则可以估计误差如下：因此可以从两次计算当中估计出每一步的截断误差,有了这个误差估计之后,通过与控制误差限比较,就可以控制步注意这个方法增加了长.计算量.Runge-Kutta-Fehlberg方法Fehlberg设计了一个更加精巧的嵌套方法如下：pp+1:在采用一个阶方法的同时,计算一个阶的结果,并由此给出误差估计11p()pnnYYthch阶的方法：#21p1()pnnYYthch阶的方法：p因此阶方法的局部截断误差可以近似为：.1#11nnnYYd制步长。可以用这个估计式来控RungeKutta,,pp+1iiijc我们记得方法中的待定系数不是唯一确定的.因此就有可能利用这个特点,选择适当的系数使得阶和阶公式中尽可能多的系数相同,从而达到减少计算量的目的.pp+1ic实际上,可以找到仅有权系数不同的阶和阶算法,如此一来,右函数的计算就可以大大减少.Runge-Kutta-Fehlberg方法Fehlberg给出的四阶、五阶公式RKF4(5)如下：552040114104185925653544227821509514104845513368082164391