如何将更多数据传递到 scipy.integrate.odeint
How to pass more data into scipy.integrate.odeint
我正在使用 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 的任何信息,或者任何其他将 f 带入方程式的方法。
谢谢
(这是使用 python 2.7,如果标题中没有暗示 scipy。)
您可以尝试从您的力数据创建一个插值对象并将其作为参数发送给 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 中所述。
我正在使用 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 的任何信息,或者任何其他将 f 带入方程式的方法。
谢谢
(这是使用 python 2.7,如果标题中没有暗示 scipy。)
您可以尝试从您的力数据创建一个插值对象并将其作为参数发送给 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 中所述。