Python实现四阶标准R-K方法求解二阶常微分方程初值问题

首先，将二阶常微分方程转化为一阶方程组：

$
\begin{cases}
y_1' = y_2 \\
y_2' = f(x, y_1, y_2)
\end{cases}
$

其中，$y_1$ 是未知函数，$y_2$ 是 $y_1$ 的导数，$f(x, y_1, y_2)$ 是已知函数。

然后，可以使用四阶标准 R-K 方法求解该方程初值问题，具体步骤如下：

1. 定义步长 $h$，初始值 $x_0$ 和 $y_{10}$，$y_{20}$；
2. 计算 $k_{11}, k_{21}, k_{12}, k_{22}$：

$
\begin{aligned}
k_{11} &= h y_{2n} \\
k_{21} &= hf(x_n, y_{1n}, y_{2n}) \\
k_{12} &= h(y_{2n} + \frac{k_{21}}{2}) \\
k_{22} &= hf(x_n + \frac{h}{2}, y_{1n} + \frac{k_{11}}{2}, y_{2n} + \frac{k_{21}}{2})
\end{aligned}
$

3. 计算 $k_{31}, k_{41}, k_{32}, k_{42}$：

$
\begin{aligned}
k_{31} &= h(y_{2n} + \frac{k_{22}}{2}) \\
k_{41} &= hf(x_n + \frac{h}{2}, y_{1n} + \frac{k_{31}}{2}, y_{2n} + \frac{k_{32}}{2}) \\
k_{32} &= h(y_{2n} + \frac{k_{31}}{2}) \\
k_{42} &= hf(x_n + \frac{h}{2}, y_{1n} + \frac{k_{32}}{2}, y_{2n} + \frac{k_{42}}{2})
\end{aligned}
$

4. 计算 $y_{1n+1}, y_{2n+1}$：

$
\begin{aligned}
y_{1n+1} &= y_{1n} + \frac{1}{6}(k_{11} + 2k_{21} + 2k_{31} + k_{41}) \\
y_{2n+1} &= y_{2n} + \frac{1}{6}(k_{12} + 2k_{22} + 2k_{32} + k_{42})
\end{aligned}
$

5. 更新 $x_n$ 为 $x_n + h$；
6. 重复步骤 2~5，直到 $x_n$ 达到所求的 $x$。

以下是 Python 代码实现：

def rk4(f, x0, y0, z0, h, x):
    n = int((x - x0) / h)
    y = [y0]
    z = [z0]
    for i in range(n):
        xn = x0 + i * h
        yn = y[-1]
        zn = z[-1]
        k1 = h * zn
        l1 = h * f(xn, yn, zn)
        k2 = h * (zn + l1 / 2)
        l2 = h * f(xn + h / 2, yn + k1 / 2, zn + l1 / 2)
        k3 = h * (zn + l2 / 2)
        l3 = h * f(xn + h / 2, yn + k2 / 2, zn + l2 / 2)
        k4 = h * (zn + l3)
        l4 = h * f(xn + h, yn + k3, zn + l3)
        y.append(yn + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4))
        z.append(zn + 1 / 6 * (l1 + 2 * l2 + 2 * l3 + l4))
    return y, z

其中，$f(x, y_1, y_2)$ 即为二阶常微分方程右端的已知函数，$x_0, y_{10}, y_{20}$ 分别为初值问题的初始值，$h$ 为步长，$x$ 为所求的 $x$。函数返回值为 $y_n$ 和 $y_n'$ 的数组。