是否有可能让 scipy.integrate.odeint 与一对向量（或向量的向量）一起工作
Is it possible to get scipy.integrate.odeint to work with a pair of vectors (or vector of vectors)
使用 scipy.integrate.odeint 解决自由落体问题
进口
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline
常数和积分区间
G = 9.8 # m / sec / sec
# First 1 second of motion
t = np.arange (0, 1.1, 0.1)
速度矢量解
这里我们return一个加速度向量在velocity_vector()中，设置初始V0，积分1秒。我们期望Vx和Vy是常数，Vz以-G的斜率线性下降。
def velocity_vector (x, t, params):
# x = (Vx, Vy, Vz)
# Ordinary differential equation - velocity of an object in frictionless free-fall.
g = params
acceleration = np.array ([0, 0, -g])
return acceleration
v0 = np.array ([1, 2, 0])
soln = odeint (velocity_vector, v0, t, args = (G,))
fig = plt.figure (1, figsize = (8,8))
ax1 = fig.add_subplot(311)
ax1.plot(t, soln [:,0])
ax1.set_xlabel ('time')
ax1.set_ylabel ('Vx')
ax1 = fig.add_subplot(312)
ax1.plot(t, soln [:,1])
ax1.set_xlabel ('time')
ax1.set_ylabel ('Vy')
ax1 = fig.add_subplot(313)
ax1.plot(t, soln [:,2])
ax1.set_xlabel ('time')
ax1.set_ylabel ('Vz')
plt.show ()
一维位置和速度的矢量解
这里我们return一个motion()中z维度的加速度和速度向量，设置初始Z0和V0，积分1秒。我们期望 Xz 是二次方的，Vz 以 -G 的斜率线性下降。
def motion (x, t, params):
# x = (Sx, Vx)
# Ordinary differential equation - velocity of an object in frictionless free-fall.
g = params
acceleration = np.array ([-g * t, -g])
return acceleration
v0 = np.array ([5.0, 0])
soln = odeint (motion, v0, t, args = (G,))
fig,axes = plt.subplots(1, 2) # one row, two columns
fig.subplots_adjust(wspace=0.6)
axes[0].plot(t, soln [:,0])
axes[0].set_xlabel ('time')
axes[0].set_ylabel ('Sz')
axes[1].plot(t, soln [:,1])
axes[1].set_xlabel ('time')
axes[1].set_ylabel ('Vz')
plt.show ()
一维位置和速度的矢量解
让我们看看是否可以通过提供向量列表和return输入向量列表来结合这两种方法。
def position_and_velocity (x, t, params):
# x = (S, V) as vectors
# Ordinary differential equation - velocity of an object in frictionless free-fall.
g = params
acceleration = np.array ([-g * t, -g])
return acceleration
s = np.array ([0, 0, 5])
v = np.array ([1, 2, 0])
SV0 = np.array ([s, v])
soln = odeint (position_and_velocity, SV0, t, args = (G,))
#fig,axes = plt.subplots(1, 2) # one row, two columns
#fig.subplots_adjust(wspace=0.6)
#axes[0].plot(t, soln [:,0])
#axes[0].set_xlabel ('time')
#axes[0].set_ylabel ('Sz')
#axes[1].plot(t, soln [:,1])
#axes[1].set_xlabel ('time')
#axes[1].set_ylabel ('Vz')
#plt.show ()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-29-a12f336fc4fc> in <module>()
9 v = np.array ([1, 2, 0])
10 SV0 = np.array ([s, v])
---> 11 soln = odeint (position_and_velocity, SV0, t, args = (G,))
12
13 #fig,axes = plt.subplots(1, 2) # one row, two columns
/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/scipy/integrate/odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg)
146 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
147 full_output, rtol, atol, tcrit, h0, hmax, hmin,
--> 148 ixpr, mxstep, mxhnil, mxordn, mxords)
149 if output[-1] < 0:
150 print(_msgs[output[-1]])
ValueError: object too deep for desired array
似乎scipy.integrate.odeint()可以处理一个向量并使用向量求解方程，但不能像本例中那样处理向量的向量。有没有办法避免 return 在这里使用六个不同的函数，而不是两个？
odeint 仅处理一维数组。要在组合系统上使用 odeint，您必须将两个 3 维向量连接成一个 6 维向量。如果您尝试使用两个现有函数来计算两个 3 维系统方程的两个右侧，则必须创建一个接受 6 维向量的新函数，将其拆分为适当的 3 维子向量，调用两个现有函数，然后将结果作为 6 维向量连接到 return.
正如 Warren 提到的，答案是 odeint 确实需要一维数组。诀窍是设置传递给 odeint 的函数，以便它将作为第二个传入的一维数组转换为所需的向量形式 - 在本例中为 2 个 3-D 向量，以向量形式执行计算，然后重塑结果返回为一维数组。
答案演示了该技术，使用方便的 numpy reshape 函数。
def position_and_velocity (x, t, params):
# x = (S, V) as vectors
# Ordinary differential equation - velocity of an object in frictionless free-fall.
G = params
g = np.array ([0, 0, -G])
# convert the 6 element vector to 2 3 element vectors of displacement and velocity
# to use vector formulation of the math
s,v = x.reshape (2,3)
acceleration = np.array ([v * t, g])
# reshape the two vector results back into one for odeint
return np.reshape (acceleration, 6)
s = np.array ([0, 0, 5])
v = np.array ([40, 10, 0])
SV0 = np.array ([s, v])
# pass reshaped displacement and velocity vector to odeint
soln = odeint (position_and_velocity, np.reshape (SV0, 6), t, args = (G,))
使用 scipy.integrate.odeint 解决自由落体问题
进口
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline
常数和积分区间
G = 9.8 # m / sec / sec
# First 1 second of motion
t = np.arange (0, 1.1, 0.1)
速度矢量解
这里我们return一个加速度向量在velocity_vector()中，设置初始V0，积分1秒。我们期望Vx和Vy是常数，Vz以-G的斜率线性下降。
def velocity_vector (x, t, params):
# x = (Vx, Vy, Vz)
# Ordinary differential equation - velocity of an object in frictionless free-fall.
g = params
acceleration = np.array ([0, 0, -g])
return acceleration
v0 = np.array ([1, 2, 0])
soln = odeint (velocity_vector, v0, t, args = (G,))
fig = plt.figure (1, figsize = (8,8))
ax1 = fig.add_subplot(311)
ax1.plot(t, soln [:,0])
ax1.set_xlabel ('time')
ax1.set_ylabel ('Vx')
ax1 = fig.add_subplot(312)
ax1.plot(t, soln [:,1])
ax1.set_xlabel ('time')
ax1.set_ylabel ('Vy')
ax1 = fig.add_subplot(313)
ax1.plot(t, soln [:,2])
ax1.set_xlabel ('time')
ax1.set_ylabel ('Vz')
plt.show ()
一维位置和速度的矢量解
这里我们return一个motion()中z维度的加速度和速度向量，设置初始Z0和V0，积分1秒。我们期望 Xz 是二次方的，Vz 以 -G 的斜率线性下降。
def motion (x, t, params):
# x = (Sx, Vx)
# Ordinary differential equation - velocity of an object in frictionless free-fall.
g = params
acceleration = np.array ([-g * t, -g])
return acceleration
v0 = np.array ([5.0, 0])
soln = odeint (motion, v0, t, args = (G,))
fig,axes = plt.subplots(1, 2) # one row, two columns
fig.subplots_adjust(wspace=0.6)
axes[0].plot(t, soln [:,0])
axes[0].set_xlabel ('time')
axes[0].set_ylabel ('Sz')
axes[1].plot(t, soln [:,1])
axes[1].set_xlabel ('time')
axes[1].set_ylabel ('Vz')
plt.show ()
一维位置和速度的矢量解
让我们看看是否可以通过提供向量列表和return输入向量列表来结合这两种方法。
def position_and_velocity (x, t, params):
# x = (S, V) as vectors
# Ordinary differential equation - velocity of an object in frictionless free-fall.
g = params
acceleration = np.array ([-g * t, -g])
return acceleration
s = np.array ([0, 0, 5])
v = np.array ([1, 2, 0])
SV0 = np.array ([s, v])
soln = odeint (position_and_velocity, SV0, t, args = (G,))
#fig,axes = plt.subplots(1, 2) # one row, two columns
#fig.subplots_adjust(wspace=0.6)
#axes[0].plot(t, soln [:,0])
#axes[0].set_xlabel ('time')
#axes[0].set_ylabel ('Sz')
#axes[1].plot(t, soln [:,1])
#axes[1].set_xlabel ('time')
#axes[1].set_ylabel ('Vz')
#plt.show ()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-29-a12f336fc4fc> in <module>()
9 v = np.array ([1, 2, 0])
10 SV0 = np.array ([s, v])
---> 11 soln = odeint (position_and_velocity, SV0, t, args = (G,))
12
13 #fig,axes = plt.subplots(1, 2) # one row, two columns
/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/scipy/integrate/odepack.py in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, rtol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg)
146 output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
147 full_output, rtol, atol, tcrit, h0, hmax, hmin,
--> 148 ixpr, mxstep, mxhnil, mxordn, mxords)
149 if output[-1] < 0:
150 print(_msgs[output[-1]])
ValueError: object too deep for desired array
似乎scipy.integrate.odeint()可以处理一个向量并使用向量求解方程，但不能像本例中那样处理向量的向量。有没有办法避免 return 在这里使用六个不同的函数，而不是两个？
odeint 仅处理一维数组。要在组合系统上使用 odeint，您必须将两个 3 维向量连接成一个 6 维向量。如果您尝试使用两个现有函数来计算两个 3 维系统方程的两个右侧，则必须创建一个接受 6 维向量的新函数，将其拆分为适当的 3 维子向量，调用两个现有函数，然后将结果作为 6 维向量连接到 return.
正如 Warren 提到的，答案是 odeint 确实需要一维数组。诀窍是设置传递给 odeint 的函数，以便它将作为第二个传入的一维数组转换为所需的向量形式 - 在本例中为 2 个 3-D 向量，以向量形式执行计算，然后重塑结果返回为一维数组。
答案演示了该技术，使用方便的 numpy reshape 函数。
def position_and_velocity (x, t, params):
# x = (S, V) as vectors
# Ordinary differential equation - velocity of an object in frictionless free-fall.
G = params
g = np.array ([0, 0, -G])
# convert the 6 element vector to 2 3 element vectors of displacement and velocity
# to use vector formulation of the math
s,v = x.reshape (2,3)
acceleration = np.array ([v * t, g])
# reshape the two vector results back into one for odeint
return np.reshape (acceleration, 6)
s = np.array ([0, 0, 5])
v = np.array ([40, 10, 0])
SV0 = np.array ([s, v])
# pass reshaped displacement and velocity vector to odeint
soln = odeint (position_and_velocity, np.reshape (SV0, 6), t, args = (G,))