四阶龙格库塔法求方程组

四阶龙格库塔法（RK4）是一种常用的求解常微分方程组的数值方法，下面给出其求解一个二阶常微分方程组的示例。

考虑方程组：

$$
\begin{cases}
y_1' = y_2 \\
y_2' = -\omega^2 y_1 - \gamma y_2 + F(t)
\end{cases}
$$

其中 $y_1$ 和 $y_2$ 分别是时间 $t$ 的函数，$\omega$、$\gamma$ 和 $F(t)$ 是常数或者已知函数。

我们将该方程组转化为向量形式：

$$
\mathbf{y}' = \begin{pmatrix} y_1' \\ y_2' \end{pmatrix} = \begin{pmatrix} y_2 \\ -\omega^2 y_1 - \gamma y_2 + F(t) \end{pmatrix} = \mathbf{f}(t,\mathbf{y})
$$

其中 $\mathbf{y} = \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}$，$\mathbf{f}$ 是一个向量值函数。

RK4 方法的基本思想是在每个步长内做四次计算，得到一个时间步长的近似解。具体步骤如下：

1. 用初始条件 $\mathbf{y}_0$，计算 $\mathbf{y}_1$：

$$
\mathbf{k}_1 = h\mathbf{f}(t_0,\mathbf{y}_0), \quad \mathbf{k}_2 = h\mathbf{f}(t_0+\frac{h}{2},\mathbf{y}_0+\frac{1}{2}\mathbf{k}_1), \quad \mathbf{k}_3 = h\mathbf{f}(t_0+\frac{h}{2},\mathbf{y}_0+\frac{1}{2}\mathbf{k}_2), \quad \mathbf{k}_4 = h\mathbf{f}(t_0+h,\mathbf{y}_0+\mathbf{k}_3)
$$

$$
\mathbf{y}_1 = \mathbf{y}_0 + \frac{1}{6}(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4)
$$

2. 用 $\mathbf{y}_1$，计算 $\mathbf{y}_2$：

$$
\mathbf{k}_1 = h\mathbf{f}(t_1,\mathbf{y}_1), \quad \mathbf{k}_2 = h\mathbf{f}(t_1+\frac{h}{2},\mathbf{y}_1+\frac{1}{2}\mathbf{k}_1), \quad \mathbf{k}_3 = h\mathbf{f}(t_1+\frac{h}{2},\mathbf{y}_1+\frac{1}{2}\mathbf{k}_2), \quad \mathbf{k}_4 = h\mathbf{f}(t_1+h,\mathbf{y}_1+\mathbf{k}_3)
$$

$$
\mathbf{y}_2 = \mathbf{y}_1 + \frac{1}{6}(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4)
$$

重复以上步骤，直到计算到所需的时间点。

下面是一个 Python 代码示例：

import numpy as np

def rk4(f, y0, t):
    """
    使用 RK4 方法求解二阶常微分方程组
    
    参数：
        f: 函数，右端函数，输入参数为 (t, y)，返回值为 y'。
        y0: 数组，初始条件。
        t: 数组，时间点。
    返回：
        y: 数组，每个时间点的近似解。
    """
    h = t[1] - t[0]
    y = np.zeros((len(t), len(y0)))
    y[0] = y0
    
    for i in range(len(t) - 1):
        k1 = h * f(t[i], y[i])
        k2 = h * f(t[i] + h/2, y[i] + k1/2)
        k3 = h * f(t[i] + h/2, y[i] + k2/2)
        k4 = h * f(t[i] + h, y[i] + k3)
        y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
        
    return y

使用方法：

# 定义右端函数
def f(t, y):
    w = 1
    gamma = 0.1
    F = 1
    return np.array([y[1], -w**2*y[0]-gamma*y[1]+F])

# 定义时间点
t = np.linspace(0, 10, 101)

# 定义初始条件
y0 = np.array([0, 1])

# 求解方程组
y = rk4(f, y0, t)

# 绘图
import matplotlib.pyplot as plt
plt.plot(t, y[:, 0], label='y1')
plt.plot(t, y[:, 1], label='y2')
plt.legend()
plt.show()

注意，RK4 方法并不一定能够保证数值解的精度，当步长 $h$ 过大时，数值解的精度可能会受到影响。因此，实际使用中需要根据具体问题来选择合适的步长。