一阶常微分方程组求解

2026年8月6日

1. 背景
在一些场景里，需要对带初值的一阶常微分方程组进行数值计算。由于龙格库塔法具有较高精度，因此这里主要记录这类问题使用龙格库塔法的处理方式。

这里约定时间变量为$x$，待求解的$m$个变量分别为$y_1,y_2,…,y_m$，因此一阶常微分方程组可以表示为如下形式
(1)
$$
\begin{cases}
y_1′ &= f_1(x,y_1,y_2,…,y_m)
\\
y_2′ &= f_2(x,y_1,y_2,…,y_m)
\\
…
\\
y_m’ &= f_m(x,y_1,y_2,…,y_m)
\end{cases}
$$
上面公式里$f_1,f_2,…,f_m$分别对应$y_1,y_2,…,y_m$一阶导函数。

往往我们知道一个启动初值，比如$x$在某个值时对应的所有$y_1,y_2,…,y_m$值，然后需要求$x$在其它位置处$y_1,y_2,…,y_m$的值。这就是本文要介绍的常微分初值问题。

2. 理论
2.1 基本过程
龙格库塔法在求解时，主要借助计算好的布彻表(Butcher Tableau)进行求解。布彻表主要包括系数$A、\mathbf{b}、\mathbf{c}$。对于$n$级的布彻表，这里$A$为$n\times n$的矩阵，$\mathbf{b}$与$\mathbf{c}$均为大小为$n$的向量。

实际处理中先根据$A、\mathbf{c}$计算出斜率$K$，具体公式如下
(2)
$$
\begin{cases}
&1级:\begin{cases}
K_{11} &=f_1(x+c_h,y_1+h\sum_{j=1}^nA_{1j}K_{j1},y_2+h\sum_{j=1}^nA_{1j}K_{j2},…,y_m+h\sum_{j=1}^nA_{1j}K_{jm})
\\
K_{12} &=f_2(x+c_h,y_1+h\sum_{j=1}^nA_{1j}K_{j1},y_2+h\sum_{j=1}^nA_{1j}K_{j2},…,y_m+h\sum_{j=1}^nA_{1j}K_{jm})
\\
…
\\
K_{1m} &=f_m(x+c_h,y_1+h\sum_{j=1}^nA_{1j}K_{j1},y_2+h\sum_{j=1}^nA_{1j}K_{j2},…,y_m+h\sum_{j=1}^nA_{1j}K_{jm})
\end{cases}
\\
\\
&2级:\begin{cases}
K_{21} &=f_1(x+c_h,y_1+h\sum_{j=1}^nA_{2j}K_{j1},y_2+h\sum_{j=1}^nA_{2j}K_{j2},…,y_m+h\sum_{j=1}^nA_{2j}K_{jm})
\\
K_{22} &=f_2(x+c_h,y_1+h\sum_{j=1}^nA_{2j}K_{j1},y_2+h\sum_{j=1}^nA_{2j}K_{j2},…,y_m+h\sum_{j=1}^nA_{2j}K_{jm})
\\
…
\\
K_{2m} &=f_m(x+c_h,y_1+h\sum_{j=1}^nA_{2j}K_{j1},y_2+h\sum_{j=1}^nA_{2j}K_{j2},…,y_m+h\sum_{j=1}^nA_{2j}K_{jm})
\end{cases}
\\
\\
…
\\
\\
&n级:\begin{cases}
K_{n1} &=f_1(x+c_h,y_1+h\sum_{j=1}^nA_{nj}K_{j1},y_2+h\sum_{j=1}^nA_{nj}K_{j2},…,y_m+h\sum_{j=1}^nA_{nj}K_{jm})
\\
K_{n2} &=f_2(x+c_h,y_1+h\sum_{j=1}^nA_{nj}K_{j1},y_2+h\sum_{j=1}^nA_{nj}K_{j2},…,y_m+h\sum_{j=1}^nA_{nj}K_{jm})
\\
…
\\
K_{nm} &=f_m(x+c_h,y_1+h\sum_{j=1}^nA_{nj}K_{j1},y_2+h\sum_{j=1}^nA_{nj}K_{j2},…,y_m+h\sum_{j=1}^nA_{nj}K_{jm})
\end{cases}
\end{cases}
$$

然后使用系数$\mathbf{b}$进行如下更新
(3)
$$
\begin{cases}
y_1 &=y_1 + h\sum_{i=1}^nb_iK_{i1}
\\
y_2 &=y_2 + h\sum_{i=1}^nb_iK_{i2}
\\
…
\\
y_m &=y_m + h\sum_{i=1}^nb_iK_{im}
\\
\\
x &= x + h
\end{cases}
$$
上面公式里$h$为$x$迭代步长。

因此主要过程为，(2)与(3)不停相互迭代直到满足条件。在迭代过程中，会根据迭代效果适当对$h$进行缩小或放大。

在实际处理过程中，(2)的计算方式的不同分为显式求解算法与隐式求解算法。

2.2 显式求解
显式求解主要适用于非刚性常微分方程组求解。所谓显式求解，指在第(2)步求解$K$的过程中，$K$可以直接根据等式右侧的表达式直接递推出来，这时求解就较快，使用这种方法就要求$A$矩阵对角线上三角必须全部为0。

考虑到动态自适应问题，实际处理时主要使用两套布彻表，在相同步长$h$下分别计算出各自的$y_1,y_2,…,y_m$，然后比较它们的误差，满足条件则进入下一个迭代，否则缩小$h$重新迭代。为了加快计算效率，一般这两套布彻表的系数$A、\mathbf{c}$一致，而仅$\mathbf{b}$不一致。

这里给出常见的Diori45算法系数
(4)
$$
\begin{cases}
A&=\begin{bmatrix}
0
\\
\\
\dfrac{1}{5} & 0
\\
\\
\dfrac{3}{40} &\dfrac{9}{40} &0
\\
\\
\dfrac{44}{45} & -\dfrac{56}{15} & \dfrac{32}{9} & 0
\\
\\
\dfrac{19372}{6561}&-\dfrac{25360}{2187}&\dfrac{64448}{6561}&-\dfrac{212}{729}&0
\\
\\
\dfrac{9017}{3168}&-\dfrac{355}{33}&\dfrac{46732}{5247}&\dfrac{49}{176}&-\dfrac{5103}{18656}&0
\\
\\
\dfrac{35}{384}&0&\dfrac{500}{1113}&\dfrac{125}{192}&-\dfrac{2187}{6784}&\dfrac{11}{84}&0
\end{bmatrix}
\\
\\
\mathbf{c}&=\begin{bmatrix}0 &\dfrac{1}{5}&\dfrac{3}{10}&\dfrac{4}{5}&\dfrac{8}{9}&1&1\end{bmatrix}^T
\\
\\
5阶:\mathbf{b}&=\begin{bmatrix}\dfrac{35}{384}&0&\dfrac{500}{1113}&\dfrac{125}{192}&-\dfrac{2187}{6784}&\dfrac{11}{84}&0\end{bmatrix}^T
\\
\\
4阶:\mathbf{b}&=\begin{bmatrix}\dfrac{5179}{57600}&0&\dfrac{7571}{16695}&\dfrac{393}{640}&-\dfrac{92097}{339200}&\dfrac{187}{2100}&\dfrac{1}{40}\end{bmatrix}^T
\end{cases}
$$

对于其它更高阶系数，可以参考Jim Verner^[1]的博客。

2.3 隐式求解
隐式求解，主要用于刚性常微分方程组的求解。所谓隐式求解，指在第(2)步求解$K$的过程中，$K$不能直接根据等式右侧的表达式直接递推出来，需要将(2)整体带入相应的求解器进行求解，这种求解效率就大打折扣。使用这种方法时$A$矩阵对角线上三角存在非0元素。

对于某个时刻，给定$x$以及$y_1,y_2,…,y_m$，这时可以得到$f_1,f_2,…,f_m$对$y_1,y_2,…,y_m$的导数，这里以一个$J$矩阵表示
(5)
$$
J=\begin{bmatrix}
\dfrac{\partial f_1}{\partial y_1} & \dfrac{\partial f_1}{\partial y_2} &…&\dfrac{\partial f_1}{\partial y_m}
\\
\\
\dfrac{\partial f_2}{\partial y_1} & \dfrac{\partial f_2}{\partial y_2} &…&\dfrac{\partial f_2}{\partial y_m}
\\
\\
…
\\
\\
\dfrac{\partial f_m}{\partial y_1} & \dfrac{\partial f_m}{\partial y_2} &…&\dfrac{\partial f_m}{\partial y_m}
\end{bmatrix}
$$
在求解(2)里的非线性方程组时，设置如下误差函数$g$
(6)
$$
\begin{cases}
&1级:\begin{cases}
g_{11} &=f_1(x+c_h,y_1+h\sum_{j=1}^nA_{1j}K_{j1},y_2+h\sum_{j=1}^nA_{1j}K_{j2},…,y_m+h\sum_{j=1}^nA_{1j}K_{jm})-K_{11}
\\
g_{12} &=f_2(x+c_h,y_1+h\sum_{j=1}^nA_{1j}K_{j1},y_2+h\sum_{j=1}^nA_{1j}K_{j2},…,y_m+h\sum_{j=1}^nA_{1j}K_{jm})-K_{12}
\\
…
\\
g_{1m} &=f_m(x+c_h,y_1+h\sum_{j=1}^nA_{1j}K_{j1},y_2+h\sum_{j=1}^nA_{1j}K_{j2},…,y_m+h\sum_{j=1}^nA_{1j}K_{jm})-K_{1m}
\end{cases}
\\
\\
&2级:\begin{cases}
g_{21} &=f_1(x+c_h,y_1+h\sum_{j=1}^nA_{2j}K_{j1},y_2+h\sum_{j=1}^nA_{2j}K_{j2},…,y_m+h\sum_{j=1}^nA_{2j}K_{jm})-K_{21}
\\
g_{22} &=f_2(x+c_h,y_1+h\sum_{j=1}^nA_{2j}K_{j1},y_2+h\sum_{j=1}^nA_{2j}K_{j2},…,y_m+h\sum_{j=1}^nA_{2j}K_{jm})-K_{22}
\\
…
\\
g_{2m} &=f_m(x+c_h,y_1+h\sum_{j=1}^nA_{2j}K_{j1},y_2+h\sum_{j=1}^nA_{2j}K_{j2},…,y_m+h\sum_{j=1}^nA_{2j}K_{jm})-K_{2m}
\end{cases}
\\
\\
…
\\
\\
&n级:\begin{cases}
g_{n1} &=f_1(x+c_h,y_1+h\sum_{j=1}^nA_{nj}K_{j1},y_2+h\sum_{j=1}^nA_{nj}K_{j2},…,y_m+h\sum_{j=1}^nA_{nj}K_{jm})-K_{n1}
\\
g_{n2} &=f_2(x+c_h,y_1+h\sum_{j=1}^nA_{nj}K_{j1},y_2+h\sum_{j=1}^nA_{nj}K_{j2},…,y_m+h\sum_{j=1}^nA_{nj}K_{jm})-K_{n2}
\\
…
\\
g_{nm} &=f_m(x+c_h,y_1+h\sum_{j=1}^nA_{nj}K_{j1},y_2+h\sum_{j=1}^nA_{nj}K_{j2},…,y_m+h\sum_{j=1}^nA_{nj}K_{jm})-K_{nm}
\end{cases}
\end{cases}
$$
观察上面公式就可以发现，在第$i$级，函数$f_1,f_2,…,f_m$它们接收的参数一致，因此，每一级都可以根据(5)计算得到一个单独的$J$，这时在第$i$里的误差函数$g$对参数求导有
(7)
$$
\dfrac{\partial g_{ij}}{\partial K_{pq}}=hA_{ip}J_{jq}-\begin{cases}
1&,i=p\text{ and }j=q
\\
0&,\text{other}
\end{cases}
$$
因此可以使用(7)组装迭代的雅可比矩阵，最终求得$K$后，带入(3)进行求解。

由于(2)的方程组一般非线性，实际求解时主要使用迭代方法求解，求解时注意以下几点:
① 迭代一般对初值敏感，$K$代表的是对应一阶导，这时可以直接使用(2)进行初始化，只是初始化用的$A,\mathbf{c}$均设置为0.
② 迭代一直不收敛时，将步长变小。

这里给出Radau IIA(3级5阶)的隐式求解系数^[2]，这个系数具有L稳定性,适合高频衰减和强刚性问题，具体系数如下
(8)
$$
\begin{cases}
A&=\begin{bmatrix}
\dfrac{88-7\sqrt{6}}{360}&\dfrac{296-169\sqrt{6}}{1800}&\dfrac{-2+3\sqrt{6}}{225}
\\
\\
\dfrac{296+169\sqrt{6}}{1800}&\dfrac{88+7\sqrt{6}}{360}&\dfrac{-2-3\sqrt{6}}{225}
\\
\\
\dfrac{16-\sqrt{6}}{36}&\dfrac{16+\sqrt{6}}{36}&\dfrac{1}{9}
\end{bmatrix}
\\
\\
\mathbf{c}&=\begin{bmatrix}\dfrac{4-\sqrt{6}}{10}&\dfrac{4+\sqrt{6}}{10}&1\end{bmatrix}^T
\\
\\
\mathbf{b}&=\begin{bmatrix}\dfrac{16-\sqrt{6}}{36}&\dfrac{16+\sqrt{6}}{36}&\dfrac{1}{9}\end{bmatrix}^T
\end{cases}
$$

另外，下面为3级6阶的高斯方法获得的系数，此系数具有辛性质(Symplectic)，适合哈密顿系统(天体力学、分子动力学)的长期积分
(9)
$$
\begin{cases}
A&=\begin{bmatrix}
\dfrac{5}{36}&\dfrac{2}{9}-\dfrac{\sqrt{15}}{15}&\dfrac{5}{36}-\dfrac{\sqrt{15}}{30}
\\
\\
\dfrac{5}{36}+\dfrac{\sqrt{15}}{24}&\dfrac{2}{9}&\dfrac{5}{36}-\dfrac{\sqrt{15}}{24}
\\
\\
\dfrac{5}{36}+\dfrac{\sqrt{15}}{30}&\dfrac{2}{9}+\dfrac{\sqrt{15}}{15}&\dfrac{5}{36}
\end{bmatrix}
\\
\\
\mathbf{c} &=\begin{bmatrix}\dfrac{1}{2}-\dfrac{\sqrt{15}}{10}&\dfrac{1}{2}&\dfrac{1}{2}+\dfrac{\sqrt{15}}{10}
\end{bmatrix}^T
\\
\\
\mathbf{b}&=\begin{bmatrix}\dfrac{5}{18}&\dfrac{4}{9}&\dfrac{5}{18}
\end{bmatrix}^T
\end{cases}
$$

2.4 辛性质
当布彻表里的系数满足如下关系时^[2]，则这套算法具有辛性质
(10)
$$
\mathbf{b}\mathbf{b}^T=\mathbf{diag}(\mathbf{b})A+A^T\mathbf{diag}(\mathbf{b})
$$

3. 应用
上面的算法已经集成到MathSword里，其中显式求解函数为ODESolverByRK45、ODESolverByRK56、ODESolverByRK67、ODESolverByRK78、ODESolverByRK89；隐式求解函数为ODESolverByRadauIIA35、ODESolverByGauss36、ODESolverByLobattoIII58、ODESolverByLobattoIIIA58、ODESolverByLobattoIIIB58、ODESolverByLobattoIIIC58。

对于显式求解，源码集成在开源框架《HeuristicSolver》里，使用时直接设置布彻表参数，就可以直接使用验证。

4. 参考
[1] Jim Verner’s Refuge for Runge-Kutta Pairs.
[2] Butcher J C. Numerical Methods for ordinary and differential equations[M]. John Wuley & Sons, 2016.

发表回复 取消回复

发表回复取消回复