使用odeint在数组中具有与时间相关的常数的微分方程组 [英] System of differential equations with time dependent constants in arrays, using odeint

问题描述

比方说，我有一个微分方程组，我想用odeint求解.系统的某些常数与时间有关，我将它们的值存储在数组中(形状为(8000，)的a，b，c和d).我希望系统在每个时间步使用这些常量的不同值.请参见简化的代码示例:

Let's say that I have a system of differential equations and I want to solve it with odeint. Some of the system's constants are time depended and I have their values stored in arrays (a,b,c and d with shape (8000, ) ). I want the system to use a different value of these constants at each time step. See the simplified code example:

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)]

到目前为止，我发现的最相关的解决方案是:具有多个参数的ODEINT(依赖于时间)我设法使代码运行没有错误，但是总时间却急剧增加，解决的系统结果完全错误.我尝试了在这里找到的任何可能的插值类型: https://docs.scipy.org/doc/scipy/reference/generation/scipy.interpolate.interp1d.html ，但结果仍然错误.这意味着我的插值方法不是最好的，或者我在数组中的点(8000个值)太多，无法在它们之间进行插值并正确地求解系统.解决这样一个系统的正确方法是什么? 我是python的新手，并且在Ubuntu 16.04 LTS上使用python 2.7.12.先感谢您.

The most relevant solution I found so far is this: Solving a system of odes (with changing constant!) using scipy.integrate.odeint? But since I have only values of my variables in arrays and not the equation of the variables that depend on time (e.g. a=f(t)), I tried to aply an interpolation between the values in my arrays, as shown here: ODEINT with multiple parameters (time-dependent) I managed to make the code running without errors, but the total time increased dramatically and the results of the system solved are completely wrong. I tried any possible type of interpolation that I found here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html but still wrong outcome. That means that my interpolation isn't the best possible, or my points in the arrays (8000 values) are too much to interpolate between them and solve the system correctly. What is the right way to solve a system like this? I am new to python and I use python 2.7.12 on Ubuntu 16.04 LTS. Thank you in advance.

推荐答案

插值器通常非常快，因此您的函数中可能还有其他内容.但是，您可以尝试使用不同的插值器(例如InterpolatedUnivariateSpline)，或减少插值节点以提高速度.但是我将目标放在您的集成上.

The interpolators are usually very fast, so probably there is something else in your function. You can however, attempt different interpolators (like InterpolatedUnivariateSpline), or reducing the interpolation nodes to increase speed. But I would aim at your integration instead.

最近，ode和odeint已被其他更灵活的功能所取代(请参见

Lately, ode and odeint have been replaced by other, more flexible functions (see here)

我将从显式方法而不是隐式方法开始(solve_ivp的默认值为runge kutta，而odeint的默认值为LSODA):

I would start with and explicit method instead of an implicit one (default for solve_ivp is runge kutta, while default for odeint is 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或更改方法.

you can increase / reduce atol or change the method.

这篇关于使用odeint在数组中具有与时间相关的常数的微分方程组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！