通过 odeint（或 solve_ivp）计算系统对时变输入的响应

Posted 2023-03-12

【中文标题】通过 odeint（或 solve_ivp）计算系统对时变输入的响应

【英文标题】：Compute system response to time-varying input by odeint (or solve_ivp)

【发布时间】：2018-06-18 23:48:20

【问题描述】：

我正在尝试解决一个简单的 ODE 以可视化时间响应，这适用于使用 SciPy 中新的 solve_ivp 集成 API 的恒定输入条件。例如：

def dN1_dt_simple(t, N1):
    return -100 * N1

sol = solve_ivp(fun=dN1_dt_simple, t_span=[0, 100e-3], y0=[N0,])

但是，我想知道是否可以绘制对随时间变化的输入的响应？例如，与其将y0 固定为N0，我可以找到对简单正弦曲线的响应吗？

是否有兼容的方式将时变输入条件传递到 API？

【问题讨论】：

你的意思是像def dN1_dt(t, N1): return -100*N1 + np.sin(t)这样的东西吗？

是的，我想是的。这看起来简单得令人尴尬！谢谢。但是，是否有可能只获得对第一个周期的响应，这样我就可以看到系统是如何衰减的？

您可以使用t_span = [0, T]，其中T 是时变输入的周期，例如2*pi 代表 sint(t)。

感谢您的建议，尽管我对长期衰减感兴趣。即，如果我用一个周期激发一个系统，那么在该周期结束后的一段时间内会发生什么。这个 API 可以做到吗？

接受的答案是一种选择。另一种选择是分两步求解系统：在一个周期后获取输出，然后再次求解系统，使用第一部分的最终状态作为系统的初始条件，关闭周期性输入。

【参考方案1】：

正是出于这个原因，您传递给solve_ivp 的函数在其签名中具有t。您可以随心所欲地使用它¹。例如，要获得平滑的一次性脉冲，您可以：

from numpy import pi, cos

def fun(t,N1):
    input = 1-cos(t) if 0<t<2*pi else 0
    return -100*N1 + input

sol = solve_ivp(fun=fun, t_span=[0,20], y0=[N0])

请注意，y0 不是您使用该术语时的输入，而是初始条件。它仅在一个时间点定义且有意义 - 您开始集成/模拟的位置。 使用 ODE，您通常将外部输入建模为力或类似的（影响系统的时间导数，如上例所示），而不是直接更改状态。 使用这种方法并在您的可兴奋系统的上下文中，N0 已经是一些外部输入的结果。

---

¹只要它对各个积分器的需求足够平滑，通常是连续可微分的 (C¹)。如果你想实现一个步骤，最好使用非常尖锐的sigmoid。