首页 文章

如何将更多数据传递到scipy.integrate.odeint

提问于
浏览
0

我正在使用odeint并且需要传递随时间变化的力以及我正在积分的位置和速度 . 力是一个已知的数据阵列,所以它不需要解决,只是插入方程式 .

这是代码:

def dr_dt(y, t):

    RHO = 1225.0          
    C_D = 0.75             
    A = 6.25e-4            
    G = 9.81               
    M_O = 100.0            
    M_P = 10.8             
    M_F = M_O - M_P        
    T = 1.86               

    M_E = (M_O - M_F) / T

    dy0 = y[1]
    dy1 = (f / (M_O - M_E * t)) - ((1.0 * RHO * C_D * A * y[1]**2)  / (2.0 * (M_O - M_E * t))) - G
    return dy0, dy1

t = np.array([0.031, 0.092, 0.139, 0.192, 0.209, 0.231, 0.248, 0.292, 0.370, 0.475, 0.671, 0.702,
          0.723, 0.850, 1.063, 1.211, 1.242, 1.303, 1.468, 1.656, 1.821, 1.834, 1.847, 1.860])
f = np.array([0.946, 4.826, 9.936, 14.090, 11.446, 7.381, 6.151, 5.489, 4.921, 4.448, 4.258, 
          4.542, 4.164, 4.448, 4.353, 4.353, 4.069, 4.258, 4.353, 4.448, 4.448, 2.933, 1.325, 0.000])

r_o = 0.0
v_o = 0.0
y = odeint(dr_dt, [r_o, v_o], t)

我知道odeint中有Dfun参数,我认为在这种情况下可以帮助我,但我找不到有关如何使用它的更多信息 . 如果有人可以传递一些信息,那就太好了 . 或者有关如何在这种情况下使用interp1d的任何信息,或者只是通过任何其他方式来获取方程式 .

谢谢

(这是使用python 2.7,如果 Headers 中没有暗示scipy . )

1 回答

  • 1

    您可以尝试从强制数据创建插值对象,并将其作为参数发送到 dr_dt

    import numpy as np
    from scipy.integrate import odeint
    from scipy.interpolate import interp1d
    
    def dr_dt(y, t, fint):
    
        RHO = 1225.0          
        C_D = 0.75             
        A = 6.25e-4            
        G = 9.81               
        M_O = 100.0            
        M_P = 10.8             
        M_F = M_O - M_P        
        T = 1.86               
    
        M_E = (M_O - M_F) / T
    
        dy0 = y[1]
        dy1 = (fint(t) / (M_O - M_E * t)) - ((1.0 * RHO * C_D * A * y[1]**2)  / (2.0 * (M_O - M_E * t))) - G
        return dy0, dy1
    
    t = np.array([0.031, 0.092, 0.139, 0.192, 0.209, 0.231, 0.248, 0.292, 0.370, 0.475, 0.671, 0.702,
              0.723, 0.850, 1.063, 1.211, 1.242, 1.303, 1.468, 1.656, 1.821, 1.834, 1.847, 1.860])
    f = np.array([0.946, 4.826, 9.936, 14.090, 11.446, 7.381, 6.151, 5.489, 4.921, 4.448, 4.258, 
              4.542, 4.164, 4.448, 4.353, 4.353, 4.069, 4.258, 4.353, 4.448, 4.448, 2.933, 1.325, 0.000])
    
    r_o = 0.0
    v_o = 0.0
    
    fint = interp1d(t, f)
    y = odeint(dr_dt, [r_o, v_o], t[:-1], args=(fint,))
    

    (我发现我不得不省略最后一个时间点,因为否则它试图插入超出原始数据的范围......我不知道这对你来说是否重要) .

    编辑:如果这对您很重要,那么还有其他 ode 函数可以整合您的微分方程,而不会超出系列中的最后一个时间点,如this question中所述 .

相关问题