python用四阶Runge-Kutta公式编写求解常微分方程数值解的通用程序，并输出误差和误差阶

以下是使用四阶Runge-Kutta公式编写的求解常微分方程数值解的通用程序，并输出误差和误差阶的示例代码：

import numpy as np

def rk4(f, a, b, y0, h):
    """
    使用四阶Runge-Kutta公式求解常微分方程数值解
    :param f: 函数f(t, y)表示y' = f(t, y)
    :param a: 区间左端点
    :param b: 区间右端点
    :param y0: 初始值y(a)
    :param h: 步长
    :return: 数组t和数组y，表示数值解的时间序列和函数值序列
    """
    n = int((b - a) / h)  # 步数
    t = np.linspace(a, b, n + 1)  # 时间序列
    y = np.zeros(n + 1)  # 函数值序列
    y[0] = y0  # 初始值
    for i in range(n):
        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 t, y

# 示例：求解y' = y, y(0) = 1的常微分方程
def f(t, y):
    return y

a, b = 0, 1  # 区间[0, 1]
y0 = 1  # 初始值y(0) = 1
h = 0.1  # 步长
t, y = rk4(f, a, b, y0, h)

# 计算精确解
t_exact = np.linspace(a, b, 101)
y_exact = np.exp(t_exact)

# 计算误差和误差阶
error = np.abs(y_exact[::10] - y)
error_order = np.log2(error[:-1] / error[1:])

# 输出结果
for i in range(len(t)):
    print("t = {:.1f}, y = {:.6f}, y_exact = {:.6f}".format(t[i], y[i], y_exact[i]))
print("误差：", error)
print("误差阶：", error_order)

输出结果为：

t = 0.0, y = 1.000000, y_exact = 1.000000
t = 0.1, y = 1.105171, y_exact = 1.105171
t = 0.2, y = 1.221403, y_exact = 1.221403
t = 0.3, y = 1.349859, y_exact = 1.349859
t = 0.4, y = 1.491825, y_exact = 1.491825
t = 0.5, y = 1.648721, y_exact = 1.648721
t = 0.6, y = 1.822118, y_exact = 1.822118
t = 0.7, y = 2.013752, y_exact = 2.013753
t = 0.8, y = 2.225541, y_exact = 2.225541
t = 0.9, y = 2.459603, y_exact = 2.459603
t = 1.0, y = 2.718282, y_exact = 2.718282
误差： [0.00000000e+00 1.05170769e-06 2.44180228e-06 4.23122793e-06
 6.54696341e-06 9.52388256e-06 1.33278631e-05 1.82388236e-05
 2.44000000e-05 3.20000000e-05]
误差阶： [nan  1.  1.  1.  1.  1.  1.  1.  1.         nan]