8.2runge-kutta方法

《8.2runge-kutta方法》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、第八章常微分方程数值解法8.2Runge-Kutta方法8.2.1Runge-Kutta方法的基本思想8.2.2几类显式Runge-Kutta方法第八章常微分方程数值解法8.2.1Runge-Kutta方法的基本思想显式Euler方法是最简单的单步法，它是一阶的，它可以看作Talylor展开后取前两项。因此，得到高阶方法的一个直接想法是用Talylor展开，如果能计算y(x)的高阶导数，则可写出p阶方法的计算方法2phh()py=y+hy′+y′′+L+y，n+1nnnn2p!(j)()j()其中yn是yxn的近似值，j=0，1，2，Lp。若将f(x，y)，∂f∂x，∂f∂y，L，分别记成f，

2、f，f，L，则对于二阶和三阶导数可表示为xyy′′=f+ff，xy22y′′′=f+2ff+ff+ff+ff。xxxyxyyyy第八章常微分方程数值解法这个方法并不实用，因为一般情况下，求f(x，y)的导数相当麻烦。从计算高阶导数的公式知道，方法的截断误差提高一阶，需要增加的计算量很大。但是由此启发我们用区间上若干个点的导数f，而不是高阶导数，将它们作线性组合得到平均斜率，将其与解的Taylor展开相比较，使前面若干项吻合，从而得到具有一定阶的方法。这就是Runge-Kutta方法的基本思想，其一般形式为Ly=y+h∑λK，n+1niii=1(8.2.1)K=f()x，y，1nn⎛i−1⎞Ki

3、=f⎜⎜xn+cih,yn+cih∑aijKj⎟⎟，i=2，3，L，L，⎝j=1⎠第八章常微分方程数值解法Li−1其中，ci≤1，∑λi=1，∑aij=1。它的局部截断误差是i=1j=1L()()*Tn+1=yxn+1−yxn−h∑λiKi，（8.2.2）i=1*Ky其中，与Kii的区别在于：用微分方程准确解y()xn代替K中的n就得i*，ca到K。参数λii和ij待定，确定它们的原则和方法是：将（8.2.2）式中iy(x)*的在n+1xn处作Taylor展开，将Ki在(xn，y()xn)处作二元Taylor展T开，将展开式按H的幂次整理后，令n+1中h的低次幂的系数为零，使Tn+1首(p+1

4、)项中h的幂次尽量高，比如使T=Oh，则称（8.2.1）式为L级p阶n+1Runge-Kutta方法（简称R-K法）。第八章常微分方程数值解法类似于显式R-K公式（8.2.1），稍加改变，就得到隐式R-K公式Lyn+1=yn+h∑λiKi，i=1⎛L⎞Ki=f⎜⎜xn+cih，yn+cih∑aijKj⎟⎟，i=1，2，L，L。⎝j=1⎠它与显式R-K公式的区别在于：显式公式中，对系数a求和的上限是i−1，从ij而aij构成的矩阵是一个严格下三角阵。而在隐式公式中，对系数aij求和的上限是L，从而aij构成的矩阵是方阵，需要用迭代法求出近似斜率Ki=(i=1，2，L，L)推导隐式公式的思路和方法

5、与显式公式法类似。第八章常微分方程数值解法8.2.2几类显式Runge-Kutta方法对于L=2，则y=y+h(λK+λK)，n+1n1122K=f()x，y，1nnK=f()x+ch，y+chK。2n2n21其局部截断误差是**T=y(x)−y(x)−h(λK+λK)（8.2.3）n+1n+1n1122第八章常微分方程数值解法将Tn+1中的各项作Taylor展开，并利用y′(xn)=f(xn，y(xn))，y′′=fx+fyf，则有23()()()h()h()()4yx=yx+hy′x+y′′x+y′′′x+Oh，n+1nnnn26*K=f()x，y()x=y′(x)，1nnn*K=f()x

6、+ch，y()x+chy′()x，2n2n2n22()''()c2h()2(3)=y′x+chyx+f+2ff+ff+Oh，n2nxxxyyy2将它们代入（8.2.3）式，整理后得()()⎛1⎞2()Tn+1=1−λ1−λ2hy′xn+⎜−λ2c2⎟hy′′xn⎝2⎠23⎡1()λ2c2()2⎤(4)+hy′′′x−f+2ff+ff⎥+Oh⎢nxxxyyy⎣62⎦第八章常微分方程数值解法c23选取λ1，λ2和2，使方法的阶尽可能高，就是使h和h的系数为零，因为h的系数一般不为零。于是得到方程组λ+λ=1，121λc=。222显然，该方程组有无穷多组解，从而得到一族二级二阶R-K方法。若以c为自

7、由参数，取c=12得中点公式22⎛hh()⎞yn+1=yn+hf⎜xn+，yn+fxn，yn⎟，（8.2.4）⎝22⎠取c=2/3得Heun公式h⎡⎛22⎞⎤yn+1=yn+⎢f()xn，yn+3f⎜xn+h，yn+hf()xn，yn⎟⎥，（8.2.5）4⎣⎝33⎠⎦第八章常微分方程数值解法取c=1得改进的Euler公式（8.1.6）。对于L=3的情形，要计算三个斜率的近似值：K=f(x，y)，1