scipy.integrate.solve_ivp事件函数使用报错及积分终止实现求助

你的代码触发TypeError且无法正确终止积分，核心问题出在事件函数的定义不符合solve_ivp的要求，以下是具体修正方案：

错误原因拆解 #

1. 事件函数返回值错误：solve_ivp要求事件函数返回单个标量，当该值穿越0时触发事件；你返回了修改后的状态数组x，完全违反规则。
2. 标量函数处理数组：用math.sqrt处理numpy数组会触发类型错误，必须使用numpy的数组兼容函数。
3. 未设置终止属性：没有给事件函数添加terminal=True，无法触发积分终止。

正确事件函数实现 #

定义一个返回标量差值的事件函数，同时添加终止属性：

import numpy as np

def position_exceed_event(t, x):
    # 计算前3个元素（位置）的模长
    r = np.linalg.norm(x[:3])
    # 返回阈值与模长的差值：当模长超过10000时，值从正变负，触发事件
    return 10000 - r

# 标记事件触发时终止积分
position_exceed_event.terminal = True
# 可选：仅当差值从正转负时触发（即模长从小于10000变为大于）
position_exceed_event.direction = -1

修正后的调用代码 #

将正确的事件函数传入events参数：

sol = solve_ivp(
    fun=lambda t, x: my_eq(t, x, param),
    t_span=integr_time,
    y0=np.concatenate((r, v)),
    method='LSODA',
    events=position_exceed_event,
    atol=1e-12,
    rtol=1e-13
)

额外说明 #

• np.linalg.norm(x[:3])是计算三维向量模长的最优方式，比手动计算平方和开根号更简洁可靠。
• direction=-1可以避免模长从大于阈值回落时误触发事件，若不需要这个限制可省略。

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