9.3显式龙格-库塔（Runge-Kutta, R-K）法

显式龙格-库塔（Runge-Kutta, R-K）法 详细讲解与推导

各位同学，今天我们用多年数值分析研究与教学的经验，从底层逻辑出发，完整讲解显式龙格-库塔法的引入、推导、计算实例与核心特性，彻底讲透这个常微分方程数值求解的核心方法。

---

一、前置基础：常微分方程初值问题与数值解法的核心思想

我们要解决的核心问题是一阶常微分方程初值问题（IVP）：

\[\begin{cases} y'(x) = f(x, y(x)) \\ y(x_0) = y_0 \end{cases} \]

其中\(f(x,y)\)在求解区域内满足Lipschitz条件，保证解的存在唯一性。

数值解法的核心是离散化：将求解区间\([a,b]\)做等距剖分，节点\(x_n = x_0 + nh\)（\(n=0,1,2,\dots,N\)），\(h\)为步长，我们的目标是求精确解\(y(x_n)\)的近似值\(y_n\)。

最基础的数值方法是欧拉法，公式为：

\[y_{n+1} = y_n + h f(x_n, y_n) \]

它用前向差商代替导数，本质是左矩形数值积分，仅具有一阶精度，局部截断误差为\(O(h^2)\)，精度太低，无法满足工程计算需求，因此我们需要更高精度的方法。

---

二、改进欧拉法：二阶R-K法的雏形

2.1 改进欧拉法的来源

对\(y'(x)=f(x,y)\)在区间\([x_n, x_{n+1}]\)上积分，得：

\[y(x_{n+1}) - y(x_n) = \int_{x_n}^{x_{n+1}} f(x,y(x))dx \]

欧拉法用左矩形公式近似积分，精度低；我们改用精度更高的梯形求积公式：

\[\int_{a}^{b} g(x)dx \approx \frac{b-a}{2}\left[ g(a) + g(b) \right] \]

代入得梯形法公式：

\[y_{n+1} = y_n + \frac{h}{2}\left[ f(x_n,y_n) + f(x_{n+1},y_{n+1}) \right] \]

但该式右侧含未知量\(y_{n+1}\)，是隐式公式，无法直接计算。因此我们用欧拉法先计算一个预测值\(\tilde{y}_{n+1}\)，再代入梯形公式做校正，得到预测-校正格式，即改进欧拉法：

\[\begin{cases} \tilde{y}_{n+1} = y_n + h f(x_n, y_n) \quad \text{预测步} \\ y_{n+1} = y_n + \frac{h}{2}\left[ f(x_n,y_n) + f(x_{n+1}, \tilde{y}_{n+1}) \right] \quad \text{校正步} \end{cases} \]

将预测步代入校正步，得到紧凑形式（即教材公式9.16）：

\[\boldsymbol{y_{n+1} = y_n + \frac{h}{2}\left[ f(x_n,y_n) + f(x_n + h, y_n + h f(x_n,y_n)) \right], \quad n=0,1,2,\dots} \]

对应的增量函数（教材公式9.17）：

\[\boldsymbol{\varphi(x_n,y_n,h) = \frac{1}{2}\left[ f(x_n,y_n) + f(x_n + h, y_n + h f(x_n,y_n)) \right]} \]

2.2 改进欧拉法的局部截断误差推导（核心）

定义：局部截断误差

假设\(x_n\)处的近似值是精确的，即\(y_n = y(x_n)\)，则数值方法计算的\(y_{n+1}\)与精确值\(y(x_{n+1})\)的差值，称为局部截断误差\(T_{n+1}\)，即：

\[T_{n+1} = y(x_{n+1}) - y_{n+1} \]

代入改进欧拉法公式，得（教材公式9.18修正版，原公式漏写系数\(\frac{1}{2}\)）：

\[\boldsymbol{T_{n+1} = y(x_{n+1}) - y(x_n) - \frac{h}{2}\left[ f(x_n,y(x_n)) + f(x_n + h, y(x_n) + h f(x_n,y(x_n))) \right]} \]

步骤1：精确解\(y(x_{n+1})\)的泰勒展开

\(x_{n+1}=x_n+h\)，对\(y(x_{n+1})\)在\(x_n\)处做泰勒展开到\(h^3\)项：

\[y(x_{n+1}) = y(x_n) + h y'(x_n) + \frac{h^2}{2!} y''(x_n) + \frac{h^3}{3!} y'''(x_n) + O(h^4) \]

利用\(y'(x)=f(x,y(x))\)，将\(y\)的各阶导数转化为\(f\)的偏导数：

1. 一阶导数：\(y'(x_n) = f(x_n,y(x_n)) = f_n\)
2. 二阶导数（复合函数全导数）：

\[y''(x_n) = \frac{d}{dx}f(x,y(x))\bigg|_{x_n} = f_x'(x_n,y_n) + f_y'(x_n,y_n) \cdot f_n \]

3. 三阶导数（再次求全导数，混合偏导连续时\(f_{xy}''=f_{yx}''\)）：

\[\begin{aligned} y'''(x_n) &= f_{xx}''(x_n,y_n) + 2f_n f_{xy}''(x_n,y_n) + f_n^2 f_{yy}''(x_n,y_n) \\ &+ f_y'(x_n,y_n)\left[ f_x'(x_n,y_n) + f_n f_y'(x_n,y_n) \right] \end{aligned} \]

步骤2：二元函数\(f(x_n+h, y_n+hf_n)\)的泰勒展开

二元函数\(f(x+\Delta x, y+\Delta y)\)在\((x_n,y_n)\)处的泰勒展开到\(h^2\)项（\(\Delta x=h, \Delta y=hf_n\)，均为\(O(h)\)）：

\[\begin{aligned} f(x_n+h, y_n+hf_n) &= f_n + \left( h \frac{\partial}{\partial x} + hf_n \frac{\partial}{\partial y} \right)f_n + \frac{1}{2!}\left( h \frac{\partial}{\partial x} + hf_n \frac{\partial}{\partial y} \right)^2 f_n + O(h^3) \\ &= \boldsymbol{f_n + h\left( f_x' + f_n f_y' \right) + \frac{h^2}{2}\left( f_{xx}'' + 2f_n f_{xy}'' + f_n^2 f_{yy}'' \right) + O(h^3)} \end{aligned} \]

步骤3：代入局部截断误差公式，合并同类项

将上述展开式代入\(T_{n+1}\)：

1. 先展开改进欧拉法的数值解：

\[\begin{aligned} y_{n+1} &= y_n + \frac{h}{2}\left[ f_n + f(x_n+h,y_n+hf_n) \right] \\ &= y_n + \frac{h}{2}\left[ f_n + f_n + h(f_x'+f_n f_y') + \frac{h^2}{2}(f_{xx}''+2f_n f_{xy}''+f_n^2 f_{yy}'') + O(h^3) \right] \\ &= y_n + h f_n + \frac{h^2}{2}(f_x' + f_n f_y') + \frac{h^3}{4}(f_{xx}''+2f_n f_{xy}''+f_n^2 f_{yy}'') + O(h^4) \end{aligned} \]

2. 与精确解的泰勒展开式相减，得到局部截断误差：

\[\begin{aligned} T_{n+1} &= \left[ y_n + h f_n + \frac{h^2}{2}y''(x_n) + \frac{h^3}{6}y'''(x_n) + O(h^4) \right] - y_{n+1} \\ &= \frac{h^3}{6}\left[ f_{xx}'' + 2f_n f_{xy}'' + f_n^2 f_{yy}'' + f_y'(f_x' + f_n f_y') \right] - \frac{h^3}{4}\left[ f_{xx}'' + 2f_n f_{xy}'' + f_n^2 f_{yy}'' \right] + O(h^4) \\ &= \boldsymbol{\frac{h^3}{3!} f_y'(x_n,y_n)\left[ f_x'(x_n,y_n) + f_n f_y'(x_n,y_n) \right] + O(h^4)} \end{aligned} \]

核心结论

改进欧拉法的局部截断误差主项为\(O(h^3)\)，因此它是二阶精度的数值方法，比欧拉法精度提升一个量级。

---

三、显式龙格-库塔法的一般形式与阶数条件

3.1 从改进欧拉法到一般显式R-K法

我们观察改进欧拉法，可以将其改写为更通用的形式：

\[\begin{cases} y_{n+1} = y_n + h\left( \frac{1}{2}K_1 + \frac{1}{2}K_2 \right) \\ K_1 = f(x_n, y_n) \\ K_2 = f(x_n + h, y_n + h K_1) \end{cases} \]

这里用2个函数值\(K_1,K_2\)的线性组合，实现了二阶精度。我们自然可以推广：用\(r\)个函数值\(K_1,K_2,\dots,K_r\)的线性组合，构造更高精度的方法，这就是r级显式龙格-库塔法。

3.2 r级显式R-K法的通用公式（教材公式9.19-9.20）

\[\boldsymbol{y_{n+1} = y_n + h \sum_{i=1}^r c_i K_i} \]

其中

\[\begin{cases} \boldsymbol{K_1 = f(x_n, y_n)} \\ \boldsymbol{K_i = f\left( x_n + \lambda_i h, y_n + h \sum_{j=1}^{i-1} \mu_{ij} K_j \right), \quad i=2,3,\dots,r} \end{cases} \]

• \(r\)：方法的级数，即每步需要计算\(f\)的次数；
• \(c_i, \lambda_i, \mu_{ij}\)：待定常数，通过匹配泰勒展开的高阶项确定，目标是让方法达到最高的精度阶数\(p\)。

3.3 核心原理：阶数条件

R-K法的核心优势是：无需计算\(f\)的高阶偏导数，仅通过多个函数值的线性组合，就能匹配泰勒展开的高阶项，实现高精度。

以二级R-K法为例，推导阶数条件：
二级显式R-K法的一般形式为：

\[\begin{cases} y_{n+1} = y_n + h(c_1 K_1 + c_2 K_2) \\ K_1 = f(x_n,y_n) \\ K_2 = f(x_n + \lambda h, y_n + h \mu K_1) \end{cases} \]

将\(K_2\)泰勒展开并代入\(y_{n+1}\)，与精确解的泰勒展开对比，令\(h、h^2\)项的系数相等，得到二级二阶R-K法的阶数条件：

\[\begin{cases} c_1 + c_2 = 1 \\ c_2 \lambda = \frac{1}{2} \\ c_2 \mu = \frac{1}{2} \end{cases} \]

该方程组有无穷多组解，对应不同的二级二阶R-K方法：

1. 取\(c_1=c_2=\frac{1}{2}, \lambda=\mu=1\)：即改进欧拉法；
2. 取\(c_1=0, c_2=1, \lambda=\mu=\frac{1}{2}\)：即中点法；
3. 取\(c_1=\frac{1}{4}, c_2=\frac{3}{4}, \lambda=\mu=\frac{2}{3}\)：即Heun法。

3.4 经典四级四阶显式R-K法

工程与科学计算中最常用的是经典四级四阶R-K法，它的性价比最高：4级计算量可达到4阶精度，而5级及以上的R-K法，精度提升远慢于计算量的增长。

经典四级四阶R-K法公式：

\[\boldsymbol{ \begin{cases} y_{n+1} = y_n + \frac{h}{6}\left( K_1 + 2K_2 + 2K_3 + K_4 \right) \\ K_1 = f(x_n, y_n) \\ K_2 = f\left( x_n + \frac{h}{2}, y_n + \frac{h}{2}K_1 \right) \\ K_3 = f\left( x_n + \frac{h}{2}, y_n + \frac{h}{2}K_2 \right) \\ K_4 = f\left( x_n + h, y_n + h K_3 \right) \end{cases} } \]

其局部截断误差为\(O(h^5)\)，具有四阶精度，是目前应用最广泛的R-K方法。

---

四、数值计算实例（迭代3次）

例题

求解常微分方程初值问题：

\[\begin{cases} y' = x + y \\ y(0) = 1 \end{cases} \]

取步长\(h=0.1\)，分别用欧拉法、改进欧拉法、经典四阶R-K法迭代3次，计算\(x=0.1,0.2,0.3\)处的数值解，并与精确解对比。

精确解推导

该方程为一阶线性非齐次微分方程，积分因子\(\mu=e^{-x}\)，解得：

\[y(x) = 2e^x - x - 1 \]

---

1. 欧拉法（一阶精度）

公式：\(y_{n+1} = y_n + h(x_n + y_n)\)，\(x_0=0,y_0=1\)

• 第1次迭代（\(x_1=0.1\)）：\(y_1 = 1 + 0.1\times(0+1) = 1.1\)
• 第2次迭代（\(x_2=0.2\)）：\(y_2 = 1.1 + 0.1\times(0.1+1.1) = 1.22\)
• 第3次迭代（\(x_3=0.3\)）：\(y_3 = 1.22 + 0.1\times(0.2+1.22) = 1.362\)

---

2. 改进欧拉法（二阶精度）

公式：

\[\begin{cases} K_1 = x_n + y_n \\ K_2 = (x_n+h) + (y_n + h K_1) \\ y_{n+1} = y_n + \frac{h}{2}(K_1+K_2) \end{cases} \]

• 第1次迭代（\(x_1=0.1\)）：\(K_1=1,K_2=1.2\)，\(y_1=1 + 0.05\times2.2=1.11\)
• 第2次迭代（\(x_2=0.2\)）：\(K_1=1.21,K_2=1.431\)，\(y_2=1.11 + 0.05\times2.641=1.24205\)
• 第3次迭代（\(x_3=0.3\)）：\(K_1=1.44205,K_2=1.686255\)，\(y_3=1.24205 + 0.05\times3.128305=1.39846525\)

---

3. 经典四级四阶R-K法（四阶精度）

公式见3.4节，迭代结果：

• 第1次迭代（\(x_1=0.1\)）：\(y_1≈1.110341667\)
• 第2次迭代（\(x_2=0.2\)）：\(y_2≈1.242805142\)
• 第3次迭代（\(x_3=0.3\)）：\(y_3≈1.399716994\)

---

五、核心知识点归纳总结

表1 常用显式R-K方法核心参数对比

方法名称	级数r	精度阶数p	每步计算f的次数	局部截断误差	核心特点
欧拉法	1	1	1	\(O(h^2)\)	计算量最小，精度最低，仅用于教学演示
改进欧拉法	2	2	2	\(O(h^3)\)	计算量适中，二阶精度，适合简单工程计算
中点法	2	2	2	\(O(h^3)\)	二级二阶，稳定性优于改进欧拉法
经典四级四阶R-K法	4	4	4	\(O(h^5)\)	精度高，通用性强，工业界标准方法

---

表2 显式R-K法核心特性总结

项目	详细说明
核心思想	通过多个节点的函数值线性组合，匹配泰勒展开高阶项，在不求高阶偏导的前提下提升精度
显式特性	\(K_i\)仅依赖前序的\(K_1\sim K_{i-1}\)，可直接递推计算，无需解方程；隐式R-K法需迭代求解，稳定性更强，适合刚性方程
阶数与级数的关系	1级→1阶、2级→2阶、3级→3阶、4级→4阶；级数≥5时，精度阶数提升远慢于计算量增长，因此四级四阶是性价比最优选择
核心优点	1. 自启动，仅需前一个节点的\(y_n\)即可计算，无需历史数据；2. 编程实现简单，仅需计算函数值，无需求导；3. 步长调整灵活，可动态适配求解需求；4. 四阶方法精度足以满足绝大多数工程计算需求
核心缺点	1. 计算量大于线性多步法，四阶方法每步需4次函数计算；2. 显式格式对刚性方程稳定性差，需极小步长，计算效率低；3. 局部截断误差估计比线性多步法复杂
适用场景	非刚性一阶常微分方程（组）初值问题，广泛应用于物理模拟、工程计算、控制系统仿真、数值分析等领域

---

表3 例题计算结果与精度对比

\(x_n\)	精确值\(y(x_n)\)	欧拉法\(y_n\)	欧拉法绝对误差	改进欧拉法\(y_n\)	改进欧拉法绝对误差	经典四阶R-K法\(y_n\)	经典四阶R-K法绝对误差
0.0	1.0	1.0	0	1.0	0	1.0	0
0.1	1.110341836	1.1	\(1.03\times10^{-2}\)	1.11	\(3.42\times10^{-4}\)	1.110341667	\(1.7\times10^{-7}\)
0.2	1.242805516	1.22	\(2.28\times10^{-2}\)	1.24205	\(7.56\times10^{-4}\)	1.242805142	\(3.7\times10^{-7}\)
0.3	1.399717616	1.362	\(3.77\times10^{-2}\)	1.39846525	\(1.25\times10^{-3}\)	1.399716994	\(6.2\times10^{-7}\)

从对比可以清晰看到：精度从低到高依次为欧拉法 < 改进欧拉法 < 经典四阶R-K法，四阶R-K法的误差比欧拉法小5个数量级，精度优势极其显著。

---

二阶显式龙格-库塔（R-K）方法 完整讲解与推导

一、二阶显式R-K方法的通用形式

对于二级（r=2）显式R-K方法，其核心是通过2个斜率值\(K_1,K_2\)的线性组合，构造数值迭代格式，通用公式为：

\[\boldsymbol{ \begin{cases} K_1 = f(x_n, y_n) \\ K_2 = f\left(x_n + \lambda_2 h, y_n + \mu_{21} h K_1\right) \\ y_{n+1} = y_n + h\left(c_1 K_1 + c_2 K_2\right) \end{cases} \quad n=0,1,2,\dots } \]

其中：

• \(h\)为求解步长，\(x_n = x_0 + nh\)为离散节点；
• \(c_1,c_2\)为斜率的权系数，\(\lambda_2\)为\(K_2\)的x方向步长系数，\(\mu_{21}\)为\(K_2\)的y方向增量系数，均为待定常数；
• 改进欧拉法是该格式的一个特例（对应\(c_1=c_2=\frac{1}{2},\lambda_2=1,\mu_{21}=1\)）。

---

二、核心推导：二阶精度的阶数条件

2.1 局部截断误差的定义

数值方法的精度由局部截断误差决定：假设前一步数值解等于精确解，即\(y_n = y(x_n)\)，则局部截断误差为：

\[T_{n+1} = y(x_{n+1}) - y_{n+1} \]

若\(T_{n+1}=O(h^{p+1})\)，则称方法具有\(p\)阶精度。我们的目标是确定参数，让方法达到二阶精度（即\(T_{n+1}=O(h^3)\)）。

2.2 泰勒展开准备

（1）精确解\(y(x_{n+1})\)的泰勒展开

将\(y(x_{n+1})=y(x_n+h)\)在\(x_n\)处泰勒展开到\(h^3\)项，结合\(y'(x)=f(x,y(x))\)，将导数转化为\(f\)的偏导数：

\[y(x_{n+1}) = y(x_n) + h y'(x_n) + \frac{h^2}{2!} y''(x_n) + \frac{h^3}{3!} y'''(x_n) + O(h^4) \]

其中：

• 一阶导数：\(y'(x_n) = f(x_n,y_n) = f_n\)
• 二阶导数（复合函数全导数）：\(y''(x_n) = f_x'(x_n,y_n) + f_y'(x_n,y_n) \cdot f_n\)
• 三阶导数：\(y'''(x_n) = f_{xx}'' + 2f_n f_{xy}'' + f_n^2 f_{yy}'' + f_y'(f_x' + f_n f_y')\)（用于后续精度上限分析）

（2）\(K_2\)的二元泰勒展开

\(K_2\)是二元函数\(f\)在\((x_n,y_n)\)处的展开，\(\Delta x=\lambda_2 h\)，\(\Delta y=\mu_{21} h f_n\)，展开到\(h^2\)项：

\[\begin{aligned} K_2 &= f(x_n + \lambda_2 h, y_n + \mu_{21} h f_n) \\ &= f_n + h\left( \lambda_2 f_x' + \mu_{21} f_n f_y' \right) + \frac{h^2}{2}\left( \lambda_2^2 f_{xx}'' + 2\lambda_2 \mu_{21} f_n f_{xy}'' + \mu_{21}^2 f_n^2 f_{yy}'' \right) + O(h^3) \end{aligned} \]

2.3 代入推导阶数条件

将\(K_1,K_2\)代入\(y_{n+1}\)的迭代公式，展开到\(h^2\)项：

\[\begin{aligned} y_{n+1} &= y_n + h\left( c_1 K_1 + c_2 K_2 \right) \\ &= y_n + h(c_1 + c_2)f_n + h^2 c_2\left( \lambda_2 f_x' + \mu_{21} f_n f_y' \right) + O(h^3) \end{aligned} \]

将其与精确解的展开式相减，得到局部截断误差：

\[\begin{aligned} T_{n+1} &= y(x_{n+1}) - y_{n+1} \\ &= h\cdot\left[1 - (c_1 + c_2)\right]f_n + h^2\cdot\left[ \left(\frac{1}{2} - c_2\lambda_2\right)f_x' + \left(\frac{1}{2} - c_2\mu_{21}\right)f_n f_y' \right] + O(h^3) \end{aligned} \]

要让方法达到二阶精度，必须让\(h\)的一次项、\(h\)的二次项系数全部为0，得到二阶显式R-K方法的阶数条件方程组：

\[\boldsymbol{ \begin{cases} 1 - c_1 - c_2 = 0 \\ \frac{1}{2} - c_2 \lambda_2 = 0 \\ \frac{1}{2} - c_2 \mu_{21} = 0 \end{cases} } \]

简化为：

\[c_1 + c_2 = 1, \quad c_2\lambda_2 = \frac{1}{2}, \quad c_2\mu_{21} = \frac{1}{2} \]

---

三、方程组的解与经典二阶R-K格式

该方程组有3个方程、4个未知量，是欠定方程组，有无穷多组解。我们取自由参数\(c_2 = a \ (a\neq0)\)，可解得通解：

\[\boldsymbol{c_1 = 1-a, \quad \lambda_2 = \mu_{21} = \frac{1}{2a}} \]

不同的\(a\)对应不同的二阶R-K格式，以下是3种最经典的特例：

特例1：改进欧拉法（Heun法）

取\(a=\frac{1}{2}\)，即\(c_2=\frac{1}{2}\)，代入通解得：

\[c_1=\frac{1}{2}, \quad \lambda_2=1, \quad \mu_{21}=1 \]

对应格式：

\[\begin{cases} K_1 = f(x_n, y_n) \\ K_2 = f(x_n + h, y_n + h K_1) \\ y_{n+1} = y_n + \frac{h}{2}(K_1 + K_2) \end{cases} \]

这是最常用的二阶R-K方法，对应数值积分的梯形公式，通用性强。

特例2：中点法（中点公式）

取\(a=1\)，即\(c_2=1\)，代入通解得：

\[c_1=0, \quad \lambda_2=\frac{1}{2}, \quad \mu_{21}=\frac{1}{2} \]

对应格式：

\[\begin{cases} K_1 = f(x_n, y_n) \\ K_2 = f\left(x_n + \frac{h}{2}, y_n + \frac{h}{2} K_1\right) \\ y_{n+1} = y_n + h K_2 \end{cases} \]

该方法用区间中点的斜率作为整体平均斜率，对应数值积分的中矩形公式，数值稳定性优于改进欧拉法。

特例3：Ralston法（最优精度二阶R-K法）

取\(a=\frac{2}{3}\)，即\(c_2=\frac{2}{3}\)，代入通解得：

\[c_1=\frac{1}{3}, \quad \lambda_2=\frac{3}{4}, \quad \mu_{21}=\frac{3}{4} \]

对应格式：

\[\begin{cases} K_1 = f(x_n, y_n) \\ K_2 = f\left(x_n + \frac{3}{4}h, y_n + \frac{3}{4}h K_1\right) \\ y_{n+1} = y_n + h\left( \frac{1}{3}K_1 + \frac{2}{3}K_2 \right) \end{cases} \]

该方法的局部截断误差主项系数最小，是所有二阶R-K方法中精度最优的格式。

---

四、关键结论：二级显式R-K方法的精度上限

二级显式R-K方法最高只能达到二阶精度，无法达到三阶精度，证明如下：

若要达到三阶精度，需要让局部截断误差的\(h^3\)项系数也为0。将\(y(x_{n+1})\)和\(K_2\)展开到\(h^3\)项，代入\(T_{n+1}\)后会发现：
\(h^3\)项中存在一项\(\frac{1}{6}f_y'(x_n,y_n)\left[f_x'(x_n,y_n) + f_n f_y'(x_n,y_n)\right]\)，该项的系数为固定常数\(\frac{1}{6}\)，无论如何选择参数都无法抵消。

因此，二级显式R-K方法的局部截断误差主项永远是\(O(h^3)\)，精度阶数只能为2，无法达到三阶。

---

五、核心知识点归纳总结

表1 二阶显式R-K方法核心信息汇总

项目	详细说明
通用形式	二级2斜率显式格式，每步仅需2次函数求值，自启动，无需历史数据
二阶精度条件	\(c_1+c_2=1\)，\(c_2\lambda_2=1/2\)，\(c_2\mu_{21}=1/2\)，有无穷多组解
经典格式	改进欧拉法、中点法、Ralston最优精度法
精度上限	最高二阶精度，局部截断误差主项为\(O(h^3)\)，无法达到三阶
核心优势	计算量远小于高阶R-K方法，精度比一阶欧拉法高一个量级，编程实现简单，步长调整灵活
适用场景	非刚性一阶常微分方程初值问题，低精度需求的工程计算、数值模拟入门、教学演示等

表2 经典二阶R-K格式参数与特性对比

方法名称	\(c_1\)	\(c_2\)	\(\lambda_2\)	\(\mu_{21}\)	核心特性
改进欧拉法	\(1/2\)	\(1/2\)	\(1\)	\(1\)	通用性最强，对应梯形求积公式，易理解易实现
中点法	\(0\)	\(1\)	\(1/2\)	\(1/2\)	数值稳定性好，对应中矩形求积公式，适合对稳定性有要求的场景
Ralston法	\(1/3\)	\(2/3\)	\(3/4\)	\(3/4\)	截断误差最小，二阶方法中精度最优，适合对精度有要求的轻量化计算

---

三阶与四阶显式龙格-库塔（R-K）方法 完整讲解与推导

一、前置核心逻辑回顾

显式R-K方法的核心是：用r个离散点的斜率值\(K_1,K_2,\dots,K_r\)的线性组合，匹配常微分方程精确解的泰勒展开高阶项，在无需计算函数高阶偏导数的前提下，实现高精度数值求解。

• 1级R-K（欧拉法）：最高1阶精度，局部截断误差\(O(h^2)\)
• 2级R-K：最高2阶精度，局部截断误差\(O(h^3)\)
• 要实现三阶精度，必须取级数\(r=3\)；要实现四阶精度，必须取级数\(r=4\)。

---

二、三阶显式R-K方法

2.1 三级显式R-K方法的通用形式

当\(r=3\)时，显式R-K方法的通用迭代公式为：

\[\boldsymbol{ \begin{cases} K_1 = f(x_n, y_n) \\ K_2 = f\left(x_n + \lambda_2 h, y_n + \mu_{21} h K_1\right) \\ K_3 = f\left(x_n + \lambda_3 h, y_n + \mu_{31} h K_1 + \mu_{32} h K_2\right) \\ y_{n+1} = y_n + h\left(c_1 K_1 + c_2 K_2 + c_3 K_3\right) \end{cases} \quad n=0,1,2,\dots } \]

其中：

• \(K_1,K_2,K_3\)为3个节点处的斜率值，每步迭代需计算3次\(f(x,y)\)；
• 待定参数共8个：权系数\(c_1,c_2,c_3\)，x方向步长系数\(\lambda_2,\lambda_3\)，y方向增量系数\(\mu_{21},\mu_{31},\mu_{32}\)；
• 显式特性：\(K_i\)仅依赖前序的\(K_1\sim K_{i-1}\)，可直接递推计算，无需解方程。

2.2 三阶精度的阶数条件推导

局部截断误差定义

假设\(y_n = y(x_n)\)（前一步数值解精确），局部截断误差为：

\[T_{n+1} = y(x_{n+1}) - y_{n+1} = y(x_n+h) - y(x_n) - h\left(c_1 K_1 + c_2 K_2 + c_3 K_3\right) \]

要实现三阶精度，需满足\(T_{n+1} = O(h^4)\)，即泰勒展开中\(h、h^2、h^3\)项的系数全部为0。

泰勒展开与系数匹配

1. 精确解的泰勒展开（到\(h^3\)项）

\[y(x_n+h) = y(x_n) + h f_n + \frac{h^2}{2}\left(f_x' + f_n f_y'\right) + \frac{h^3}{6}\left(f_{xx}'' + 2f_n f_{xy}'' + f_n^2 f_{yy}'' + f_y'(f_x' + f_n f_y')\right) + O(h^4) \]

其中\(f_n = f(x_n,y_n)\)，\(f_x',f_y'\)为\(f\)的一阶偏导数，\(f_{xx}'',f_{xy}'',f_{yy}''\)为二阶偏导数。

2. \(K_2,K_3\)的二元泰勒展开

• \(K_2\)展开到\(h^2\)项：

\[K_2 = f_n + h\left(\lambda_2 f_x' + \mu_{21} f_n f_y'\right) + \frac{h^2}{2}\left(\lambda_2^2 f_{xx}'' + 2\lambda_2 \mu_{21} f_n f_{xy}'' + \mu_{21}^2 f_n^2 f_{yy}''\right) + O(h^3) \]

• \(K_3\)展开到\(h^2\)项（需代入\(K_2\)的展开式）：

\[\begin{aligned} K_3 &= f_n + h\left(\lambda_3 f_x' + (\mu_{31}+\mu_{32}) f_n f_y'\right) + h^2\left[ \lambda_3 \mu_{32} f_y'(f_x' + f_n f_y') + \frac{1}{2}\lambda_3^2 f_{xx}'' + \lambda_3 (\mu_{31}+\mu_{32}) f_n f_{xy}'' + \frac{1}{2}(\mu_{31}+\mu_{32})^2 f_n^2 f_{yy}'' \right] + O(h^3) \end{aligned} \]

3. 代入匹配，得到阶数方程组
   将\(K_1,K_2,K_3\)代入\(y_{n+1}\)，与精确解展开式逐项对比，令\(h、h^2、h^3\)项系数为0，得到三阶显式R-K方法的阶数条件：

\[\boldsymbol{ \begin{cases} c_1 + c_2 + c_3 = 1 \quad \text{匹配h项系数} \\ \lambda_2 = \mu_{21} \\ \lambda_3 = \mu_{31} + \mu_{32} \quad \text{节点一致性条件} \\ c_2 \lambda_2 + c_3 \lambda_3 = \frac{1}{2} \quad \text{匹配h²项系数} \\ c_2 \lambda_2^2 + c_3 \lambda_3^2 = \frac{1}{3} \quad \text{匹配h³项二阶偏导系数} \\ c_3 \lambda_2 \mu_{32} = \frac{1}{6} \quad \text{匹配h³项混合偏导系数} \end{cases} } \]

该方程组为8个未知量、6个方程的欠定非线性方程组，有无穷多组解，所有满足该条件的格式均为三阶显式R-K方法。

2.3 经典三阶R-K格式：库塔三阶方法

最常用的三阶R-K格式为库塔三阶方法，其参数取值为：

\[c_1=\frac{1}{6}, c_2=\frac{4}{6}, c_3=\frac{1}{6}, \lambda_2=\frac{1}{2}, \lambda_3=1, \mu_{21}=\frac{1}{2}, \mu_{31}=-1, \mu_{32}=2 \]

对应迭代公式：

\[\boldsymbol{ \begin{cases} K_1 = f(x_n, y_n) \\ K_2 = f\left(x_n + \frac{h}{2}, y_n + \frac{h}{2} K_1\right) \\ K_3 = f\left(x_n + h, y_n - h K_1 + 2h K_2\right) \\ y_{n+1} = y_n + \frac{h}{6}\left(K_1 + 4K_2 + K_3\right) \end{cases} } \]

核心特性

1. 该格式对应数值积分中的辛普森（Simpson）公式，辛普森公式的代数精度为3，与三阶R-K方法的精度匹配；
2. 局部截断误差为\(O(h^4)\)，每步仅需3次函数求值，计算量适中，精度比二阶R-K方法提升一个量级；
3. 自启动，步长调整灵活，适合中等精度需求的工程计算。

---

三、四阶显式R-K方法

3.1 核心背景

三级R-K方法最高只能达到三阶精度，要实现四阶精度，必须取级数\(r=4\)。四阶显式R-K方法是工程与科学计算中应用最广泛的常微分方程数值解法，核心原因是：

• 四级四阶R-K方法每步4次函数求值，可达到四阶精度，局部截断误差\(O(h^5)\)；
• 当级数\(r\geq5\)时，精度阶数的提升远慢于计算量的增长（如5级R-K最高仅能达到4阶精度），因此四级四阶R-K方法是性价比最优的选择。

3.2 经典四级四阶R-K方法

四阶R-K方法的阶数条件更为复杂（11个方程、13个未知量），有无穷多组解，其中最经典、最常用的格式为经典龙格-库塔方法，迭代公式为：

\[\boldsymbol{ \begin{cases} K_1 = f(x_n, y_n) \\ K_2 = f\left(x_n + \frac{h}{2}, y_n + \frac{h}{2} K_1\right) \\ K_3 = f\left(x_n + \frac{h}{2}, y_n + \frac{h}{2} K_2\right) \\ K_4 = f\left(x_n + h, y_n + h K_3\right) \\ y_{n+1} = y_n + \frac{h}{6}\left(K_1 + 2K_2 + 2K_3 + K_4\right) \end{cases} } \]

每个斜率的物理意义

1. \(K_1\)：区间起点\(x_n\)处的斜率（欧拉法的斜率）；
2. \(K_2\)：用\(K_1\)预测区间中点\(x_n+h/2\)处的斜率；
3. \(K_3\)：用\(K_2\)修正后，重新计算区间中点\(x_n+h/2\)处的斜率，进一步提升精度；
4. \(K_4\)：用\(K_3\)预测区间终点\(x_n+h\)处的斜率；
5. 最终迭代值：对4个斜率做加权平均，中点斜率权重更高（2倍），起点和终点权重为1，符合数值积分的精度优化逻辑。

精度特性

经典四阶R-K方法的局部截断误差为\(O(h^5)\)，具有四阶精度，在步长\(h\)足够小时，数值解的误差会随步长缩小以\(h^4\)的速度下降，精度远高于二阶、三阶R-K方法。

---

四、数值计算实例（迭代3次）

我们沿用之前的常微分方程初值问题，对比不同阶数R-K方法的精度：

\[\begin{cases} y' = x + y \\ y(0) = 1 \end{cases} \]

精确解：\(y(x) = 2e^x - x - 1\)，取步长\(h=0.1\)，迭代3次，计算\(x=0.1,0.2,0.3\)处的数值解。

1. 库塔三阶方法计算结果

\(x_n\)	三阶R-K数值解	精确值	绝对误差
0.0	1.0	1.0	0
0.1	1.110333333	1.110341836	\(8.5\times10^{-6}\)
0.2	1.242792222	1.242805516	\(1.33\times10^{-5}\)
0.3	1.399698778	1.399717616	\(1.88\times10^{-5}\)

2. 经典四阶R-K方法计算结果

\(x_n\)	四阶R-K数值解	精确值	绝对误差
0.0	1.0	1.0	0
0.1	1.110341667	1.110341836	\(1.7\times10^{-7}\)
0.2	1.242805142	1.242805516	\(3.7\times10^{-7}\)
0.3	1.399716994	1.399717616	\(6.2\times10^{-7}\)

精度对比结论

四阶R-K方法的误差比三阶R-K方法小2个数量级，比二阶R-K方法小4个数量级，精度优势极其显著。

---

五、核心知识点归纳总结

表1 各阶显式R-K方法核心特性对比

方法类型	级数r	精度阶数p	每步函数求值次数	局部截断误差	核心特点
欧拉法	1	1	1	\(O(h^2)\)	计算量最小，精度最低，仅用于教学演示
二阶R-K（改进欧拉法）	2	2	2	\(O(h^3)\)	计算量适中，低精度需求场景适用
三阶R-K（库塔三阶）	3	3	3	\(O(h^4)\)	中等精度，计算量与精度平衡较好
四阶R-K（经典格式）	4	4	4	\(O(h^5)\)	精度高，通用性强，工业界标准方法，性价比最优

表2 三阶与四阶R-K方法核心参数与适用场景

项目	三阶显式R-K方法	四阶显式R-K方法
阶数条件	6个方程，8个未知量，无穷多解	11个方程，13个未知量，无穷多解
经典格式权重	\(K_1:1, K_2:4, K_3:1\)，总权重6	\(K_1:1, K_2:2, K_3:2, K_4:1\)，总权重6
计算效率	每步3次函数求值，适合大规模轻量化计算	每步4次函数求值，适合高精度工程计算
数值稳定性	稳定性优于二阶R-K，弱于四阶R-K	稳定性好，适用范围广
典型适用场景	中等精度的物理模拟、控制系统仿真、教学演示	高精度科学计算、工程仿真、有限元分析、航天航空计算

核心结论

1. 显式R-K方法的精度阶数与级数直接相关：\(r\)级显式R-K方法最高可达到\(r\)阶精度（\(r\leq4\)）；
2. 三阶、四阶R-K方法均为自启动方法，步长调整灵活，编程实现简单，无需历史数据，远优于线性多步法；
3. 经典四阶R-K方法是目前应用最广泛的格式，在精度、计算量、稳定性之间达到了完美平衡，是常微分方程初值问题数值求解的首选方法。

---

四阶龙格-库塔法求解常微分方程初值问题 完整解析

一、例题完整还原

原初值问题

求解一阶常微分方程初值问题：

\[\begin{cases} y' = f(x,y) = y - \frac{2x}{y} \\ y(0) = 1 \end{cases} \]

求解要求：取步长\(h=0.2\)，在区间\(x\in[0,1]\)上，用经典四阶龙格-库塔（R-K）方法求解数值解，并与精确解对比。

精确解推导

该方程为伯努利方程，通过变量代换可求得解析解：
令\(u=y^2\)，则\(u'=2yy'\)，原方程转化为一阶线性非齐次微分方程：

\[u' - 2u = -4x \]

结合初值\(u(0)=y^2(0)=1\)，解得\(u=1+2x\)，因此原方程的精确解为：

\[\boldsymbol{y(x) = \sqrt{1+2x}} \]

---

二、适配的经典四阶龙格-库塔公式

经典四阶R-K方法的通用格式，代入本题的\(f(x,y)=y-\frac{2x}{y}\)、步长\(h=0.2\)，得到迭代公式：

\[\boldsymbol{ \begin{cases} K_1 = y_n - \frac{2x_n}{y_n} \\ K_2 = \left(y_n + \frac{h}{2}K_1\right) - \frac{2x_n + h}{y_n + \frac{h}{2}K_1} \\ K_3 = \left(y_n + \frac{h}{2}K_2\right) - \frac{2x_n + h}{y_n + \frac{h}{2}K_2} \\ K_4 = \left(y_n + hK_3\right) - \frac{2(x_n + h)}{y_n + hK_3} \\ y_{n+1} = y_n + \frac{h}{6}\left(K_1 + 2K_2 + 2K_3 + K_4\right) \end{cases} } \]

其中\(h=0.2\)，\(\frac{h}{2}=0.1\)，\(\frac{h}{6}=\frac{1}{30}\approx0.0333333\)。

---

三、完整分步迭代计算过程

初始条件：\(n=0\)，\(x_0=0\)，\(y_0=1\)，共迭代5步，计算\(x=0.2,0.4,0.6,0.8,1.0\)处的数值解。

第1步：计算\(x_1=0.2\)处的\(y_1\)

1. \(K_1 = 1 - \frac{0}{1} = 1.0\)
2. \(y_0 + 0.1K_1 = 1.1\)，\(K_2 = 1.1 - \frac{0.2}{1.1} \approx 0.91818182\)
3. \(y_0 + 0.1K_2 = 1.09181818\)，\(K_3 = 1.09181818 - \frac{0.2}{1.09181818} \approx 0.90863816\)
4. \(y_0 + 0.2K_3 = 1.18172763\)，\(K_4 = 1.18172763 - \frac{0.4}{1.18172763} \approx 0.84324016\)
5. \(y_1 = 1 + \frac{0.2}{6}\times(1 + 2\times0.91818182 + 2\times0.90863816 + 0.84324016) \approx \boldsymbol{1.1832}\)

第2步：计算\(x_2=0.4\)处的\(y_2\)

1. \(K_1 = 1.1832 - \frac{0.4}{1.1832} \approx 0.84513374\)
2. \(y_1 + 0.1K_1 = 1.26771337\)，\(K_2 = 1.26771337 - \frac{0.6}{1.26771337} \approx 0.79442028\)
3. \(y_1 + 0.1K_2 = 1.26264203\)，\(K_3 = 1.26264203 - \frac{0.6}{1.26264203} \approx 0.78744796\)
4. \(y_1 + 0.2K_3 = 1.34068959\)，\(K_4 = 1.34068959 - \frac{0.8}{1.34068959} \approx 0.74398177\)
5. \(y_2 = 1.1832 + \frac{0.2}{6}\times(0.84513374 + 2\times0.79442028 + 2\times0.78744796 + 0.74398177) \approx \boldsymbol{1.3417}\)

第3步：计算\(x_3=0.6\)处的\(y_3\)

1. \(K_1 = 1.3417 - \frac{0.8}{1.3417} \approx 0.74544152\)
2. \(y_2 + 0.1K_1 = 1.41624415\)，\(K_2 = 1.41624415 - \frac{1.0}{1.41624415} \approx 0.71015113\)
3. \(y_2 + 0.1K_2 = 1.41271511\)，\(K_3 = 1.41271511 - \frac{1.0}{1.41271511} \approx 0.70485831\)
4. \(y_2 + 0.2K_3 = 1.48267166\)，\(K_4 = 1.48267166 - \frac{1.2}{1.48267166} \approx 0.67332186\)
5. \(y_3 = 1.3417 + \frac{0.2}{6}\times(0.74544152 + 2\times0.71015113 + 2\times0.70485831 + 0.67332186) \approx \boldsymbol{1.4833}\)

第4步：计算\(x_4=0.8\)处的\(y_4\)

1. \(K_1 = 1.4833 - \frac{1.2}{1.4833} \approx 0.67429306\)
2. \(y_3 + 0.1K_1 = 1.55072931\)，\(K_2 = 1.55072931 - \frac{1.4}{1.55072931} \approx 0.64792828\)
3. \(y_3 + 0.1K_2 = 1.54809283\)，\(K_3 = 1.54809283 - \frac{1.4}{1.54809283} \approx 0.64375433\)
4. \(y_3 + 0.2K_3 = 1.61205087\)，\(K_4 = 1.61205087 - \frac{1.6}{1.61205087} \approx 0.61952637\)
5. \(y_4 = 1.4833 + \frac{0.2}{6}\times(0.67429306 + 2\times0.64792828 + 2\times0.64375433 + 0.61952637) \approx \boldsymbol{1.6125}\)

第5步：计算\(x_5=1.0\)处的\(y_5\)

1. \(K_1 = 1.6125 - \frac{1.6}{1.6125} \approx 0.62025194\)
2. \(y_4 + 0.1K_1 = 1.67452519\)，\(K_2 = 1.67452519 - \frac{1.8}{1.67452519} \approx 0.59959364\)
3. \(y_4 + 0.1K_2 = 1.67245936\)，\(K_3 = 1.67245936 - \frac{1.8}{1.67245936} \approx 0.59620006\)
4. \(y_4 + 0.2K_3 = 1.73174001\)，\(K_4 = 1.73174001 - \frac{2.0}{1.73174001} \approx 0.57683223\)
5. \(y_5 = 1.6125 + \frac{0.2}{6}\times(0.62025194 + 2\times0.59959364 + 2\times0.59620006 + 0.57683223) \approx \boldsymbol{1.7321}\)

---

四、计算结果汇总与精度分析

数值解与精确解对比表

| \(x_n\) | 四阶R-K数值解\(y_n\) | 精确解\(y(x_n)=\sqrt{1+2x_n}\) | 绝对误差\(|y_n - y(x_n)|\) |
|-------|---------------------|--------------------------------|---------------------------|
| 0.0 | 1.0000 | 1.0000 | 0 |
| 0.2 | 1.1832 | 1.1832 | \(1.3\times10^{-5}\) |
| 0.4 | 1.3417 | 1.3416 | \(1.1\times10^{-4}\) |
| 0.6 | 1.4833 | 1.4832 | \(9.0\times10^{-5}\) |
| 0.8 | 1.6125 | 1.6125 | \(4.8\times10^{-5}\) |
| 1.0 | 1.7321 | 1.7321 | \(7.2\times10^{-5}\) |

精度分析

1. 四阶R-K方法的数值解与精确解高度吻合，最大绝对误差仅为\(1.1\times10^{-4}\)，在步长\(h=0.2\)的情况下，精度远高于二阶改进欧拉法。
2. 四阶R-K方法每步需4次函数求值，改进欧拉法每步仅需2次；但本题中步长放大一倍，总计算步数减少一半，整体计算量与改进欧拉法几乎持平，却获得了数量级更高的精度，充分体现了算法选择的重要性。

---

五、例题核心结论与知识点拓展

1. 四阶R-K方法的核心优势

   • 自启动特性：仅需前一个节点的数值解即可迭代，无需历史数据，步长调整灵活；
   • 精度与计算量平衡：四级四阶格式是常微分方程数值求解中性价比最高的格式，广泛应用于工程与科学计算；
   • 无需计算函数的高阶偏导数，仅通过函数值的线性组合即可实现四阶精度，编程实现简单。

2. 算法适用边界
   龙格-库塔方法的推导基于泰勒展开，要求待求解的函数具有良好的光滑性；若解的光滑性较差（如存在间断、尖点），四阶R-K方法的精度可能反而不如二阶改进欧拉法，实际计算中需根据问题特性选择合适的算法。

3. 工程应用价值
   本题的微分方程是典型的可分离变量方程，用于验证数值方法的精度；四阶R-K方法可直接推广到一阶微分方程组、高阶微分方程的求解，是控制系统仿真、物理过程模拟、有限元分析等领域的基础算法。

---

变步长的龙格-库塔（R-K）方法 完整讲解与推导

一、方法引入：固定步长R-K法的核心痛点

显式龙格-库塔法是常微分方程初值问题的主流求解方法，但固定步长的格式存在天然的局限性：

1. 精度与计算量的矛盾：单步来看，步长越小，局部截断误差越小；但步长缩小会导致固定求解区间内的总步数激增，不仅大幅提升计算量，还会引发计算机舍入误差的严重积累，最终反而降低数值解的可靠性。
2. 无法适配解的局部特性：微分方程的解在不同区间的变化特性差异极大——解变化平缓的区域，用大步长即可满足精度；解变化剧烈（如瞬态、脉冲、尖点）的区域，必须用小步长才能控制误差。固定步长只能按最苛刻的区域选择极小步长，造成大量不必要的计算开销。

因此，变步长R-K法的核心目标是：实现步长的自适应调整，在满足预设精度的前提下，最大化计算效率。要实现这一目标，必须解决两个核心问题：

1. 如何量化、检验当前步长下数值解的精度（误差估计）；
2. 如何根据精度检验结果，动态调整步长。

---

二、核心原理：基于步长折半的误差事后估计

我们以工程最常用的经典四阶显式R-K法为基础，推导变步长的核心逻辑。
四阶R-K法的核心特性是：单步局部截断误差为\(O(h^5)\)，即步长为\(h\)时，局部截断误差与步长的5次方成正比。

2.1 两次计算的误差表达式

从节点\(x_n\)出发，对同一目标节点\(x_{n+1}=x_n+h\)，做两次不同步长的计算：

1. 大步长单步计算：以\(h\)为步长，一步计算到\(x_{n+1}\)，得到数值解\(y_{n+1}^{(h)}\)。
   由于四阶R-K法的局部截断误差主项为\(ch^5\)（\(c\)是与\(f(x,y)\)的五阶偏导数相关的常数，在\(x_n\)的小邻域内近似为定值），因此有：

   \[\boldsymbol{y(x_{n+1}) - y_{n+1}^{(h)} \approx c h^5} \tag{1} \]

   其中\(y(x_{n+1})\)是\(x_{n+1}\)处的精确解。

2. 小步长两步计算：将步长折半为\(\frac{h}{2}\)，分两步走到\(x_{n+1}\)：第一步从\(x_n\)到\(x_{n+1/2}=x_n+\frac{h}{2}\)，第二步从\(x_{n+1/2}\)到\(x_{n+1}\)，最终得到数值解\(y_{n+1}^{(h/2)}\)。
   每一步\(\frac{h}{2}\)的截断误差主项为\(c\left(\frac{h}{2}\right)^5\)，两步的总截断误差为两步误差之和，因此有：

   \[\boldsymbol{y(x_{n+1}) - y_{n+1}^{(h/2)} \approx 2c \cdot \left(\frac{h}{2}\right)^5 = \frac{c h^5}{16}} \tag{2} \]

2.2 误差比与事后估计式推导

对比式(1)和式(2)，可得到核心结论：步长折半后，局部截断误差约减小到原来的\(\frac{1}{16}\)，即：

\[\frac{y(x_{n+1}) - y_{n+1}^{(h/2)}}{y(x_{n+1}) - y_{n+1}^{(h)}} \approx \frac{1}{16} \]

将式(1)减去式(2)，消去未知的精确解\(y(x_{n+1})\)和常数\(c\)，得到：

\[y_{n+1}^{(h/2)} - y_{n+1}^{(h)} \approx c h^5 - \frac{c h^5}{16} = \frac{15}{16} c h^5 \]

整理得：

\[c h^5 \approx \frac{16}{15} \left( y_{n+1}^{(h/2)} - y_{n+1}^{(h)} \right) \]

将上式代入式(2)，最终得到误差事后估计式：

\[\boldsymbol{y(x_{n+1}) - y_{n+1}^{(h/2)} \approx \frac{1}{15} \left| y_{n+1}^{(h/2)} - y_{n+1}^{(h)} \right|} \]

该式的核心价值

无需知道方程的精确解，也无需计算\(f(x,y)\)的高阶偏导数，仅通过步长折半前后两次数值解的偏差，即可估计出当前步长下数值解的真实误差，实现了精度的可量化、可检验。

---

三、变步长R-K法的具体执行流程

我们定义两次计算结果的偏差为：

\[\boldsymbol{\Delta = \left| y_{n+1}^{(h/2)} - y_{n+1}^{(h)} \right|} \]

结合误差估计式，\(\frac{\Delta}{15}\)即为步长\(\frac{h}{2}\)下数值解的真实误差近似值。给定预设的精度容许值\(\varepsilon\)，按以下两种情况动态调整步长：

情况1：误差超标，步长折半

若\(\boldsymbol{\Delta > \varepsilon}\)，说明当前步长\(h\)过大，截断误差超出了精度要求。

• 操作：将步长折半为\(\frac{h}{2}\)，重复上述“大步长+小步长”的计算与检验流程，直到\(\Delta < \varepsilon\)为止；
• 结果选取：取最终满足精度的折半步长计算结果\(y_{n+1}^{(h/2)}\)作为\(x_{n+1}\)处的数值解。

情况2：精度冗余，步长加倍

若\(\boldsymbol{\Delta < \varepsilon}\)，说明当前步长\(h\)过小，计算存在冗余，效率偏低。

• 操作：将步长加倍为\(2h\)，重复计算与检验流程，直到\(\Delta > \varepsilon\)为止；
• 结果选取：将最后一次超标的步长折半一次，得到既满足精度、又最大化步长的最优值，以此计算数值解。

补充说明

这种通过步长加倍/折半实现自适应调整的方法，是变步长R-K法最基础、最易实现的格式。表面上看，单步的计算量有所增加（四阶R-K法步长\(h\)单步需4次函数求值，步长\(\frac{h}{2}\)两步需8次函数求值，单步检验共需12次函数求值），但由于在解平缓的区域可以用极大的步长，总步数大幅减少，整体计算效率通常远高于固定步长格式。

---

四、拓展：更高效的变步长R-K法——嵌入型R-K格式

步长折半法虽然逻辑简单，但单步计算开销大，工程上更常用的是嵌入型（嵌入式）龙格-库塔法，也叫R-K-Fehlberg（RKF）法。

核心原理

嵌入型R-K法的核心创新是：构造一对\(p\)阶和\(p+1\)阶的R-K公式，共享全部斜率\(K_i\)，一次计算即可同时得到\(p\)阶和\(p+1\)阶两个数值解，用两个解的差值直接估计局部截断误差，无需重复计算，计算效率大幅提升。

最经典的是RKF45格式（四阶-五阶嵌入R-K法）：

• 用5阶R-K公式的结果作为参考真值，4阶R-K公式的结果作为数值解；
• 两个结果的差值即为局部截断误差的估计值；
• 一次计算仅需6次函数求值，即可完成误差估计与步长调整，效率远高于步长折半法。

MATLAB中最常用的常微分方程求解器ode45，核心就是RKF45格式，它是目前工业界自适应求解非刚性常微分方程的标准方法。

---

五、方法特性与适用场景总结

表1 变步长R-K法核心特性

项目	详细说明
核心优势	1. 自适应适配解的局部特性，平缓区用大步长降开销，剧烈区用小步长保精度； 2. 误差事后估计，实现精度可控，无需预知解的高阶特性； 3. 鲁棒性强，无需提前预判解的变化规律，通用型极强； 4. 整体计算效率远高于固定步长R-K法，尤其适合解变化差异大的问题。
局限性	1. 基础折半法单步计算开销高于固定步长； 2. 显式格式的本质限制，对强刚性方程稳定性差，无法通过变步长解决； 3. 步长频繁切换可能引入少量额外舍入误差，影响可忽略。
经典格式	1. 四阶R-K步长折半法（易实现，教学与入门常用）； 2. RKF45嵌入型R-K法（高效，工程与工业仿真主流）； 3. 其他嵌入格式：RKF23（二阶-三阶）、Dormand-Prince法（ode45的优化格式）。
适用场景	1. 解的变化率差异大的非刚性常微分方程初值问题； 2. 对数值解精度有明确要求，需严格控制误差的工程计算； 3. 无法提前预判解的变化特性，需要通用自适应求解的场景； 4. 物理过程仿真、控制系统瞬态分析、电路仿真等领域。

表2 固定步长与变步长R-K法对比

对比维度	固定步长R-K法	变步长R-K法
精度控制	无法实时控制，只能靠提前设置步长保证，易出现精度不足或冗余	实时误差估计，严格满足预设精度，无冗余
计算效率	解平缓区存在大量无效计算，整体效率低	自适应调整步长，总步数最少，整体效率高
实现难度	逻辑简单，极易编程	基础折半法实现简单，嵌入型格式实现稍复杂
通用性	仅适合解变化均匀的问题，通用性差	适配绝大多数非刚性常微分方程，通用性极强

---

六、核心结论

变步长龙格-库塔法解决了固定步长格式“精度与效率不可兼得”的核心痛点，通过误差事后估计和步长自适应调整，在保证精度的前提下最大化计算效率。

• 入门与教学场景，优先使用四阶R-K步长折半法，逻辑清晰，易理解、易实现；
• 工程与工业仿真场景，优先使用RKF45嵌入型R-K法，计算效率更高，是目前的行业标准。
  同时需注意，变步长仅能优化显式R-K法的精度与效率，无法解决其对强刚性方程的稳定性问题，刚性方程需采用隐式变步长R-K法求解。