python用四阶Runge-Kutta公式编写求解常微分方程数值解的通用程序

可以按照以下步骤，使用 Python 编写求解常微分方程数值解的通用程序：

1. 导入必要的库和函数：

import numpy as np
import matplotlib.pyplot as plt

2. 定义四阶 Runge-Kutta 方法的函数：

def rk4(f, t0, y0, h, n):
    t = np.zeros(n+1)
    y = np.zeros((n+1, len(y0)))
    t[0] = t0
    y[0] = y0
    for i in range(n):
        k1 = h*f(t[i], y[i])
        k2 = h*f(t[i]+0.5*h, y[i]+0.5*k1)
        k3 = h*f(t[i]+0.5*h, y[i]+0.5*k2)
        k4 = h*f(t[i]+h, y[i]+k3)
        y[i+1] = y[i] + (k1+2*k2+2*k3+k4)/6
        t[i+1] = t[i] + h
    return t, y

其中，参数 `f` 是常微分方程的右端函数，即 $f(t,y)$；参数 `t0` 和 `y0` 分别是初始时刻和初始值；参数 `h` 是步长；参数 `n` 是时间步数。

3. 定义需要求解的常微分方程的右端函数，例如：

def f(t, y):
    return np.array([y[1], -y[0]])

4. 设置初始时刻、初始值、步长和时间步数：

t0 = 0
y0 = np.array([0, 1])
h = 0.1
n = 100

5. 调用 `rk4` 函数求解常微分方程的数值解：

t, y = rk4(f, t0, y0, h, n)

6. 绘制数值解的图像：

plt.plot(t, y[:, 0], label='y(t)')
plt.plot(t, y[:, 1], label='y\'(t)')
plt.legend()
plt.show()

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt

def rk4(f, t0, y0, h, n):
    t = np.zeros(n+1)
    y = np.zeros((n+1, len(y0)))
    t[0] = t0
    y[0] = y0
    for i in range(n):
        k1 = h*f(t[i], y[i])
        k2 = h*f(t[i]+0.5*h, y[i]+0.5*k1)
        k3 = h*f(t[i]+0.5*h, y[i]+0.5*k2)
        k4 = h*f(t[i]+h, y[i]+k3)
        y[i+1] = y[i] + (k1+2*k2+2*k3+k4)/6
        t[i+1] = t[i] + h
    return t, y

def f(t, y):
    return np.array([y[1], -y[0]])

t0 = 0
y0 = np.array([0, 1])
h = 0.1
n = 100

t, y = rk4(f, t0, y0, h, n)

plt.plot(t, y[:, 0], label='y(t)')
plt.plot(t, y[:, 1], label='y\'(t)')
plt.legend()
plt.show()