Как я могу контролировать odeint, чтобы остановить интеграцию, когда результат достигнет порога?

Question

Вот мой код.

import numpy as np
from scipy.integrate import odeint

#Constant
R0=1.475
gamma=2.
ScaleMeVfm3toEskm3 = 8.92*np.power(10.,-7.)

def EOSe(p):
    return np.power((p/450.785),(1./gamma))

def M(m,r):
    return (4./3.)*np.pi*np.power(r,3.)*p

# function that returns dz/dt
def model(z,r):
    p, m = z
    dpdr = -((R0*EOSe(p)*m)/(np.power(r,2.)))*(1+(p/EOSe(p)))*(1+((4*math.pi*(np.power(r,3))*p)/(m)))*((1-((2*R0)*m)/(r))**(-1.))
    dmdr = 4.*math.pi*(r**2.)*EOSe(p)
    dzdr = [dpdr,dmdr]
    return dzdr

# initial condition
r0=10.**-12.
p0=10**-6.
z0 = [p0, M(r0, p0)]

# radius
r = np.linspace(r0, 15, 100000)

# solve ODE
z = odeint(model,z0,r)

Результат z[:,0] продолжает уменьшаться, как я и ожидал. Но я хочу только положительных ценностей. Можно запустить код и попробовать print(z[69306]), и он покажет [2.89636405e-11 5.46983202e-01]. Это последнее, что я хочу, чтобы odeint прекратил интеграцию.

Конечно, приведенный код показывает

RuntimeWarning: invalid value encountered in power
  return np.power((p/450.785),(1./gamma))

потому что результат p начинает быть отрицательным. Для любых дальнейших пунктов odeint дает результат [nan nan].

Однако я мог бы использовать np.nanmin (), чтобы найти минимум z [:, 0], который не является nan. Но у меня есть набор p0 значений для моей работы. Мне нужно будет позвонить odeint в цикле, как

P=np.linspace(10**-8.,10**-2.,10000)
for p0 in P:
#the code for solving ode provided above.

, что занимает больше времени.

Я думаю, что это уменьшило бы время выполнения, если бы я мог просто остановиться до того, как z [:, 0] станет отрицательным значением?

Answer 1

Вот модифицированный код, использующий solve_ivp:

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

# Constants
R0 = 1.475
gamma = 2.

def EOSe(p):
    return np.power(np.abs(p)/450.785, 1./gamma)

def M(m, r):
    return (4./3.)*np.pi*np.power(r,3.)*p

# function that returns dz/dt
# note: the argument order is reversed compared to `odeint`
def model(r, z):
    p, m = z
    dpdr = -R0*EOSe(p)*m/r**2*(1 + p/EOSe(p))*(1 + 4*np.pi*r**3*p/m)*(1 - 2*R0*m/r)**(-1)
    dmdr = 4*np.pi * r**2 * EOSe(p)
    dzdr = [dpdr, dmdr]
    return dzdr

# initial condition
r0 = 1e-3
r_max = 50
p0 = 1e-6
z0 = [p0, M(r0, p0)]

# Define the event function
# from the doc: "The solver will find an accurate value
# of t at which event(t, y(t)) = 0 using a root-finding algorithm. "
def stop_condition(r, z):
    return z[0]

stop_condition.terminal = True

# solve ODE
r_span = (r0, r_max)
sol = solve_ivp(model, r_span, z0,
                events=stop_condition)


print(sol.message)
print('last p, m = ', sol.y[:, -1], 'for r_event=', sol.t_events[0][0])

r_sol = sol.t
p_sol = sol.y[0, :]
m_sol = sol.y[1, :]

# Graph
plt.subplot(2, 1, 1);
plt.plot(r_sol, p_sol, '.-b')
plt.xlabel('r'); plt.ylabel('p');

plt.subplot(2, 1, 2);
plt.plot(r_sol, m_sol, '.-r')
plt.xlabel('r'); plt.ylabel('m');

На самом деле, использование событий в этом случае не предотвращает предупреждение из-за отрицательного значения p.Причина в том, что решатель все равно собирается оценить модель на p<O.Решением является получение абсолютного значения p в квадратном корне (как в коде выше).Использование np.sign(p)*np.power(np.abs(p)/450.785, 1./gamma) также дает интересный результат.