数组中具有时间相关常数的微分方程组，使用 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或者改变方法。