scipy 中具有时变系数的 ODE

问题描述

我正在评估一组具有时变系数的 ODE

def deriv(y,t,N,coefficients):
    S,I,R = y
    dSdt = coefficients['beta'](t) * S * I / N * -1
    dIdt = coefficients['beta'](t) * S * I / N - coefficients['gamma']* I
    dRdt = coefficients['gamma'] * I
    return dSdt,dIdt,dRdt

特别是，我在预先计算的数组中有“beta”值，其大小等于 int(max(t))。

coefficients = {'beta' : beta_f,'gamma':0.1}
def beta_f(t):
     return  mybetas.iloc[int(t)]

# Initial conditions vector
y0 = (S0,I0,R0)
# Integrate the SIR equations over the time grid,t.
ret = odeint(deriv,y0,args=(N,coefficients))

当我运行 odeint 时，它也会评估超出 max(t) 的值，从而引发 beta_f 中的索引超出范围错误。如何限制 odeint 的评估跨度？

解决方法

自 len(mybetas) == int(max(t)) 起，即使 t 的值不超过 max(t)，您也会得到越界错误。例如，mybetas.iloc[int(max(t))] 会给你越界错误，即使 int(max(t)) <= max(t) 对于 t 的正值也是如此。

但就您而言，odeint 确实会检查集成域之外的一些值。几周前我不得不处理一个与您类似的问题，以下关于 stackoverflow 的两个讨论真的很有帮助：

integrate.ode sets t0 values outside of my data range

Solve ODEs with discontinuous input/forcing data

第二个链接解释了为什么在 odeint 循环中的每个整数时间步上一个接一个地使用 for 求解 ODE，而不是让 odeint处理由 Beta 值的跳跃引起的导数的不连续性。

否则，如果这适合您的研究案例，您可以内插您的 beta，并让函数 beta_f 返回 beta 的内插值。当然，您必须将插值域稍微扩展到您的积分域之外，因为 odeint 可能想要评估某个 t 大于 max(t) 的导数。

numpy ode python scipy scipy