如何在NumPy中实现欧拉法、龙格-库塔法、中点法的全区间积分函数？

看起来你已经有了清晰的框架思路，我来帮你补全并梳理这个数值积分实现，完全符合你要求的单步计算函数+全区间统一积分函数的结构：

核心需求回顾 #

• 三个独立的单步计算函数：euler、rk4、middlepoint，仅完成单次步长的积分计算并返回结果
• 一个统一的integrateall函数，根据指定的方法，完成整个时间区间的积分过程
• 基于NumPy实现，保证数值计算的高效性

完整代码实现 #

import numpy as np
import matplotlib.pyplot as plt  # 修正原代码的导入错误
from scipy.integrate import odeint

def euler(y, t, h, ode):
    """欧拉法单步计算"""
    return y + h * ode(y, t)

def middlepoint(y, t, h, ode):
    """中点法（改进欧拉法的一种）单步计算"""
    k1 = ode(y, t)
    y_mid = y + (h / 2) * k1
    t_mid = t + (h / 2)
    k2 = ode(y_mid, t_mid)
    return y + h * k2

def rk4(y, t, h, ode):
    """经典龙格-库塔4阶法单步计算"""
    k1 = ode(y, t)
    k2 = ode(y + (h / 2) * k1, t + (h / 2))
    k3 = ode(y + (h / 2) * k2, t + (h / 2))
    k4 = ode(y + h * k3, t + h)
    return y + (h / 6) * (k1 + 2*k2 + 2*k3 + k4)

def integrateall(y0, t_span, h, ode, method='rk4'):
    """
    全区间数值积分函数
    参数:
        y0: 初始值（标量或NumPy数组）
        t_span: 时间区间，格式为[start, end]
        h: 积分步长
        ode: 待积分的常微分方程，格式为dy/dt = ode(y, t)
        method: 积分方法，可选'euler'、'middlepoint'、'rk4'，默认rk4
    返回:
        t: 时间点数组
        y: 对应时间点的积分结果数组
    """
    # 生成时间点数组
    t_start, t_end = t_span
    t = np.arange(t_start, t_end + h, h)
    # 初始化结果数组
    y = np.zeros((len(t), len(np.atleast_1d(y0))))
    y[0] = y0
    # 根据选择的方法执行单步循环
    step_func = {
        'euler': euler,
        'middlepoint': middlepoint,
        'rk4': rk4
    }[method]
    
    for i in range(len(t) - 1):
        y[i+1] = step_func(y[i], t[i], h, ode)
    
    # 如果初始值是标量，返回一维数组
    if np.isscalar(y0):
        y = y.flatten()
    
    return t, y

函数细节说明 #

单步计算函数 #

每个单步函数都严格遵循输入输出规范：接受当前状态y、当前时间t、步长h、ODE函数ode，返回下一步的状态值：

• 欧拉法：最简单的一阶数值方法，计算量小但精度较低
• 中点法：通过计算区间中点的导数来更新状态，精度比欧拉法高，属于二阶方法
• RK4法：经典四阶龙格-库塔法，通过四个不同位置的导数加权平均更新状态，精度高，是工程中常用的方法

全区间积分函数integrateall #

• 自动生成均匀步长的时间点数组
• 根据指定的方法映射到对应的单步函数
• 循环执行单步计算，收集所有时间点的结果
• 兼容标量和数组形式的初始值，返回对应维度的结果

示例用法：对比不同方法的精度 #

我们用一个简单的常微分方程$\frac{dy}{dt} = -y$（初始值$y(0)=1$）来测试，同时对比scipy.odeint的结果：

# 定义待积分的ODE
def ode(y, t):
    return -y

# 计算参数
y0 = 1
t_span = [0, 5]
h = 0.1

# 用三种方法计算
t_euler, y_euler = integrateall(y0, t_span, h, ode, method='euler')
t_mid, y_mid = integrateall(y0, t_span, h, ode, method='middlepoint')
t_rk4, y_rk4 = integrateall(y0, t_span, h, ode, method='rk4')

# 用scipy的odeint作为参考
t_ref = np.linspace(t_span[0], t_span[1], 100)
y_ref = odeint(ode, y0, t_ref)

# 绘图对比
plt.figure(figsize=(10, 6))
plt.plot(t_ref, y_ref, 'k--', label='Exact (odeint)')
plt.plot(t_euler, y_euler, 'r-', label='Euler Method')
plt.plot(t_mid, y_mid, 'g-', label='Middlepoint Method')
plt.plot(t_rk4, y_rk4, 'b-', label='RK4 Method')
plt.xlabel('Time t')
plt.ylabel('y(t)')
plt.title('Comparison of Numerical Integration Methods')
plt.legend()
plt.grid(True)
plt.show()

运行这段代码后，你可以清晰看到不同方法的精度差异：欧拉法误差最大，中点法次之，RK4法几乎和精确解重合。

内容的提问来源于stack exchange，提问作者ZelelB