是否有方法访问scipy.integrate.odeint中的当前时间步长？
我试图解决一个ode系统，其中ode的形式取决于种群是否会耗尽。基本上，如果x不低于阈值，我从群体x中获得。如果我需要采取这个时间步长的量大于该阈值，则我将采取所有x到该点，其余的从z。
我正试图通过检查我将在这一时间步长中花费多少来实现这一点，然后在DE中的种群x和z之间进行分配。
要做到这一点，我需要能够访问ODE求解器中的步长，以计算这个时间步长将采取的措施。我正在使用scipy.integrate.odeint-是否有方法访问定义ode的函数中的时间步长？
或者，您可以访问解算器中的最后一次是什么吗？我知道这不一定是下一个时间步长，但如果这是我能做的最好的，对我来说可能是一个足够好的近似值。或者还有其他我没有想过的选择吗？
下面的MWE不是我的方程组，而是我可以想出的来说明我在做什么。问题是，在第一个时间步长上，如果时间步长为1，则总体将过低，但由于时间步长较小，最初可以从x中获取所有值。
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
plt.interactive(False)
tend = 5
tspan = np.linspace(0.0, tend, 1000)
A = 3
B = 4.09
C = 1.96
D = 2.29
def odefunc(P,t):
x = P[0]
y = P[1]
z = P[2]
if A * x - B * x * y < 0.6:
dxdt = A/5 * x
dydt = -C * y + D * x * y
dzdt = - B * z * y
else:
dxdt = A * x - B * x * y
dydt = -C * y + D * x * y
dzdt = 0
dPdt = np.ravel([dxdt, dydt, dzdt])
return dPdt
init = ([0.75,0.95,100])
sol = odeint(odefunc, init, tspan, hmax = 0.01)
x = sol[:, 0]
y = sol[:, 1]
z = sol[:, 2]
plt.figure(1)
plt.plot(tspan,x)
plt.plot(tspan,y)
plt.plot(tspan,z)
当然，您可以一起破解一些可能有效的东西。您可以记录t，但必须注意可能不会不断增加。这取决于ODE算法及其工作方式(前向、后向和中心有限差分(。
但它会让你知道你在哪里。
logger = [] # visible in odefunc
def odefunc(P,t):
x = P[0]
y = P[1]
z = P[2]
print(t)
logger.append(t)
if logger: # if the list is not empty
if logger[-1] > 2.5: # then read the last value
print('hua!')
if A * x - B * x * y < 0.6:
dxdt = A/5 * x
dydt = -C * y + D * x * y
dzdt = - B * z * y
else:
dxdt = A * x - B * x * y
dydt = -C * y + D * x * y
dzdt = 0
dPdt = np.ravel([dxdt, dydt, dzdt])
return dPdt
print(logger)
正如另一个答案中所指出的，在每次调用odeint中的ODE函数时，时间可能不会严格增加，尤其是对于棘手的问题。
在ode函数中处理这种不连续性的最稳健方法是使用事件来查找示例中(A * x - B * x * y) - 0.6的零的位置。对于不连续的解决方案，使用终端事件精确地在零处停止计算，然后更改ode函数。在solve_ivp中，可以使用events参数执行此操作。请参阅解决ivp文档，特别是与炮弹轨迹相关的示例。odeint不支持事件，solve_ivp有一个LSODA方法，它调用与odeint相同的Fortran库。
这里是一个简短的示例，但您可能需要在解决sol2之前，额外检查sol1是否到达了终端事件。
from scipy.integrate import solve_ivp
tend = 10
def discontinuity_zero(t, y):
return y[0] - 10
discontinuity_zero.terminal = True
def ode_func1(t, y):
return y
def ode_func2 (t, y):
return -y**2
sol1 = solve_ivp(ode_func1, t_span=[0, tend], y0=[1], events=discontinuity_zero, rtol=1e-8)
t1 = sol1.t[-1]
y1 = [sol1.y[0, -1]]
print(f'time={t1} y={y1} discontinuity_zero={discontinuity_zero(t1, y1)}')
sol2 = solve_ivp(ode_func2, t_span=[t1, tend], y0=y1, rtol=1e-8)
plt.plot(sol1.t, sol1.y[0,:])
plt.plot(sol2.t, sol2.y[0,:])
plt.show()
这将打印以下内容，其中不连续的时间精确到7位。
time=2.302584885712467 y=[10.000000000000002] discontinuity_zero=1.7763568394002505e-15