如何使用scipy solve_ivp的events功能求解X=0.5对应的W、y值

solve_ivp的events参数要求传入的函数返回一个数值，求解器会自动检测该数值过零点的时刻，也就是你要的事件触发时刻。要捕获状态量X=0.5的事件，只需要让事件函数返回当前X值 - 0.5即可。

原代码问题 #

你写的事件函数固定返回0.5，没有和当前求解的X值绑定，所以永远不会触发过零点检测，事件不会生效。另外原代码中kprime、Fa0、epsilon、alpha四个参数没有定义，运行前需要先赋值为你的实际参数值。

修改后的完整代码 #

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 请先补充你的实际参数值
kprime = 1.0  # 示例值，替换为你的实际值
Fa0 = 1.0     # 示例值，替换为你的实际值
epsilon = 0.0 # 示例值，替换为你的实际值
alpha = 0.01  # 示例值，替换为你的实际值

def dYdW(W, Y):
    X, y = Y
    dXdW = (kprime/Fa0)*(1-X)*y*(1/(1+epsilon*X))
    dydW = -1*alpha*(1+epsilon*X)*(1/2*y)
    return np.array([dXdW, dydW])

# 正确的事件函数
def event(W, Y): 
    X, y = Y
    # 返回X-0.5，当X=0.5时返回值为0，触发事件
    return X - 0.5

# 可选配置：触发事件后是否终止求解，默认不终止，需要的话取消注释下面这句
# event.terminal = True

X0 = np.array([0, 1])
Wspan = np.array([0, 100])
Weval, h = np.linspace(*Wspan, 500, retstep=True)

sol = solve_ivp(dYdW, Wspan, X0, max_step=h, events=event)

# 提取X=0.5时的结果
if len(sol.t_events[0]) > 0:
    W_x05 = sol.t_events[0][0]
    X_x05, y_x05 = sol.y_events[0][0]
    print(f"X=0.5时，对应W={W_x05:.4f}，对应y={y_x05:.4f}")
else:
    print("在指定的Wspan范围内X未达到0.5")

# 绘图修正
plt.plot(sol.t, sol.y.T)
plt.xlabel("W")
plt.ylabel("取值")
plt.legend(["转化率X", "无量纲压力y"])
plt.show()

结果说明 #

• 如果配置了event.terminal = True，求解器会在X第一次达到0.5时就停止计算，适合只需要阈值点结果的场景
• sol.t_events是一个列表，每个元素对应一个事件函数的触发时刻数组，这里只有一个事件，所以取sol.t_events[0]获取所有X=0.5的时刻
• sol.y_events结构和t_events对应，存储触发时刻的状态变量值

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