如何让Python中solve_ivp采用均匀时间步长计算？

问题描述 #

需要按dt=0.001的固定时间步长求解抛体运动微分方程，虽然已将t_eval设为均匀分布的时间数组，但solve_ivp内部仍使用大小不一的步长（如1e-8秒）计算，而非指定的固定步长。简化代码如下：

t0 = 0
tf = 10
dt = 0.001

# def projectile_motion(t, y):
    # 微分方程实现
projectile_motion.t_prev = t0
sol = solve_ivp(projectile_motion, (t0, tf), [x0, y0, vx0, vy0], t_eval = np.linspace(t0, tf, int((tf-t0)/dt)+1))

核心问题 #

1. t_eval的作用局限：它仅控制输出结果的时间点，不影响求解器内部的计算步长。solve_ivp默认使用自适应步长算法（如RK45），会根据误差调整步长以满足精度要求。
2. 原函数逻辑错误：projectile_motion中手动用dt修改燃料剩余量、总质量等状态（如projectile_motion.m_fuel -= dt*m_fuel_duration），这种操作依赖固定步长，但solve_ivp的自适应步长会多次调用该函数，每次的t间隔不一定等于dt，导致状态计算混乱，甚至触发求解器缩小步长来捕捉异常变化。

解决方案 #

步骤1：修正微分方程函数 #

移除所有依赖固定dt的手动状态修改，改为基于当前时间t计算状态：

def projectile_motion(t, y):
    global check, x_of_fuel_out, x_hit, check_hit
    # 初始化状态（仅第一次调用时执行）
    if 't_start' not in projectile_motion.__dict__:
        projectile_motion.t_start = t
        projectile_motion.initial_m_total = m_total
        projectile_motion.initial_m_fuel = m_fuel
    
    x, y_pos, vx, vy = y
    v = np.sqrt(vx**2 + vy**2)
    
    # 计算燃料燃烧情况
    burn_time = min(t - projectile_motion.t_start, fuel_duration)
    m_fuel_used = m_fuel_duration * burn_time
    m_fuel_remaining = max(projectile_motion.initial_m_fuel - m_fuel_used, 0)
    current_m_total = projectile_motion.initial_m_total - (projectile_motion.initial_m_fuel - m_fuel_remaining)
    
    # 计算受力
    F_gravity = np.array([0, -current_m_total * g])
    if v > 1e-6:  # 避免除以0
        F_drag = -0.5 * rho * v**2 * Cd * A * np.array([vx, vy])/v 
    else:
        F_drag = np.array([0, 0])
    
    if m_fuel_remaining > 0:
        # 计算推力（基于当前燃料消耗率）
        F_thrust = v_e * m_fuel_duration * np.array([np.cos(theta), np.sin(theta)])
    else:
        if check == 0:
            check = 1
            x_of_fuel_out = x
        F_thrust = np.array([0, 0])
    
    # 检测落地
    if y_pos < 0 and check_hit == 0:
        check_hit = 1
        x_hit = x
    
    # 计算加速度
    F_total = F_gravity + F_drag + F_thrust
    a = F_total / current_m_total
    
    # 打印信息（可选）
    print(f"t: {t:.6f}      x: {round(x,2)},     v: {round(vx,2)}, {round(vy,2)} ,     thrust: {round(F_thrust[0],2)}, {round(F_thrust[1],2)},     drag: {round(F_drag[0],2)}, {round(F_drag[1],2)}     fuel left: {m_fuel_remaining:.4f}")
    
    return [vx, vy, a[0], a[1]]

步骤2：强制solve_ivp使用固定步长 #

调用solve_ivp时，添加max_step=dt参数，限制求解器的内部步长不超过dt，结合t_eval的均匀数组，即可实现固定步长计算：

t_eval = np.linspace(t0, tf, int((tf-t0)/dt)+1)
sol = solve_ivp(projectile_motion, (t0, tf), [x0, y0, vx0, vy0], 
                t_eval=t_eval, 
                max_step=dt,
                method='RK45')  # 可选，默认就是RK45

说明 #

• max_step=dt强制求解器的内部计算步长不超过dt，配合t_eval的均匀时间点，求解器会在每个t_eval点输出结果，且内部计算严格按等于dt的步长推进。
• 修正后的微分方程函数是纯函数，仅依赖输入的t和y计算导数，符合solve_ivp的调用规范，避免了状态混乱。

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