在函数内的 scipy.integrate.odeint 中访问当前时间步长 [英] Access current time step in scipy.integrate.odeint within the function

问题描述

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

Is there a way to access what the current time step is in scipy.integrate.odeint?

我正在尝试解决一个 ODE 系统，其中 ODE 的形式取决于种群是否会耗尽.基本上，我从人口 x 中获取，前提是 x 不会低于阈值.如果我需要采取这个时间步长的数量大于该阈值，我会将所有 x 带到那个点，其余的来自 z.

I am trying to solve a system of ODEs where the form of the ode depends on whether or not a population will be depleted. Basically I take from population x provided x doesn't go below a threshold. If the amount I need to take this timestep is greater than that threshold I will take all of x to that point and the rest from z.

我试图通过检查我将采取多少时间步长，然后在 DE 中的群体 x 和 z 之间进行分配来做到这一点.

I am trying to do this by checking how much I will take this time step, and then allocating between populations x and z in the DEs.

为此，我需要能够访问 ODE 求解器中的步长，以计算此时间步长将采取的措施.我正在使用 scipy.integrate.odeint - 有没有办法在定义 odes 的函数中访问时间步?

To do this I need to be able to access the step size within the ODE solver to calculate what will be taken this time step. I am using scipy.integrate.odeint - is there a way to access the time step within the function defining the odes?

或者，您可以访问求解器中最后一次的内容吗?我知道它不一定是下一个时间步骤，但如果这是我能做的最好的，它对我来说可能是一个足够好的近似值.或者还有其他我没有想到的选择吗?

Alternatively, can you access what the last time was in the solver? I know it won't necessarily be the next time step, but it's likely a good enough approximation for me if that is the best I can do. Or is there another option I've not thought of to do this?

下面的 MWE 不是我的方程组，而是我可以想出的来试图说明我在做什么.问题是在第一个时间步长，如果时间步长为 1，那么人口会过低，但由于时间步长很小，最初你可以从 x 中获取所有内容.

The below MWE is not my system of equations but what I could come up with to try to illustrate what I'm doing. The problem is that on the first time step, if the time step were 1 then the population will go too low, but since the timestep will be small, initially you can take all from 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)

推荐答案

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

As pointed out in the another answer, time may not be strictly increasing at each call to the ODE function in odeint, especially for stiff problems.

在 ode 函数中处理这种不连续性的最健壮的方法是使用一个事件来找到 (A * x - B * x * y) - 0.6 的零位置在你的例子中.对于不连续解，使用终止事件将计算精确地停在零处，然后更改 ode 函数.在 solve_ivp 中，您可以使用 events 参数执行此操作.请参阅解决 ivp 文档，特别是与炮弹轨迹相关的示例.odeint 不支持事件，并且 solve_ivp 有一个 LSODA 方法可用，它调用与 odeint 相同的 Fortran 库.

The most robust way to handle this kind of discontinuity in the ode function is to use an event to find the location of the zero of (A * x - B * x * y) - 0.6 in your example. For a discontinuous solution, use a terminal event to stop the computation precisely at the zero, and then change the ode function. In solve_ivp you can do this with the events parameter. See the solve ivp documentation and specifically the examples related to the cannonball trajectories. odeint does not support events, and solve_ivp has an LSODA method available that calls the same Fortran library as odeint.

这是一个简短的示例，但您可能需要在求解 sol2 之前另外检查 sol1 是否到达了终端事件.

Here is a short example, but you may want to additionally check that sol1 reached the terminal event before solving for sol2.

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 位数字.

This prints the following, where the time of the discontinuity is accurate to 7 digits.

time=2.302584885712467  y=[10.000000000000002]   discontinuity_zero=1.7763568394002505e-15

这篇关于在函数内的 scipy.integrate.odeint 中访问当前时间步长的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！