scipy solve_ivp events：导数积分中的事件

Posted 2023-03-24

【中文标题】scipy solve_ivp events：导数积分中的事件

【英文标题】：scipy solve_ivp events : events in the integral of the derivative

【发布时间】：2021-11-18 01:54:58

【问题描述】：

import numpy as np
from scipy.optimize import fsolve
from scipy.integrate import quad
from scipy.integrate import solve_ivp

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
  return(0.5)

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)

import matplotlib.pyplot as plt
plt.plot(sol.t, sol.y.T)
plt.xlabel("w")
plt.ylabel("X (conversion, y(dimensionless pressure)")
plt.legend("Conversion","pressure")

sol.t_events

我正在尝试求解函数 X = 0.5 的时间，并确定发生这种情况时 W 和 y 的值。我真的不太清楚如何让事件函数为我正在解决的导数之外的事情工作

【参考方案1】：

如果您要查找 y 的第二个值何时等于 0.5，则必须在该值为真时返回 0。尝试这样的事情，我假设你的所有参数都是 1。

import numpy as np
from scipy.integrate import solve_ivp
kprime=1
Fa0=1
epsilon=1
alpha=1
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):
    return Y[1] - 0.5
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)
print(sol.t_events)
print(sol.y_events)

产量

[array([1.06945742])]
[array([[0.46528842, 0.5        ]])]