如何在scipy.solve_ivp中去除重复过零点并匹配向量符号

你遇到的是数值求解中常见的精度误差问题：solve_ivp检测到phi_dot=0事件时，返回的是接近0的极小值（如1.905e-13），这个值会同时出现在事件前后的结果数组里，导致后续数组符号不一致。可以通过两种方式调整：

方式1：事后修正结果数组 #

直接修改后续数组中那个过零点的符号，匹配数组内其他元素的符号：

import numpy as np
from scipy.integrate import solve_ivp

# 假设已完成求解，得到sol对象
# 定位事件发生的索引
event_time = sol.t_events[0][0]
event_idx = np.where(np.isclose(sol.t, event_time))[0][0]

# 提取事件后的phi_dot数组并修改
phi_dot_after = sol.y[1][event_idx:].copy()
# 取后续第一个非零元素的符号
sign = np.sign(phi_dot_after[1])
# 调整过零点的符号
phi_dot_after[0] = sign * abs(phi_dot_after[0])

# 替换原数组
sol.y[1][event_idx:] = phi_dot_after

方式2：调整事件函数与求解流程 #

让求解器在事件发生时终止，再手动以带符号的极小值作为初始条件继续积分，从根源避免符号不一致：

# 定义事件函数，设置终止规则与检测方向
def event(t, y):
    return y[1]
event.terminal = True  # 事件发生时终止求解
event.direction = -1   # 仅检测phi_dot从正到负的穿越（可按需调整方向）

# 第一次求解至事件点
sol1 = solve_ivp(your_ode_func, t_span, y0, events=event)

# 修改事件点的phi_dot为带符号的极小值，匹配后续趋势
new_y0 = sol1.y[:, -1].copy()
new_y0[1] = -1e-13  # 对应后续的负号

# 继续求解后续区间
sol2 = solve_ivp(your_ode_func, [sol1.t[-1], t_end], new_y0)

# 合并两次求解的结果
combined_t = np.concatenate([sol1.t, sol2.t[1:]])
combined_y = np.concatenate([sol1.y, sol2.y[:, 1:]], axis=1)

说明：数值计算中无法得到严格的0，那个极小值是精度限制导致的，调整符号符合物理上的趋势（事件后phi_dot为负），是合理的处理方式。

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