在函数中访问scipy.integrate.odeint中的当前时间步长

Question

有没有办法访问scipy.integrate.odeint中的当前时间步长

我试图解决一个颂歌系统，其中颂歌的形式取决于人口是否会枯竭。基本上，只要{}不低于阈值，我就从总体{}中获取。如果我需要采取此时间步的数量大于该阈值，我将采取所有x到该点，其余的从z开始

我正试图通过检查这个时间步长，然后在DEs中的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)

Answer 1

正如在另一个答案中指出的，每次调用odeint中的ODE函数时，时间可能不会严格增加，特别是对于僵硬的问题

处理ode函数中这种不连续性的最可靠的方法是使用事件来查找示例中(A * x - B * x * y) - 0.6的零的位置。对于不连续解，使用终端事件精确地在零处停止计算，然后更改ode函数。在solve_ivp中，可以使用events参数执行此操作。参见solve ivp documentation，特别是与炮弹轨迹相关的示例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

Answer 2

当然，你可以一起破解一些可能有用的东西。 您可以记录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)