Scipy solve_ivp：能否在运行过程中修改max_step参数？

可以实现动态调整max_step，核心思路是分阶段积分或使用更灵活的ode类，以下是具体方案：

方法1：分阶段积分（基于solve_ivp事件机制） #

solve_ivp本身不支持在积分过程中直接修改参数，但可以通过事件函数检测弹性体与表面的距离，触发积分终止，再切换步长重新启动积分，循环完成整个过程。

步骤： #

1. 定义两个事件函数：一个检测弹性体接近表面，另一个检测离开表面

def near_surface(t, y):
    # y[0]为弹性体到表面的距离，threshold是你定义的"靠近"阈值
    return y[0] - 0.1  # 示例阈值0.1
near_surface.terminal = True
near_surface.direction = -1  # 仅当距离从大于阈值降到小于时触发

def far_from_surface(t, y):
    return y[0] - 0.2  # 示例离开阈值0.2
far_from_surface.terminal = True
far_from_surface.direction = 1  # 仅当距离从小于阈值升到大于时触发

2. 分阶段执行积分：

from scipy.integrate import solve_ivp

t_span = [0, 10]  # 总积分时间
y0 = [10, 0]  # 初始距离、初始速度
t_eval = []
y_eval = []

current_t = t_span[0]
current_y = y0
large_step = 0.01  # 远离表面时的大步长
small_step = 0.0001  # 靠近表面时的小步长

while current_t < t_span[1]:
    # 第一阶段：大步长积分到靠近表面
    sol1 = solve_ivp(your_ode_func, [current_t, t_span[1]], current_y, 
                     method='Radau', max_step=large_step, events=near_surface)
    # 保存第一阶段结果
    t_eval.extend(sol1.t)
    y_eval.extend(sol1.y.T)
    
    if sol1.status == 1:  # 触发了靠近表面事件
        current_t = sol1.t_events[0][0]
        current_y = sol1.y_events[0][0]
        # 第二阶段：小步长积分到离开表面
        sol2 = solve_ivp(your_ode_func, [current_t, t_span[1]], current_y, 
                         method='Radau', max_step=small_step, events=far_from_surface)
        t_eval.extend(sol2.t[1:])  # 避免重复保存触发点
        y_eval.extend(sol2.y.T[1:])
        current_t = sol2.t[-1]
        current_y = sol2.y[:, -1]
    else:
        # 未触发事件，直接完成积分
        break

# 最终结果整理为数组
import numpy as np
t_eval = np.array(t_eval)
y_eval = np.array(y_eval)

方法2：使用scipy.integrate.ode类（更灵活的步长控制） #

如果不想分阶段调用solve_ivp，可以使用更底层的ode类，它支持在每一步积分前动态调整max_step：

from scipy.integrate import ode

r = ode(your_ode_func).set_integrator('radau')
r.set_initial_value(y0, t_span[0])

threshold_near = 0.1
threshold_far = 0.2
max_step_large = 0.01
max_step_small = 0.0001

t_eval = [t_span[0]]
y_eval = [y0]

while r.successful() and r.t < t_span[1]:
    current_distance = r.y[0]
    # 根据距离动态设置max_step
    if current_distance < threshold_near:
        r.set_integrator('radau', max_step=max_step_small)
    elif current_distance > threshold_far:
        r.set_integrator('radau', max_step=max_step_large)
    # 执行一步积分，步长不超过当前设置的max_step和剩余时间
    next_t = min(r.t + (max_step_small if current_distance < threshold_near else max_step_large), t_span[1])
    r.integrate(next_t)
    t_eval.append(r.t)
    y_eval.append(r.y.copy())

t_eval = np.array(t_eval)
y_eval = np.array(y_eval)

注意事项 #

• 阈值的选择需要根据你的物理模型调整，确保在弹性体进入碰撞交互区域前切换到小步长
• 使用事件函数时，direction参数要正确设置，避免在阈值附近反复触发终止
• 分阶段积分时，注意不要重复保存事件触发点的状态（比如示例中sol2.t[1:]）

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