我只是在玩 python 的科学图书馆，特别是这里的例子:http://faculty1.coloradocollege.edu/~sburns/toolbox/ODE_II.html
我通过定义“力函数”对其进行了修改，并尝试制作一个简单的(但非物理示例)，其中力取决于对象的 x 和 y 位置。另外我想用 quiever 绘制力场。问题是我不明白如何使力函数正确地依赖于对象位置的 x 和 y 分量(函数)，尤其是我得到下面引用的错误。
from pylab import *
from scipy.integrate import odeint
import numpy as np
## set initial conditions and parameters
g = 9.81 # acceleration due to gravity
th = 91 # set launch angle
th = th * pi/180. # convert launch angle to radians
v0 = 100 # set speed
x0=0 # specify initial conditions
y0=0
vx0 = v0*sin(th)
vy0 = v0*cos(th)
## define force Funktion
def F_func(x):
F = zeros(2)
F[0] = x[1] #
F[1] = x[0] #
return F
## define function to compute f(X,t)
def f_func(state,time):
f = zeros(4) # create array to hold f vector
f[0] = state[2] # f[0] = x component of velocity
f[1] = state[3] # f[1] = x component of velocity
f[2] = F_func(state[:2])[0] # f[2] = acceleration in x direction
f[3] = F_func(state[:2])[1] # f[3] = acceleration in y direction
return f
## set initial state vector and time array
X0 = [ x0, y0, vx0, vy0] # set initial state of the system
t0 = 0.
tf = 10
tau = 0.1
#tf = input("Enter final time: ")
#tau = input("Enter time step: ")
# create time array starting at t0, ending at tf with a spacing tau
t = arange(t0,tf,tau)
## solve ODE using odeint
X = odeint(f_func,X0,t) # returns an 2-dimensional array with the
# first index specifying the time and the
# second index specifying the component of
# the state vector
print X
# putting ':' as an index specifies all of the elements for
# that index so x, y, vx, and vy are arrays at times specified
# in the time array
x = X[:,0]
y = X[:,1]
vx = X[:,2]
vy = X[:,3]
## plot the trajectory
fig = figure()
ax = fig.add_subplot(1,1,1)
## Enlarge Limits by en Percent
en = 0.05
#xMin,xMax = 0,10
xMin,xMax = min(x),max(x)
yMin,yMax = min(y),max(y)
xMin,xMax = xMin - (xMax-xMin)*en,xMax + (xMax-xMin)*en
yMin,yMax = yMin - (yMax-yMin)*en,yMax + (yMax-yMin)*en
#plot(x,y,[xMin,xMax],[yMin,yMax])
#plot(x,y,[0,10],[0,10])
ax.plot(x,y)
ax.axis('tight')
xlim([xMin,xMax])
ylim([yMin,yMax])
xG,yG = meshgrid(linspace(xMin,xMax,10),linspace(yMin,yMax,5))
ax.quiver(xG,yG,F_func(zip(xG,yG))[0],F_func(zip(xG,yG))[1],pivot='middle',minshaft=5,minlength=1,alpha=0.1)
xlabel('x')
ylabel('y')
show()
使用此代码我得到以下错误:
ValueError Traceback (most recent call last)
/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
176 else:
177 filename = fname
--> 178 __builtin__.execfile(filename, *where)
/home/myuser/python/test.py in <module>()
87 xG,yG = meshgrid(linspace(xMin,xMax,10),linspace(yMin,yMax,5))
88
---> 89 ax.quiver(xG,yG,F_func(zip(xG,yG))[0],F_func(zip(xG,yG))[1],pivot='middle',minshaft=5,minlength=1,alpha=0.1)
90
91
/home/myuser/python/test.py in F_func(x)
20 def F_func(x):
21 F = zeros(2)
---> 22 F[0] = x[1] #
23 F[1] = x[0] #
24 return F
ValueError: setting an array element with a sequence.
有人可以解释一下这个问题以及如何解决吗？
请您参考如下方法:
尝试使用
def F_func(x):
F1 = x[1]
F2 = x[0]
return array([F1, F2])
而不是你的 F_func。 并将 quiver 的调用替换为
ax.quiver(xG,yG,F_func(array([xG,yG]))[0],F_func(array([xG,yG]))[1],pivot='middle',minshaft=5,minlength=1,alpha=0.1)
您可以通过在 F_func 的第一行放置一个 print x 来找出 numpy 给您错误的原因。
F 只是一个双元素数组，您正试图将数组分配给它的两个元素中的每一个(因此出现错误“setting an array element with a sequence”)。
发生这种情况是因为 zip(xG, yG) 做了一些您可能不希望它做的事情。 (尝试在 ipython 中使用它)。
这是我在运行您的代码时得到的 ODE 解的结果图:
