三阶Runge-Kutta方法

2012-2013（1）专业课程实践论文三阶Runge-Kutta方法王曹旭，0818180106，R数学08-1班常微分方程初值问题的数值解法是求方程（1）的解在点列1(0,1,)n n n x x h n -=+= 上的近似值n y ，这里n h 是1n x -到n x 的步长，一般略去下标记为h 。

00(,)()dyf x y dxy x y⎧=⎪⎨⎪=⎩（1）三阶Runge-Kutta 方法，它的计算公式是：⎪⎪⎩⎪⎪⎨⎧+++=++==+++=+))((,(),,(),,(),(2131213322111sK rK qh y qh x f K phK y ph x f K y x f K K K K h y y n n n n n n n n λλλ R K-方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算三次函数值f 。

在用龙格库塔方法时,要注意N 的选择要合适，N 太大，会使计算量加大，N 太小，h 较大，可能会使误差增大。

因此选择合适的N 很重要。

我们要在考虑精度的基础上，选择合适的N 。

在此，用c 语言实现了三阶Runge-Kutta 方法。

三阶Runge-Kutta流程图#include<stdlib.h>#include<stdio.h>/*n表示几等分，n+1表示他输出的个数*/int RungeKutta(double y0,double a,double b,int n,double *x,double *y,double (*function)(double,double)){double h=(b-a)/n,k1,k2,k3;int i;// x=(double*)malloc((n+1)*sizeof(double));// y=(double*)malloc((n+1)*sizeof(double));x[0]=a;y[0]=y0;for(i=0;i<n;i++){x[i+1]=x[i]+h;k1=function(x[i],y[i]);k2=function(x[i]+h/2,y[i]+h*k1/2);k3=function(x[i]+h,y[i]-h*k1+2*h*k2);y[i+1]=y[i]+h*(k1+4*k2+k3)/6;}return 1;}double function(double x,double y){return y-2*x/y;}//例子求y'=y-2*x/y(0<x<1);y0=1;int main(){ int i;double x[7],y[7];printf("用三阶龙格-库塔方法\n");RungeKutta(1,0,1,6,x,y,function);for(i=0;i<7;i++)printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]);return 1;}例1．求解方程⎩⎨⎧=<<+='1,10,0y x y x y 步长2.0=h解：（1）将程序序中rerurn y-2*x/y 改为rerurn x+y（2）输入,,,,0n b a y 的值：;5;1;0;1（3）显示输出结果：435019.3y 650070.2043616.2583310.1242667.1000000.1000000.1800000.0600000.0400000.0200000.0000000.05432054320============y y y y y x x x x x x 时，时，时，时，时，时，例2．求解方程⎩⎨⎧=<<+='110),1/(30y x x y y 步长2.0=h解：（1）将程序序中rerurn y-2*x/y 改为rerurn 3*y/(1+x)(2) 输入n b a y ,,,0的值：;5;1;0;1(3) 显示输出结果：961755.7805800.5079391.4734787.2724242.1000000.1000000.1800000.0600000.0x 400000.0x 200000.0000000.0543210543210============y y y y y y x x x x 时，时，时，时，时，时，。

龙格-库塔方法〔Runge-Kutta〕3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法，若用显式单步法(3.2.2)当，即数值求积用左矩形公式，它就是Euler法(3.1.2)，方法只有一阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶精度的.若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)一样，用前面已算得的f值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值(即由前面(i-1)个已算出的表示，故公式是显式的.例)，而ki如当r=2时，公式可表示为(3.2.5)其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值，可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5)，要求选择参数，使公式的精度阶p 尽量高，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6)令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯一.因r=2，由〔3.2.5〕式可知，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只要都可得到二阶精度R-K方法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)与中点公式(3.2.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即其中将按二元函数在处按Taylor公式展开，然后代入局部截断误差表达式，可得可得三阶方法，其系数共有8个，所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组，解也是不唯一的，通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法：(3.2.11)附：R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法与步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。

Runge-Kutta方法1.(,,)(:,(,)3),.m m m m u u h t u h t u h u ϕϕ+-=构造单步高精度方法单步方法的一般形式如下为了提高截断误差的阶可以是的非线性函数()()001(12,)()(,.)()m m m m m m u f t u u t u u u hf Eule t u r h u hf o +'=⎧⎨=⎩=+=+方法：截断误差的阶常微分方程()()()1111(())111()1(,,)!!(!(4),())()q j j q q j j j j j q j j j d f t u t u Tayl h t u h h h j j or u t h u t h u t o h j t dt ϕ--===+-===+∑∑∑展开：Runge Kutta -方法1(3),(),(5)1().,,m m m m u u h t q u t u h q ϕ+-=单步方法称为阶的如果对于微分方程的精确解是使成立的最定大整数义()00()()(,,)(,,)(,),()()(,,5).()()q u t h u t h t u t u h u f t u u t u u t h t u h u t h t Tayl h o h or q ϕϕϕ+-=+'==++当满足(4)形式时,满足条件:与在点的级数次数不超过阶的项重合0010,0,()()lim lim (,,),(,)lim (,,)2(3)(,)(,,)(,,0)=(,(,,)).m m m h m h h u h u t h u t t u h u f t u t u h f t u ht u h ht u f t u u u h t u h ϕϕϕϕϕ→→+→→+-'==⇒==- 若是微分方程的解则当时此时称单步法为相容的定义.当关于连续时,相容性等价于()111111(6)(,,),(,)(7),,2,3,,(8),2,3,,Nn n n n l l n n nl n n n l l t u h k k f t u k f t h u N Runge Kutta h k N n Nc b n b ϕαα=-=-===⎧⎪⎛⎫⎨-=++= ⎪⎪⎝⎭⎩==∑∑∑级方法:这里()(,,)(6)-(8),.3Runge Kutta t u h N ϕ-由定义时称算法为级当方法1()-()(,,)(,,(,,)(3))m m m m u t h u u u h t h t u h t t u h u h ϕϕϕ++=-=问题关键的构造()()()()3()31''22231()(()(),)!()()(),,(.))(9)(j j j u t u u t t t u uu t tu u t u t f t u u t h u t u t f t u Ru u nge K t utta f ff f ff h u t f f f ff f f t o f h j u ==+=+==+=+++++-+∑三级方法推导首先注意满足微分方程：1222133332132231(,)(,)(,()(,,)(1))1n nn k f t u f k f t h u h t u k k f t h u h b k hb k h k c ααααϕ=⎧⎪⎪⎪==⎨⎪=++⎪=++-+⎩=⎪∑()233211()()()()1026,2u t u tt tu uuu t h u t hf h F h Ff G o h F f ff G f ff f f ++++++==++=+多元函数的Taylor 公式()()()()0000000(1),,Taylor ,,,mm m i i m i m i m i i f x y h k f x y C h kx y f x y x y n x y --=∂⎛⎫∂∂+= ⎪∂∂∂∂⎝⎭∑式称为在点的阶公式用二项展开形式表述为()()()()()()()()()()00210000000000,1,=1,1!0,1Lag ,,,2range !1,!1n n n n n n P h k f x y h k f x y h k f x y x y f x h y k R R R h k f x h y k n x y h k f x y n x y x y θθθ+⎛⎫⎛⎫∂∂∂∂++++ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎛⎫∂∂+++ ⎪∂∂⎝⎭+++=⎛⎫∂∂=+++ ⎪+∂∂⎝+⎭∈其中称为余项.()()()()()0000000,,;1,;+(Tayl ) ,or ,f x y P x y U P n U P x h y k δδ=++设函数在点的某个邻域内有直到阶的连续偏导数定理1公式则对中的任意一点成立：()122213333213222222221211222222(,)(,)(11)(,)(,())112()(2)()21()22t u tt tu uu k f t u f k f t h u h k k f t h u h b k hb k k f h f k f h f k f k t u Ta k f o h f h F h G y o h r o l αααααααα==⎧⎪=++⎨⎪=++-+⎩=++++++=+++将在点展成次多项式:()333321322223333213222222233232332132213{[()]}1{2[()]2[(1()()2)]}()t u tt tu u u u k f h f b k b k f h f b k b k f b k b k f o f h F h b F h h f G o ααααααααα=++-+++-++-++++=++类似得到:()()()()123223322223232223133(,,)()()141[2()]()26(,,)1213,(,,)n n n u t u h c c c f h c c Fh c b F t u h c f c c u G t o h k h ϕααϕαααϕ==++++++++-=∑将代入式得到的展开式()()1,2,3149N Runge Kutta =-比较与即可导出时的方法()()223110:(,,114)15)(1c c t u h c f o h c N Euler ϕ=+====当，式化为令,是时:方法23311(9)()()()()2()(=(,,)6)u u t h u t h u t hf h F h F u t h t u h f G o h ϕ+-=-+++++()()()()()3222122222312222211()()()()()()=(26915.0141(,,)()(215(9).11,2)2,1,6)u u t h u t h N c c t u h c c t u u t h u t hf h F h Ff G o h c c f hc F h c o h h G αϕϕαα+-+-==++=+==++=++++当时令，式化为比较与，则导出方程组这里有三个未知数，只有两个方程,可以得到含一个参数的解族可以得到无:穷多个二级方法两个较常见的二级方法是：()122111,,1221(,)(,(,))217m m m m m m m m c c u u h f t u f t h u hf t u α+===-=+++⎡⎤⎣⎦导出的算法是()()23,914,:N h =当时要求与的所有次数不超过阶的项均相等导出下列决定参数方程组()12322332222332232112181316c c c c c c c c b ααααα++=⎧⎪⎪+=⎪⎪⎨+=⎪⎪⎪=⎪⎩,,.(18)含两个参数的解族所对应的均为三阶方法下面是两个常见的三阶算法()1232233222232322223332(,,)()()1[2()]1411(9)()()()()()=(,,()2))6(2u u u t h u t hf h F h Ff G t u h c c c f h c c o h Fh c b Ff c c u t h u t h t o h u h G ϕαααααϕ+-=+++=+++-+++++++()123233213122I ,0,,,,,44333c c c ααα======计算公式为Runge-Kutta 方法的导出()11312132(3)4(,),1119,3322,33.m m m m m m m m h u u k k k f t u k f t h u hk k f t h He u n hk u +⎧-=+⎪⎪=⎪⎪⎛⎫⎨=++ ⎪⎪⎝⎭⎪⎛⎫⎪=++ ⎪⎪⎝⎭⎩它叫做三阶方法()()()()()123233*********121211II 1,2,636246,2011,22,,,,,.,2m m m m m m m m c c c a b h u u k k k k f t u k f t h Kutta Runge u hk k f t h u h Kutt k h a k α+======⎧-=++⎪⎪=⎪⎨⎛⎫⎪=++ ⎪⎪⎝⎭⎪=+--+⎩此时算式为称为三阶算法这种算法曾是较为流行的方法()()()()112341213243226,11,212211,22,m m m m m m m m m m h u u k k k k k f t u k f t h u hk k f t h u hk k f t h u hk +⎧-=+++⎪⎪=⎪⎪⎛⎫⎪=++ ⎪⎨⎝⎭⎪⎪⎛⎫=++ ⎪⎪⎝⎭⎪=++⎪⎩,类似可以导出带两个参数的四级方法常用者有下列两种算法：()()()()112341213124123338,11,223321,33,m m m m m m m m m m h u u k k k k k f t u k f t h u hk k f t h u hk hk k f t h u hk hk hk +⎧-=+++⎪⎪=⎪⎪⎛⎫⎪=++ ⎪⎨⎝⎭⎪⎪⎛⎫=+-+ ⎪⎪⎝⎭⎪=++-+⎪⎩()22-,-.Runge Kutta Runge Kutta 是在所有方法最常用者称为古典方法。

runge-kutta方法Runge-Kutta方法是求解常微分方程组的一种常用数值方法，它是由德国数学家卡尔·龙格和马丁·库塔（Martin Kutta）发明的。

该方法是一种高阶精度的方法，可以求解一阶或高阶的常微分方程组。

$$ y_{n+1} = y_n + \sum_{i=1}^s b_i k_i $$$y_n$表示第$n$个步骤中的输出，$y_{n+1}$是下一个步骤的输出。

$k_i$表示第$i$个中间变量，$b_i$是权值，它们是根据所采用的方法计算得出的。

我们将在下文中详细介绍这些变量和权值的计算方法。

在数值计算中，我们通常采用自适应步长的方法，这需要我们给出一个误差容限值$\epsilon$，根据误差限定步长大小。

如果当前步长不满足误差要求，我们会将步长缩小再进行计算，或者选择更高阶的方法进行计算，以提高精度。

接下来我们将详细介绍如何使用Runge-Kutta方法求解常微分方程组。

1. Runge-Kutta方法的中间变量$k_i$的计算我们需要计算中间变量$k_i$的值。

它们的计算公式如下：$$ k_1 = hf(t_n, y_n) $$$$ k_2 = hf(t_n + c_2 h, y_n + a_{21}k_1) $$$$ ... $$$$ k_s = hf(t_n + c_s h, y_n + \sum_{j=1}^{s-1}a_{sj}k_j) $$$h$表示步长大小，$t_n$和$y_n$分别表示当前时间和状态，$c_i$和$a_{ij}$是与所采用的方法相关的常数。

$$ c_2 = \frac{1}{2} , a_{21} = \frac{1}{2}$$- 改进的Runge-Kutta方法（RK4）：根据具体选择的方法，我们可以计算出$k_i$的值。

$b_{ij}$也是与所采用的方法相关的常数，对于不同的方法而言，它们的计算公式不同。

根据计算出的中间变量$k_i$和权值$b_i$，我们可以进行Runge-Kutta方法的迭代，求解常微分方程组。

《微分方程数值解法》课程实验报告1 实验内容要求：用经典三级三阶R-K 格式求解微分方程初值问题，并给出误差分析。

2算法描述先用C 的方法写出一个算法动态库,里面封装龙格库塔算法函数采用迭代原理，用递归实现,生成并模块化导出dll,lib 文件,两个文件中均包含了该函数偏移地址,在在源文件中隐式链接该库里的龙格库塔函数,从而得出结果.3 实验数据与实验结果 ⎩⎨⎧=-+-==1|1'0x y x y y1.0)1,0(=∈h x4 程序代码清单：Win32动态库Algorithm.dll 代码: #include <stdio.h> #include <math.h>#define e 2.718281828459045double x[11]={0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}; double y[12]={1.0,0}; int x1=0,x2=1; int count=1;void RungeKutta(double h){ double K1,K2,K3; double ExactSolution; double Error; if(count==12){ return; } K1=-y[count-1]+x[count-1]+1; K2=-(y[count-1]+h*K1/3)+(x[count-1]+h*K1/3)+1; K3=-(y[count-1]+2*h*K2/3)+(x[count-1]+2*h/3)+1; y[count]=y[count-1]+h*(K1+3*K3)/4; ExactSolution=x[count-1]+pow(e,-x[count-1]); Error=y[count-1]-ExactSolution; printf("%lf %lf %lf %lf\n",x[count-1],ExactSolution,y[count-1],Error);count++; RungeKutta(h); }模块化导出文件Algorithm.def代码:LIBRARY AlgorithmEXPORTSRungeKutta @1入口函数实现功能代码:#include <stdio.h>#pragma comment(lib,"../lib/Algorithm.lib")int main(){double h;printf("用三级三阶龙贝格库塔方法解微分方程y'=x-y+1,x属于区间(0,1)\n");printf("请输入系数h的值:\n");scanf("%lf",&h);printf("-----------------------------------------------------\nx 精确解 R-k解yn 误差:\n");RungeKutta(h);return 0;}5:运行结果:1 实验内容要求：用三阶Admas 预报修正格式求解差分方程初边值问题。

具体步骤：
（1）以四阶实时连续RK公式（1-1）为基础，取，首先按右端连续函数预测出并保存条件函数的值构成插值
函数公式（1-2）；
（2）判别的符号，若，则在区间
上出现间断，转步（3）进行间断处理，若，则
表示在该区间未出现间断，转步（5）；
（3）调用插值公式计算间断点，算法采用混合反置连分公式（1-3）；
（4）根据在中的位置，确立不同的三阶龙格库塔法公式：比如，将间断区间分为两部分[tn,]和
[]，在[tn,]中间用经典三阶龙格库塔法求解，然
后转换数学模型，在[]用新构造的三阶龙格库塔法
求解；
（5）继续按照原步长积分，转（1）直至结束。

三阶龙格库塔法改变：
三阶龙格库塔系数公式：
解释第四步：
《1》这个算法的目的就是遇到间断后，改变龙格库塔公式去计算，为了达到更好的精度
《2》怎样去改变龙格库塔公式：在上面中三阶龙格库塔系数公式已给出，得到,后得到别的系数, 之后将这些系数代入得到新的三阶龙格库塔法公式。

《3》怎样得到系数,：遇到间断点后，比如间断点为t*，间断区间为[tn,tn+1]。

根据t*在间断区间[tn,tn+1]的位置,比如那么，等于的任意一个值。

《4》等到遇到了下一个间断点，继续改变方程。

公式
四阶实时连续RK公式为：（1-1）：
公式（1-2）：
公式（1-3）：
经典三阶龙格库塔法公式：
这个模型可以求解一般间断模型：
比如：
例3：以小球自由落体为数学模型，该程序要求能实时仿真该小球自由落体弹起又落体的过程。