Блог

 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, когда это произойдет. Я действительно не слишком хорошо знаю, как заставить функцию событий работать для вещей, выходящих за рамки производной, для которой я решаю

Если вы пытаетесь определить, когда второе значение 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        ]])]