四阶龙格库塔法求解一次常微分方程组(python实现)

一、前言

之前在博客发布了关于使用四阶龙格库塔方法求解一次常微分方程组的文章，由于代码缺少具体的验证，部分朋友可能存在疑问，因此这里打算再重新写一篇博客来验证一下程序的正确性，另外，这里是使用python语言来实现的。

二、RK4求解方程组的要点

使用RK4求解一元方程的过程是非常容易的，但是当转变成多变量的情况下，如何求解方程组，可能有部分朋友会出现问题，这里我总结出来主要有以下要点。

1. 将方程组转化为RK4求解要求的标准形式

RK4求解要求方程具有形如$\dot{y}=f(y,t)$的形式，即是在方程的左边只有变量的一阶微分项，方程的右边包含变量项和自变量项，对于n元情况，我们需要有n个方程，每个方程对应一个自变量项，形式如下：
$$\{\begin{array}{l}\dot{y_1}=f_1(y_1,\cdots ,y_n,t)\\ ⋮\\ \dot{y_n}=f_n(y_1,\cdots ,y_n,t)\end{array}$$
转化为标准形式之后就可以进行求解了。

2. 注意区分每个方程的独立性

RK4求解是通过指定一个较小的步进距离，来逐步求解前进一步之后的函数值，每一步下的函数值求解都需要用到前一步的结果，属于递推过程。对于方程组的求解过程，独立性是指针对某个特定变量时，递推公式中只改变特定变量的递推关系，而其它变量不变，例如，方程组中$y_i$的第m+1项的RK4递推关系可以写作：
{ h m = t m + 1 − t m k 1 i m = f i ( y 1 m , ⋯   , y i m , ⋯   , y n m , t m ) k 2 i m = f i ( y 1 m , ⋯   , y i m + h m 2 k 1 i m , ⋯   , y n m , t m + h m 2 ) k 3 i m = f i ( y 1 m , ⋯   , y i m + h m 2 k 2 i m , ⋯   , y n m , t m + h m 2 ) k 4 i m = f i ( y 1 m , ⋯   , y i m + h m k 3 i m , ⋯   , y n m , t m + h m ) y i m + 1 = y i m + h m 6 ( k 1 i m + 2 k 2 i m + 2 k 3 i m + k 4 i m ) \begin{cases} h_m = t_{m+1}-t_m \\ k_{1i}^{m} = f_i(y_1^m, \cdots, y_i^m,\cdots, y_n^m, t_m)\\ k_{2i}^{m} = f_i(y_1^m, \cdots,y_i^m+\frac{h_m}{2}k_{1i}^m,\cdots,y_n^m, t_m+\frac{h_m}{2})\\ k_{3i}^{m} = f_i(y_1^m,\cdots,y_i^m+\frac{h_m}{2}k_{2i}^m,\cdots,y_n^m,t_m+\frac{h_m}{2})\\ k_{4i}^{m} = f_i(y_1^m,\cdots,y_i^m+h_mk_{3i}^m,\cdots,y_n^m,t_m+h_m)\\ y_i^{m+1} = y_i^m+\frac{h_m}{6}(k_{1i}^m+2k_{2i}^m+2k_{3i}^m+k_{4i}^m) \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪