2016-02-23 24 views
0

我使用的是odeint,需要传递随时间变化的力以及我正在整合的位置和速度。这个力是一个已知的数据数组,所以它不需要解决,只是插入到等式中。如何将更多数据传入scipy.integrate.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) 

我知道有是Dfun参数odeint,我相信可以帮助我在这种情况下,但我无法找到如何使用它的许多信息。如果有人可以传递一些信息,那就太好了。或者有关如何在这种情况下使用interp1d的任何信息,或者简单地通过任何其他方式将f带入等式。

谢谢

(这是使用Python 2.7,如果不是由它的标题与SciPy的暗示。)

回答

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中所述。

+0

甜,谢谢你帮助很多。一切正常,但我仍然有一个问题。我只是意识到我仍然收到不正确的值,并意识到我在f上的单位是错误的。我需要将它们缩放一个因子或1000,并且现在接收到将所有t值全部保留时收到的错误。'x_new中的值高于插值范围'。有什么建议么? –

+0

因此,从我观察到的情况来看,t的值运行并最终达到1.86672,超过了提供的值,然后引发错误。唯一的问题是我无法理解如何防止这种情况/为什么会发生。 –

+0

另一个更新。看起来好像f在4448停止了。所以就我所见,odeint力量跳过它的极限,反过来,停止f short(不能插入)并且引起错误。仍然没有解决。 –

相关问题