Runge-Kutta法的历史发展与应用

Runge-Kutta法的历史发展与应用

摘要Runge-Kutta法是极其重要的常微分方程数值解法，本文仅就其起源及发展脉络加以简要研究。对Runge、Heun以及Kutta等人的贡献做出适当评述，指出Runge-Kutta 方法起源于Euler折线法。同时对Runge-Kutta法的应用做简要研究。

关键词 Euler折线法 标准四阶Runge-Kutta法 应用

一、发展历史[1]

1.1 Euler折线法

在微分方程研究之初，瑞士数学家L.Euler(1707.4—1783.9)做出了开创性的工作。他和其他一些数学家在解决力学、物理学问题的过程中创立了微分方程这门学科。在常微分方程方面，Euler 在1743年发表的论文中，用代换y ekx给出了任意阶常系数线性微分方程的古典解法，最早引入了“通解”和“特解”的概念。

1768年，Euler在其有关月球运行理论的著作中，创立了广泛用于求初值问题

y f(x,y), x0 x X (1.1) y(x) a (1.2) 0

的数值解的方法，次年又把它推广到二阶方程。欧拉的想法如下：我们选择h 0，然后在x0 x x0 h情况下用解函数的切线

l(x) y0 (x x0)f(x0,y0)

代替解函数。这样对于点

x1 x0 h

就得到

y1 y0 hf(x0,y0)。

在(x1,y1)重复如上的程序再次计算新的方向就会得到所谓的递推公式：

xm 1 xm h, ym 1 ym hf(xm,ym),

Runge-Kutta法的历史发展与应用

这就是Euler 方法。通过连接所有这些切线得到的函数被称为Euler 折线。如果我们令h 0， 这些折线就会越来越接近解函数。

Euler 折线法是最早出现的，虽然它亦是常微分方程初值问题的最简单的数值解法， 但它的一些特性和研究方法对于更复杂的方法却具有普遍意义。几十年后，法国数学家A .L .Cauchy(1789.8—1857.5)在历史上首次研究了常微分方程的局部性态。对于给定的初值问题(1.1)和(1.1)，在f(x,y)连续可微的假设下， 他用类似于欧拉折线的方法构造逼近解， 利用微分中值定理估计逼近解之间差的上界，严格证明了以x0为中心的一个小区间上逼近解收敛， 其极限函数即为所提问题的解。同时Cauchy 指出，这种方法也适用于常微分方程组，所以欧拉方法有时又称Euler-Cauchy 折线法。

2.2 Runge-Kutta方法

德国数学家C.D.T.Runge(1856—1927)是数值方法发展史上具有里程碑作用的人物。1895年，他在Hanover发表了关于微分方程数值解法的经典论文《常微分方程数值解法》， 此文成为常微分方程Runge-Kutta方法的发端。此后，Runge结合教学活动积极投身于发展一般的数值分析特别是各种实际应用中的Runge-Kutta方法（严格来说，此方法在Kutta做出工作后才能称作Runge-Kutta方法）。

Runge伟大创造的思想是什么呢？他的灵感来自于初值问题（1.1）和（1.2）与积分问题

x0 h

y(x0 h) y(x0)

x0 f(x)dx（此时f与y无关） （1.3）

之间的对比，显然，等式（1,3）右侧数值积分的精确度决定y(x0 h)的精确度，Runge发现， Euler方法采用的是左矩形公式

x0 h

x0 f(x)dx hf(x0)，

即用高为f(x0)宽为h的矩形代替数值积分， 而这个公式的精确度并不高。因此他说：最好通过插入上述Euler步骤的结果来代替未出现的y值， 把精度更高的中点法则和梯形法则拓展到微分方程。

Runge-Kutta法的历史发展与应用

11M: y1 y0 hf(x0 h,y0 hf(x0,y0))， 22

T: y1 y0 h(f(x0,y0) f(x0 h,y0 hf(x0,y0)))。 其中M和T分别表示用中点法和梯形法算得的数值积分。与他的后继者一样，Runge 用Taylor 展开式表明上述两方法的局部误差是，方法的阶为2。不过他的梦想却是使用具有更高精度的Simpson法则。但是众所周知，1212R M (T M)/3的微小变化往往易产生假象，令人误以为可以获得更高的阶。Taylor 级数展开表明，如果f依赖于y，事实上这个表达形式只是2 阶的。接着Runge发现，通过重复使用Euler步骤对梯形法则做些许调整，会使形式R M (T M)/3成为3 阶方法。Runge还把他的方法及方法的展开式拓展到微分方程组。

1900 年，Runge的同胞K .Heun评论说，Runge获得的上述方法是归纳性的而且是令人费解的，他提出采用“更具一般性”的Gauss方法。于是一般的Gauss求积公式

y1 y0 h bif(x0 cih)，

i 1s

被扩充为

k1 f(x0,y0),

k2 f(x0 c2h,y0 c2hk1), k3 f(x0

c3h,y0 c3hk2),

y1 y0 h biki （1.4）

i 1s

把（1.4）式的右端进行二元Taylor展开后与y(x h)的Taylor 展开式的对应项的系数比较，适当选取参数使方法具有尽可能高的精度。

Runge的另一个同胞W.M.Kutta(1867—1944)，1894到1909年在Munich做助教，在那里受到Runge文章的吸引并在Heun论文的激励下发表了他的研究结果。他认为，为什么不让已经求得的导数值全部进入到新的求值点的计算中呢？

Runge-Kutta法的历史发展与应用

基于这样的想法，Heun格式就被Kutta代替为如下格式：

k1 f(x0,y0),

k2 f(x0 c2h,y0 a21hk1),

k3 f(x0 c3h,y0 a31hk1 a32hk2),

k4 f(x0 c4h,y0 a41hk1

a42hk2 a43hk3)

y1 y0 h biki

i 1s

(其中s称为级)这个格式在满足所需的阶条件上能够允许更多的自由度。在古典的Runge-Kutta方法中，对系数的选择极大地取决于由这些系数构成的方法是否方便进行桌上计算，而所谓的古典方法是指在前计算机时代得到的方法。而这些方法对于在自动计算机上使用则未必是最适合的方法。

显然Runge-Kutta方法是一种特殊的单步方法，事实上这个方法可以看作在(xm,xm 1)上取若干条积分曲线的若干个点的切线斜率，再进行一次（或多次）算术（或加权）平均后产生的新斜率，再按这个斜率从(xm,ym)出发，以直线代曲线向前推进一步的过程。与Taylor展开方法相比，Runge-Kutta方法不用增加微商f(x,y)的次数就可以得到较高的阶。

在Runge于1895年提出Runge-Kutta方法的雏形的时候，他给出的方法是2级2阶和4级3阶的；Heun在1900年获得的两个方法分别是3级3阶和8级4阶的；而Kutta在1901年得到的两个方法则分别是4级4阶和6级5阶的。因为Runge是提出这个方法雏形的人而Kutta则是得到了此方法在前计算时代最高的阶，所以这个方法被称为Runge-Kutta方法。到1925 年，另外一位数学家Nyström也得到了6级5阶的方法。此后直到计算机诞生也未产生新的重要成果，所以Nyström在1925年得到6级5阶方法的论文有时也被称为古典Runge-Kutta方法与现代Runge-Kutta方法的分水岭。

2.3现代Runge-Kutta方法

Runge-Kutta方法在创立之初并未达到完善，后来又经过大量数学家的共同努力才日趋成熟。20世纪60年代，新西兰数学家J.C.Butcher（1933—）给出了现代Runge-Kutta方法的一般形式：

Runge-Kutta法的历史发展与应用

s yn 1 yn h biki;

i 1

s ki f(xn cih,yn h aijkj), i,j 1,2,

j 1

s ci aij, i 1,2,,s. j 1 ,s; （1.5）

为了展示（1.5）式中出现的系数，Butcher给出了后来以其名字命名的Butcher点阵：

c1

c2

csa11a21as1

b1a12a22as2b

2a1sa2s assbs

我们分别把s维向量c、b以及s s矩阵A定义为

c c1,c2,,cs , b b1,b2,T,bs , A aij .

显然， 一s级Runge-Kutta方法被其Butcher点阵完全确定，这样，研究

某一Runge-Kutta方法的性质就转化为研究与其相对应的矩阵A的性质。此点阵为研究Runge-Kutta方法的性质提供了方便的途径。在计算机诞生之前，Runge-Kutta方法被禁锢在只能用手进行计算的实际问题上，但是随着计算机的出现，Runge-Kutta方法呈现出新的史无前例的重要性。能够解决的问题变得越来越大、越来越复杂，自动化的误差检测和步长控制不但变得可行而且变得必要了，Runge-Kutta方法的发展不但表现在理论上而且表现在技术上。

除了Runge-Kutta方法在微分方程求解中扮演的传统角色外，人们发现相关

类型的初值问题可以用Runge-Kutta方法或适合于更一般问题的Runge-Kutta方法求解，比如Runge-Kutta方法被应用到了Hamilton系统中。

二、实际应用[2]

2.1 经典四阶Runge-Kutta方法

经典的四阶Runge-Kutta方法是：

Runge-Kutta法的历史发展与应用

h y y (k1 2k2 2k3 k4)n n 16 k1 f(xn,yn)

11 k f(x h,y hk1) 2nn22 11 k f(x h,y hk2)nn 322 k4 f(xn h,yn hk3)

它的局部截断误差Tn 1 O(h5)，所以为四阶方法，这是最常用的四阶Runge-Kutta方法，数学库中都有用此方法求解初值问题的软件，这种方法的优点是精度较高，缺点是每步要计算四个函数值，计算量较大。

2.2 两种群共存模型的Runge-Kutta法解答

2.2.1 问题提出

在自然界中处于同一环境下有多个种群生存，它们之间存在着一定的关系，它们之间可能是食饵与捕食者，或是相互依存，或是相互竞争或其它的关系。简单的食饵与捕食者模型是由两位数学生态学的先驱者A.J.Lotka与意大利数学家V.Volterra各自建立，因此模型被命名为Lotka-Volterra模型。在该文献中，在Lotka-Volterra模型的基础上建立考虑捕食者另有食物来源，并且把被捕食种群作为生存资源，限制捕食种群增长的两种群共存的数学模型，并对模型进行分析和模型解释。两种群相互共生存的现象是普遍存在的,研究其种群的数量演变规律具有重要意义，它有助于我们对种群变化情况的了解和正确解释，可以帮助我们合理解决这些生态问题。下面是该文献的两种群相互共存的微分方程模型。

2.2.2 模型建立与分析

①假设有种群A和种群B，种群A是种群B的食物来源，但不是唯一的食物来源，种群B还另有食物来源。

②假设M1是环境容许的种群A的最大数量，M2是环境容许的种群B的最大数量。

③假设种群A能独立生存，设其固有增长率为 1，由于种群A是种群B的食物来源，故设 1 0为单位数量的种群B（相对M2而言）掠取 1倍的单位种群B的量（相对M1而言），还考虑到种群A的自身的阻滞作用，因此设：

Runge-Kutta法的历史发展与应用

y2y1 y1(t) 1y1 1 1 ，y1(t)表示在t时刻种群A的数量。 M2M1

④因为种群B另有食物来源，设其增长率为 2，假设种群A作为生存资源限制种群B的增长，设对其增长率的影响率为 2y2，（ 2 0）（表示种群A对y1

种群B的影响，若种群A充足，y2小，则对种群B的增长影响较小；若种群Ay1

不充足，y2大，则对种群B的增长影响较大），不考虑种群B自身内部竞争，因y1

y 此设y2(t) 2y2 1 22 ，y2(t)表示在t时刻种群B的数量。 y1

根据上面的假设可以得到相互共存的两种群的微分方程模型。

dy1 y2y1 y1 11 1dtMM 21 （2.1） dy2 y 1 y2 22 2 dty 1

其初始条件为y1(0) y10，y2(0) y20。

将（2.1）可归结为如下形式的一阶常微分方程组

dy1 F1(t,y1,y2) dt （2.2） dy 2 F(t,y,y)212 dt

因为此方程组没有解析解，故采用标准的四阶Runge-Kutta方法求出数值解，求解公式如下：

Runge-Kutta法的历史发展与应用

h y y (k11 2k21 2k31 k41)1,n 1,n 16 h y y (k12 2k22 2k32 k42)2,n 2,n 16 k1j Fj(tn,y1,n,y2,n)

k F t 1h,y 1hk,y 1hk j 1,2

j n1,n112,n12 2j222 k F t 1h,y 1hk,y 1hk

j n1,n212,n22 3j222 11 k Ft h,y hk,y hkj n1,n312,n32 4j22

F1y 2 F y 1 1 12 11,1 111

y1M2 M1 y2M2

F2 2 2y2 F22 2 2y2 , 2 y1y12 y2y12因为

存在且连续，由定理可知采用标准的四阶Runge-Kutta方法是收敛的。

为了对微分方程模型进行数值计算，假设初始值为

y1(0) 100,y2(0) 200，

方程组（2.1）中 1 06.， 2 0.4， 1 0.5， 2 0.2，M1 1000，M2 3000，用Matlab编程如下：

function dq=myfile_2(t,y)

dq=[0.5*y(1)*(l-0.5*y(2)/3000-y(l)/1000);0.4*y(2)*(l-0.2*y(2)/y(l))];

ode45(@myfile_2,[0,100],[100;200]),grid

Runge-Kutta法的历史发展与应用

于是得y1(t)、y2(t)的图形1如下：

为了体现模型中所设参数 1， 2， 1， 2在生态学上的意义，我们再设一组参数 1 0.5， 2 0.3， 1 0.2， 2 1.5，用Matlab编程如下： function dq=myfile_3(t,y)

dq=[0.5*y(1)*(l-0.2*y(2)/3000-y(l)/1000);0.3*y(2)*(l-1.5*y(2)/y(l))];

ode45(@myfile_3,[0,100],[100:200]),grid

Runge-Kutta法的历史发展与应用

于是得图形2：

比较上述两个图形可以看出种群A与种群B的相互影响，种群B的数量与对其增长率的影响为 2y2， 2 0 有关，系数 2可以反映种群A对种群B的影y1

响程度，系数 2越小，种群A对种群B的增长影响较小，系数 2越大，种群A对种群B的增长影响较大，种群B的数量随着种群A的数量的变化而出现不同的情况变化。系数 1可以反映种群B对种群A的影响程度，系数 1越小，种群B对种群A的增长影响较小，系数 1越大，种群B对种群A的增长影响较大，从而种群A的数量变化出现不同的情况。从上述图形分析比较可以知道，通过改变种群A的数量可以使种群B的数量发生变化，但最终两种群趋向于稳

定而共存。

三、总结

通过在网上查阅资料，本文介绍了常微分方程初值解法——Runge-Kutta法的发展历程和其中的主要人物，并将其方法应用到求解建立的常微分方程模型，通过数值计算作出图形，并对模型作出分析，从而帮助我们分析解决实际问题和实际现象的解释，有利于我们采取实际有效的措施解决问题。

由于我们在解决现

Runge-Kutta法的历史发展与应用

实中的许多问题往往通过建立数学模型来解决，而建立的模型很多是微分方程模型，并且很难或无法求出解析解，因此，运用数值解法对于求解常微分方程模型的数值解有着重要的意义，研究微分方程的数值求解具有很大的潜力，是值得我们去探讨和研究的。

参考文献：

[1]林立军，常微分方程数值解法——Runge_Kutta法的历史浅析[J]，辽宁师范大学学报（自然科学版），2003

[2]秦军，Rung-Kutta法在求解微分方程模型中的应用[J]，2010

[3]李庆扬,王能超,易大义.数值分析[M]，清华大学出版社，2001