使用odeint的数组中具有时间相关常数的微分方程组

Question

假设我有一个微分方程系统，我想用odeint解决它。一些系统的常量是时间依赖的，我将它们的值存储在数组中（a，b，c和d的形状（8000，））。我希望系统在每个时间步使用这些常量的不同值。请参阅简化的代码示例：

t = np.linspace(0,100,8000)

a = np.array([0, 5, 6, 12, 1.254, ..., 0.145])     # shape (8000, )
b = np.array([1.45, 5.9125, 1.367, ..., 3.1458])
c = np.array([0.124, 0.258, 0.369, ..., 0.147])
d = np.array([7.145, 5.123, 6.321, ..., 0.125])

def system(k,t):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)

    return [vcx_i_dot, vcy_i_dot, psi_i_dot wz_i_dot]

k0 = [0.1257, 0, 0, 0]

k = odeint(system, k0, t)

vcx_i = k[:,0]
vcy_i = k[:,1]
psi_i = k[:,2]
wz_i = k[:,3]

psi_i = [system(t_i, k_i)[2] for k_i, t_i in zip(k,t)]
wz_i = [system(t_i, k_i)[3] for k_i, t_i in zip(k,t)]

到目前为止我找到的最相关的解决方案是：Solving a system of odes (with changing constant!) using scipy.integrate.odeint?但是因为我只有数组中变量的值而不是依赖于时间的变量方程（例如a = f（t）），所以我尝试了插值在我的数组中的值之间，如下所示：ODEINT with multiple parameters (time-dependent)我设法使代码运行没有错误，但总时间急剧增加，系统解决的结果完全错误。我尝试了我在这里找到的任何可能类型的插值：https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html但仍然是错误的结果。这意味着我的插值不是最好的，或者我的数组中的点（8000个值）太多而无法在它们之间进行插值并正确地解决系统问题。解决这样一个系统的正确方法是什么？我是python的新手，我在Ubuntu 16.04 LTS上使用python 2.7.12。先感谢您。

Answer 1

插值器通常非常快，所以你的函数可能还有其他内容。但是，您可以尝试不同的插值器（如InterpolatedUnivariateSpline），或减少插值节点以提高速度。但我会瞄准你的整合。

最近，ode和odeint已经被其他更灵活的功能所取代（参见here）

我会先用显式方法而不是隐式方法（solve_ivp的默认值是runge kutta，而odeint的默认值是LSODA）：

interp = scipy.interpolate.interp1d(t,(a,b,c,d))

def system(t,k):
    vcx_i = k[0]
    vcy_i = k[1]
    psi_i = k[2]
    wz_i = k[3]
    a, b, c, d = interp(t)

    vcx_i_dot = a*np.cos(psi_i)-b*np.sin(psi_i)
    vcy_i_dot = b*np.cos(psi_i)+a*np.sin(psi_i)
    psi_i_dot = wz_i
    wz_i_dot = c*vcx_i-a*np.sin(psi_i)-d*np.sin(psi_i)-b*np.cos(psi_i)
    return [vcx_i_dot, vcy_i_dot, psi_i_dot, wz_i_dot]

k0 = [0.1257, 0, 0, 0]
steps = 1
method = 'RK23'
atol = 1e-3
s = solve_ivp(dydt, (0, 100), k0, method=method, t_eval=t, atol=atol, vectorized=True)

你可以增加/减少atol或改变方法。