数值分析_数值计算小论文

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年发表的论文中，用代换kx y e =给出了任意阶常系数线性微分方程的古典解法，最早引入了“通解”和“特解”的概念。

1768年，Euler 在其有关月球运行理论的著作中，创立了广泛用于求初值问题
00
(,), (1.1)() (1.2)y f x y x x X y x a '=<≤⎧⎨=⎩ 的数值解的方法，次年又把它推广到二阶方程。

欧拉的想法如下：我们选择0h >，然后在00x x x h ≤≤+情况下用解函数的切线
0000()()(,)l x y x x f x y =+-
代替解函数。

这样对于点
10x x h =+
就得到
1000(,)y y hf x y =+。

在11(,)x y 重复如上的程序再次计算新的方向就会得到所谓的递推公式：
11, (,),m m m m m m x x h y y hf x y ++=+=+
这就是Euler 方法。

通过连接所有这些切线得到的函数被称为Euler 折线。

如果我们令0h →， 这些折线就会越来越接近解函数。

Euler 折线法是最早出现的，虽然它亦是常微分方程初值问题的最简单的数值解法， 但它的一些特性和研究方法对于更复杂的方法却具有普遍意义。

几十年后，法国数学家A .L .Cauchy(1789.8—1857.5)在历史上首次研究了常微分方程的局部性态。

对于给定的初值问题(1.1)和(1.1)，在(,)f x y 连续可微的假设下， 他用类似于欧拉折线的方法构造逼近解， 利用微分中值定理估计逼近解之间差的上界，严格证明了以0x 为中心的一个小区间上逼近解收敛， 其极限函数即为所提问题的解。

同时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）与积分问题
0000()()()x h
x y x h y x f x dx f y ++=+
⎰（此时与无关）
（1.3） 之间的对比，显然，等式（1,3）右侧数值积分的精确度决定0()y x h +的精确度，Runge 发现， Euler 方法采用的是左矩形公式
000()()x h
x f x dx hf x +≈⎰，
即用高为0()f x 宽为h 的矩形代替数值积分， 而这个公式的精确度并不高。

因此他说：最好通过插入上述Euler 步骤的结果来代替未出现的y 值， 把精度更高的中点法则和梯形法则拓展到微分方程。

10000011:(,(,))22
M y y hf x h y hf x y ==+++， 1000000011:((,)(,(,)))22T y y h f x y f x h y hf x y ==++++。

其中M 和T 分别表示用中点法和梯形法算得的数值积分。

与他的后继者一样，Runge 用Taylor 展开式表明上述两方法的局部误差是，方法的阶为2。

不过他的梦想却是使用具有更高精度的Simpson 法则。

但是众所周知，()/3R M T M =+-的微小变化往往易产生假象，令人误以为可以获得更高的阶。

Taylor 级数展开表明，如果f 依赖于y ，事实上这个表达形式只是2 阶的。

接着Runge 发现，通过重复使用Euler 步骤对梯形法则做些许调整，会使形式()/3R M T M =+-'成为3 阶方法。

Runge 还把他的方法及方法的展开式拓展到微分方程组。

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

于是一般的Gauss 求积公式
1001()s
i i i y y h b f x c h ==++∑，
被扩充为
100202021303032(,),
(,),(,),
k f x y k f x c h y c hk k f x c h y c hk ==++=++ 101s
i i i y y h b k ==+∑ （1.4）
把（1.4）式的右端进行二元Taylor 展开后与()y x h +的Taylor 展开式的对应项的系数比较，适当选取参数使方法具有尽可能高的精度。

Runge 的另一个同胞W.M.Kutta(1867—1944)，1894到1909年在Munich 做助教，在那里受到Runge 文章的吸引并在Heun 论文的激励下发表了他的研究结果。

他认为，为什么不让已经求得的导数值全部进入到新的求值点的计算中呢？
基于这样的想法，Heun 格式就被Kutta 代替为如下格式：
100202021130303113224040411422433(,),
(,),
(,),
(,)
k f x y k f x c h y a hk k f x c h y a hk a hk k f x c h y a hk a hk a hk ==++=+++=++++
101 s i i i y y h b k ==+∑
(其中s 称为级)这个格式在满足所需的阶条件上能够允许更多的自由度。

在古典的Runge-Kutta 方法中，对系数的选择极大地取决于由这些系数构成的方法是否方便进行桌上计算，而所谓的古典方法是指在前计算机时代得到的方法。

而这些方法对于在自动计算机上使用则未必是最适合的方法。

显然Runge-Kutta 方法是一种特殊的单步方法，事实上这个方法可以看作在1(,)m m x x +上取若干条积分曲线的若干个点的切线斜率，再进行一次（或多次）算术（或加权）平均后产生的新斜率，再按这个斜率从(,)m m x y 出发，以直线代曲线向前推进一步的过程。

与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 方法的一般形式：
1111;
(,), ,1,2,
,;, 1,2,,.s n n i i i s i n i n ij j j s i ij j y y h b k k f x c h y h a k i j s c a i s +===⎧=+⎪⎪
⎪⎪=++=⎨⎪
⎪==⎪⎪⎩∑∑∑ （1.5）
为了展示（1.5）式中出现的系数，Butcher 给出了后来以其名字命名的Butcher 点阵：
1
111212
212221
212s s s s s ss s
c a a a c a a a c a a a b b b 我们分别把s 维向量c 、b 以及s s ⨯矩阵A 定义为
[][]1212,,,, ,,,, .T s s ij c c c c b b b b A a ⎡⎤===⎣⎦
显然， 一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 方法是：
112341213243(22)6(,)
11(,)2211(,)22(,)
n n n n n n n n n n h y y k k k k k f x y k f x h y hk k f x h y hk k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ 它的局部截断误差51()n T O h +=，所以为四阶方法，这是最常用的四阶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 还另有食物来源。

②假设1M 是环境容许的种群A 的最大数量，2M 是环境容许的种群B 的最大数量。

③假设种群A 能独立生存，设其固有增长率为1μ，由于种群A 是种群B 的食物来源，故设10λ>为单位数量的种群B （相对2M 而言）掠取1λ倍的单位种群
B 的量（相对1M 而言），还考虑到种群A 的自身的阻滞作用，因此设：
21111121()1y y y t y M M μλ⎛⎫=-- ⎪⎝
⎭，1()y t 表示在t 时刻种群A 的数量。

④因为种群B 另有食物来源，设其增长率为2μ，假设种群A 作为生存资源限制种群B 的增长，设对其增长率的影响率为221
y y λ-，（20λ>）（表示种群A 对种群B 的影响，若种群A 充足，21
y y 小，则对种群B 的增长影响较小；若种群A 不充足，21
y y 大，则对种群B 的增长影响较大），不考虑种群B 自身内部竞争，因此设222221()1y y t y y μλ⎛⎫=- ⎪⎝
⎭，2()y t 表示在t 时刻种群B 的数量。

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

1211112122222111dy y y y dt M M dy y y dt y μλμλ⎧⎛⎫=--⎪ ⎪⎪⎝⎭⎨⎛⎫⎪=- ⎪⎪⎝⎭⎩
（2.1） 其初始条件为110(0)y y =，220(0)y y =。

将（2.1）可归结为如下形式的一阶常微分方程组
11122212(,,)(,,)dy F t y y dt dy F t y y dt
⎧=⎪⎪⎨⎪=⎪⎩ （2.2） 因为此方程组没有解析解，故采用标准的四阶Runge-Kutta 方法求出数值解，求解公式如下：
1,11,112131412,12,1222324211,2,21,112,1231,212,2241,312,(22)6(22)6(,,)
111,,222111,,22211,,22n n n n j j n n n j j n n n j j n n n j j n n n h y y k k k k h y y k k k k k F t y y k F t h y hk y hk k F t h y hk y hk k F t h y hk y ++=++++=++++=⎛⎫=+++ ⎪⎝⎭⎛⎫
=+++ ⎪⎝⎭=+++()321,2j hk ⎧⎪⎪⎪⎪⎪⎪⎪=⎨⎪⎪⎪⎪⎪⎛⎫⎪ ⎪⎪⎝⎭⎩
因为 12111111
111212222222222222112121,2,F y F y y M M y M F y F y y y y y μλμλμλμλμλμ⎛⎫∂∂=--=- ⎪∂∂⎝
⎭∂∂==-∂∂
存在且连续，由定理可知采用标准的四阶Runge-Kutta 方法是收敛的。

为了对微分方程模型进行数值计算，假设初始值为
12(0)100,(0)200y y ==，
方程组（2.1）中10.6μ=，20.4μ=，10.5λ=，20.2λ=，11000M =，23000M =，用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
于是得1()y t 、2()y t 的图形1如下：
为了体现模型中所设参数1μ，2μ，1λ，2λ在生态学上的意义，我们再设一组参数10.5μ=，20.3μ=，10.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
于是得图形2：
比较上述两个图形可以看出种群A 与种群B 的相互影响，种群B 的数量与对
其增长率的影响为221
y y λ-，()20λ>有关，系数2λ可以反映种群A 对种群B 的影响程度，系数2λ越小，种群A 对种群B 的增长影响较小，系数2λ越大，种群A 对种群B 的增长影响较大，种群B 的数量随着种群A 的数量的变化而出现不同的情况变化。

系数1λ可以反映种群B 对种群A 的影响程度，系数1λ越小，种群B 对种群A 的增长影响较小，系数1λ越大，种群B 对种群A 的增长影响较大，从而种群A 的数量变化出现不同的情况。

从上述图形分析比较可以知道，通过改变种群A 的数量可以使种群B 的数量发生变化，但最终两种群趋向于稳
定而共存。

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

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

参考文献：
[1]林立军，常微分方程数值解法——Runge_Kutta 法的历史浅析[J]，辽宁师范大学学报（自然科学版），2003
[2]秦军，Rung-Kutta 法在求解微分方程模型中的应用[J]，2010
[3]李庆扬,王能超,易大义.数值分析[M]，清华大学出版社，2001。