最后更新: 2026-06-28 21:02 查看: 2 次纠错评论(0) 课件开VIP

龙格-库塔法

8.3 龙格—库塔法

由上节知道，截断误差的阶是衡量一个方法精度高低的主要依据。能否用提高截断误差阶来提高方法的精确度呢？回答是肯定的．比较有代表性的是泰勒级数法和龙格一库塔（Runge－Kutta）法，下面介绍经典的龙格－库塔法．关于泰勒级数法，感兴趣的读者可以参考相关书籍。

从理论上讲，只要函数 $y=f(x)$ 在区间 $[a, b]$ 上充分光滑，那么它的各阶导数值 $y^{(k)}\left(x_n\right)$ 与函数 $y(x)$ 在区间 $[a, b]$ 上某些点的值就相互有联系，就是说，函数值可用各阶导数值近似地表示出来，反之，各阶导数值也可用函数在一些点上值的线性组合近似地表示出来．事实上，欧拉法和改进的欧拉法也可以看成是导数值用函数值的线性组合表示的特例，例如，改进欧拉法可以写成

$$
\left\{\begin{array}{l}
y_{n+1}=y_n+\frac{h}{2}\left(k_1+k_2\right), \\
k_1=f\left(x_n, y_n\right), \\
k_2=f\left(x_n+h, y_n+h k_1\right) .
\end{array}\right.
$$

此公式也可称为二阶龙格一库塔公式。
为了导出龙格一库塔法的一般公式，我们取如下的线性组合形式：

$$
\begin{equation*}
y_{n+1}=y_n+h \sum_{i=1}^\gamma w_i k_i, \tag{8.31}
\end{equation*}
$$

其中

$$
k_i=f\left(x_n+b_i h, y_n+h \sum_{j=1}^{i-1} a_{i j} k_j\right),
$$

即

$$
\left\{\begin{array}{l}
k_1=f\left(x_n, y_n\right), \\
k_2=f\left(x_n+b_2 h, y_n+h a_{21} k_1\right), \\
k_3=f\left(x_n+b_3 h, y a_{31} k_1+h a_{32} k_2\right), \\
\quad \quad \cdots \cdots
\end{array}\right.
$$

而 $w_1, w_2, \cdots, w_v ; b_1=0, b_2, b_3, \cdots, b_v ; a_{21}, a_{31}, \cdots, a_{v, v-1}$ 中除 $b_1=0$ 外均为待定系数。适当选取这些系数，使得局部截断误差的阶尽可能高即可．显然，当 $\gamma=1$ 时，式（8．31）就是欧拉公式．

下面我们先导出 $\gamma=2$ 时的公式．将 $k_1, k_2$ 在同一点 $\left(x_n, y_n\right)$ 泰勒展开，则有

$$
\begin{align*}
k_1 & =f\left(x_n, y_n\right) \\
k_2 & =f\left(x_n+b_2 h, y_n+h a_{21} k_1\right) \\
& =f\left(x_n, y_n\right)+h\left(\left.b_2 \frac{\partial f}{\partial x}\right|_{\left(x_n, y_n\right)}+\left.a_{21} k_1 \frac{\partial f}{\partial y}\right|_{\left(x_n, y_n\right)}\right)+O\left(h^2\right) \tag{8.32}
\end{align*}
$$

将（8．32）代入（8．31）并与 $y\left(x_n+h\right)$ 在 $x_n$ 点的泰勒展开式

$$
\begin{aligned}
y\left(x_n+h\right)= & y\left(x_n\right)+h f\left(x_n, y_n\right)+\frac{h^2}{2}\left(\left.\frac{\partial f}{\partial x}\right|_{\left(x_n, y_n\right)}+\left.f\left(x_n, y_n\right) \frac{\partial f}{\partial y}\right|_{\left(x_n, y_n\right)}\right) \\
& +\frac{h^3}{3!}\left(\left.\frac{\partial^2 f}{\partial x^2}\right|_{\left(x_n, y_n\right)}+\left.2 f\left(x_n, y_n\right) \frac{\partial^2 f}{\partial x \partial y}\right|_{\left(x_n, y_n\right)}+\left.f^2\left(x_n, y_n\right) \frac{\partial^2 f}{\partial y^2}\right|_{\left(x_n, y_n\right)}\right. \\
& \left.+\left.\frac{\partial f}{\partial y}\right|_{\left(x_n, y_n\right)}\left(\left.\frac{\partial f}{\partial x}\right|_{\left(x_n, y_n\right)}+\left.f\left(x_n, y_n\right) \frac{\partial f}{\partial y}\right|_{\left(x_n, y_n\right)}\right)\right)+\cdots
\end{aligned}
$$

逐项比较，若令 $h, h^2$ 项的系数相等，便得到

$$
w_1+w_2=1, \quad w_2 b_2=\frac{1}{2}, \quad w_2 a_{21}=\frac{1}{2}
$$

把 $b_2$ 作为自由参数来确定 $w_1$ 和 $w_2$ ，如取 $b_2=1$ ，则 $w_1=w_2=\frac{1}{2}, a_{21}=1$ ，这时（8．32）正好就是改进的欧拉法，截断误差的阶为 $O\left(h^3\right)$ ．

这样一来，二阶龙格－库塔法的具体迭代关系为

$$
\left\{\begin{array}{l}
y_{n+1}=y_n+h k_2,  \tag{8.33}\\
k_1=f\left(x_n, y_n\right), \\
k_2=f\left(x_n+\frac{h}{2}, y_n+\frac{h}{2} k_1\right) .
\end{array}\right.
$$

公式（8．33）称为二阶龙格一库塔公式。
对于 $\gamma=3$ 的情形，我们也可以完全仿上述方法推导出三阶龙格一库塔公式．这时参数满足下列条件

$$
\left\{\begin{array} { l } 
{ w _ { 1 } + w _ { 2 } + w _ { 3 } = 1 , }  \tag{8.34}\\
{ b _ { 2 } w _ { 2 } + b _ { 3 } w _ { 3 } = \frac { 1 } { 2 } , } \\
{ a _ { 2 1 } w _ { 2 } + ( a _ { 3 1 } + a _ { 3 2 } ) w _ { 3 } = \frac { 1 } { 2 } , } \\
{ b _ { 2 } ^ { 2 } w _ { 2 } + b _ { 3 } ^ { 2 } w _ { 3 } = \frac { 1 } { 3 } . }
\end{array} \left\{\begin{array}{l}
b_2 a_{21} w_2+b_2\left(a_{31}+a_{32}\right) w_3=\frac{1}{3}, \\
a_{21}^2 w_2+\left(a_{31}+a_{32}\right)^2 w_3=\frac{1}{3}, \\
b_2 a_{32} w_2=\frac{1}{6}, \\
a_{21} a_{32} w_3=\frac{1}{6} .
\end{array}\right.\right.
$$

（8．34）比较简单的一组解为

$$
b_2=\frac{1}{2}, \quad b_3=1, \quad a_{21}=\frac{1}{2}, \quad a_{31}=-1, \quad a_{32}=2, \quad w_1=\frac{1}{6}, \quad w_2=\frac{4}{6}, \quad w_3=\frac{1}{6} .
$$

将它代入（8．31）得

$$
\left\{\begin{array}{l}
y_{n+1}=y_n+\frac{h}{6}\left(k_1+4 k_2+k_3\right),  \tag{8.35}\\
k_1=f\left(x_n, y_n\right), \\
k_2=f\left(x_n+\frac{h}{2}, y_n+\frac{k_1}{2}\right), \\
k_3=f\left(x_n+h, y_n-k_1+2 k_2\right) .
\end{array}\right.
$$

（8．35）称为三阶龙格一库塔公式．当然，参数的不同选取公式（8．35）就有不一样的形式，但它们的截断误差阶都是 $O\left(h^4\right)$ 。

通常人们所说的龙格一库塔法是指四阶的．我们可以仿照二阶的情形推导出此公式，不过太繁杂，此处从略，常用的四阶龙格一库塔公式是

$$
\left\{\begin{array}{l}
y_{n+1}=y_n+\frac{h}{6}\left(k_1+2 k_2+2 k_3+k_4\right),  \tag{8.36}\\
k_1=f\left(x_n, y_n\right), \\
k_2=f\left(x_n+\frac{1}{2} h, y_n+\frac{h}{2} k_1\right), \\
k_3=f\left(x_n+\frac{1}{2} h, y_n+\frac{h}{2} k_2\right), \\
k_4=f\left(x_n+h, y_n+h k_3\right) .
\end{array}\right.
$$

公式（8．36）的截断误差阶为 $O\left(h^5\right)$ ．
例8．6 用二阶方法和四阶龙格一库塔法求下面微分方程

$$
\left\{\begin{array}{l}
y^{\prime}=-y \ln y, \quad x \in[0,1], \\
y(0)=\frac{1}{2}
\end{array}\right.
$$

在步长为 $h=1 / 2$ ，条件下 $t=1$ 的值．
令 $h=1 / 2$ ，要计算出 $t=1$ 的值，每种方法都需要两个步骤．对于二阶龙格一库塔法，具体计算如下：

第一步

$$
\begin{aligned}
& \bar{f}=f\left(t_0, y_0\right)=f(0,1 / 2)=-(1 / 2) \ln (1 / 2)=0.3465735903, \\
& \bar{y}=y_0+(2 / 3) h \bar{f}=0.5+(1 / 3) \bar{f}=0.6155245301,
\end{aligned}
$$

$$
\begin{aligned}
y_1 & =y_0+\frac{1 / 2}{4}\left[\bar{f}+3 f\left(t_0+(2 / 3) h, \bar{y}\right)\right] \\
& =0.5+0.125 \times[0.3465735903+3 \times 0.2987020395] \\
& =0.6553349636
\end{aligned}
$$

第二步

$$
\begin{aligned}
\bar{f} & =f\left(t_1, y_1\right)=f(0.5,0.6553349636)=0.276950309, \\
\bar{y} & =y_1+(2 / 3) h \bar{f}=0.6553349636+1 / 3 \times 0.276950309=0.7476517333, \\
y_2 & =y_1+\frac{1 / 2}{4}\left[\bar{f}+f\left(t_1+(2 / 3) h, \bar{y}\right)\right] \\
& =0.6553349636+0.125 \times[0.276950309+3 \times 0.2174305868] \\
& =0.7714902223 .
\end{aligned}
$$

对于四阶龙格一库塔法有
第一步

$$
\begin{aligned}
& k_1=h f\left(t_0, y_0\right)=0.17328679513999 \\
& z_1=y_0+(1 / 2) k_1=0.58664339756999 \\
& k_2=h f\left(t_0+(1 / 2) h, z_1\right)=0.15643965031862, \\
& z_2=y_0+(1 / 2) k_2=0.57821982515931, \\
& k_3=h f\left(t_0+(1 / 2) h, z_2\right)=0.15837474613445, \\
& z_3=y_0+k_3=0.65837474613445, \\
& k_4=h f\left(t_0+h, z_3\right)=0.13759406302779, \\
& y_1=y_0+(1 / 6)\left(k_1+2 k_2+2 k_3+k_4\right)=0.65675160851232 .
\end{aligned}
$$

第二步

$$
\begin{aligned}
& k_1=h f\left(t_1, y_1\right)=0.13806541027436, \\
& z_1=y_1+(1 / 2) k_1=0.72578431364950, \\
& k_2=h f\left(t_1+(1 / 2) h, z_1\right)=0.11630780609087, \\
& z_2=y_1+(1 / 2) k_2=0.71490551155775, \\
& k_3=h f\left(t_1+(1 / 2) h, z_2\right)=0.11996289517416, \\
& z_3=y_1+k_3=0.77671450368648, \\
& k_4=h f\left(t_1+h, z_3\right)=9.8131054204024 \mathrm{E}-2, \\
& y_2=y_1+(1 / 6)\left(k_1+2 k_2+2 k_3+k_4\right)=0.77487458634706 .
\end{aligned}
$$

在 $t=1$ 时，解的精确值是 $y(1)=0.7749206845$ ，因此我们可以看到，即使用比较粗的网格，四阶方法的效果也很好，而二阶方法的效果则不佳——这也是高阶方法的一个好处．通常来说，如果我们对一个问题应用某种方法，在相同网格条件下，

高阶的方法一般会有更好的准确性。
例 8.7 在区间 $[0,1.5]$ 上，取 $h=0.1$ ，用四阶龙格－库塔法求解下面的微分方程：

$$
\left\{\begin{array}{l}
y^{\prime}=y-\frac{2 x}{y}, \\
y(0)=1 .
\end{array}\right.
$$

解 取 $h=0.1$ ，由（8．36）得

$$
\left\{\begin{array}{l}
k_1=x_n^2-y_n, \\
k_2=\left(x_n+0.05\right)^2-\left(y_n+0.05 k_1\right), \\
k_3=\left(x_n+0.05\right)^2-\left(y_n+0.05 k_2\right), \\
k_4=\left(x_n+0.1\right)^2-\left(y_n+0.1 k_3\right) .
\end{array}\right.
$$

把初始条件 $x_0=0, y_0=1$ ，代入上式得 $k_1=-1, k_2=-0.9475, k_3=-0.9501, k_4=0.8950$ ．将这些 $k$ 值代入（8．36）得

$$
y_1=1+\frac{0.1}{6}[-1+2(-0.9475-0.9501)-0.8950] \approx 0.9052 .
$$

重复上述步骤可算出 $y_2, y_3, \cdots, y_{10}$ 等。
值得注意的是，虽然龙格－库塔法比较高效，但是该方法有一个不足：与其他方法相比，它们需要对函数 $f(x)$ 进行计算．对于像我们的示例这样的简单方程，这没什么大不了的，但是对于较大的系统，这是一个费时的问题。为了获得与龙格一库塔法一样高精度的方法，通常采取多步法．

最后简单说明一下单步法的收敛性．所谓单步法，就是在计算 $y_{n+1}$ 时只用到它前一步的信息 $y_n$ ．显式单步法的共同特征是，它们都是将 $y_n$ 加上某种形式的增量得出 $y_{n+1}$ ，其计算公式为

$$
y_{n+1}=y_n+h \varphi\left(x_n, y_n, h\right)
$$

$\varphi\left(x_n, y_n, h\right)$ 称为增量函数，不同的单步法，对应不同的增量函数．这里不加证明地引用如下定理．

定理 8.3 若单步法的增量函数 $\varphi(x, y, h)$ 满足条件 $|\varphi(x, y, h)-\varphi(x, \bar{y}, h)| \leqslant L_{\varphi}|y-\bar{y}|$（利普希茨条件），且初值 $y_0$ 是准确的，即 $y_0=y\left(x_0\right)$ ，则该单步法是收敛的．

科数题库(单机版) 会议室预约系统(book)

今日还可看 0 篇未注册用户每天查看4篇，注册用户每天8篇，开通VIP会员无限制查看。

免费注册

《高等数学》难点解析

高数教程 泰勒公式切线与法线切平面与法平面驻点·拐点·极值点·零点间断点渐进线瑕积分欧拉方程伯努利方程 Abel 收敛定理偏导数的几何意义偏导数的几何意义梯度数量场与向量场多元函数极值拉格朗日算子通量与散度环流量与旋度格林公式高斯公式斯托克斯公式三大公式比较傅里叶级数极坐标微元点法式方程变上限定积分 X型计算面积 Y型计算面积微分的意义渐近线间断点 y''+py'+qy=f(x)方程高斯黎曼傅里叶变换(复数) 拉普拉斯变换(复数)

《线性代数》难点解析

线代教程 近世代数对数学的整体思考线性的意义矩阵乘法(列视角) 矩阵乘法(行视角) 矩阵左乘矩阵右乘逆矩阵求解方程组阶梯形矩阵的求法方程组解的判定四阶行列式的计算线性变换的意义线性空间向量组的等价线性空间的几何意义基础解系的求法施密特正交化特征值与特征向量的意义矩阵相似的几何意义矩阵可对角化的理解秩的意义(向量版) 秩的意义(方程版) 二次型的意义

《概率论与数理统计》难点解析

概率教程 置信区间与上a分位数概率中的“取”与“放” 贝叶斯公式全概率公式泊松分布指数分布伽玛分布二维密度图的意义卷积的意义相关系数的意义 k阶矩是与矩母函数卡方分布的作用单正态区间估计理解假设检验理解切比雪夫不等式中心极限定理

上一篇：改进欧拉法

下一篇：线性多步法

本文对您是否有用？有用 (0) 无用 (0)