【发布时间】：2017-10-25 03:07:03
【问题描述】：
我正在尝试解决一个基本上类似于这个的颂歌系统，但多了一个弹簧和阻尼器 ==> http://scipy-cookbook.readthedocs.io/items/CoupledSpringMassSystem.html
不过，我有一个小问题，因为我要实现的参数之一是时间相关的。我的第一次尝试如下：
import scipy as sci
import numpy as np
import matplotlib.pyplot as plt
def bump(t):
if t <= (0.25 / 6.9):
return (0.075 * (1 - np.cos(np.pi * 8 * 6.9 * t)))
else:
return 0
def membre_droite(w, t, p):
x1,y1,x2,y2,x3,y3 = w
m1,m2,m3,k1,k2,k3,l1,l2,l3,c2,c3 = p
f = [y1,
(-k1 * (x1 - l1 - bump(t)) + k2 * (x2 - x1 - l2) + c2 * (y2 - y1)) / m1,
y2,
(-c2 * (y2 - y1) - k2 * (x2 - x1 - l2) + k3 * (x3 - x2 - l3) + c3 * (y3 - y2)) / m2,
y3,
(-c3 * (y3 - y2) - k3 * (x3 - x2 - l3)) / m3]
return f
# Initial values
x11 = 0.08
y11 = 0
x22 = 0.35
y22 = 0
x33 = 0.6
y33 = 0
# Parameters
m1 = 90
m2 = 4000
m3 = 105
k1 = 250000
k2 = 25000
k3 = 30000
l1 = 0.08
l2 = x22-x11
l3 = x33-x22
c2 = 2500
c3 = 850
# Initial paramaters regrouped + time array
time=np.linspace(0.0, 5, 1001)
w0 = [x11,y11,x22,y22,x33,y33]
p0 = [m1,m2,m3,k1,k2,k3,l1,l2,l3,c2,c3]
x1,y1,x2,y2,x3,y3 = sci.integrate.odeint(membre_droite, w0, time, args=(p0,)).T
plt.plot(time,x1,'b')
plt.plot(time,x2,'g')
plt.plot(time,x3,'r')
plt.plot(time,y2,'yellow')
plt.plot(time,y3,'black')
plt.xlabel('t')
plt.grid(True)
plt.legend((r'$x_1$', r'$x_2$', r'$x_3$', r'$y_2$', r'$y_3$'))
plt.show()
我得到的错误是：
if t <= (0.25 / 6.9):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
我找过类似的案例，发现了这个话题 ==> Solving a system of odes (with changing constant!) using scipy.integrate.odeint?
然后我尝试将我的代码调整为这种格式：
import scipy as sci
import numpy as np
import matplotlib.pyplot as plt
def bump(t):
if t <= (0.25 / 6.9):
return (0.075 * (1 - np.cos(np.pi * 8 * 6.9 * t)))
else:
return 0
def membre_droite(w, t, bump):
x1,y1,x2,y2,x3,y3 = w
f = [y1,
(-250000 * (x1 - x11 - bump(t)) + 25000 * (x2 - x1 - x22 + x11) + 2500 * (y2-y1)) / 90,
y2,
(-2500 * (y2 - y1) - 25000 * (x2 - x1 - x22 + x11) + 30000 * (x3 - x2 - x33 + x22) + 850 * (y3 - y2)) / 4000,
y3,
(-850 * (y3 - y2) - 30000 * (x3 - x2 - x33 + x22)) / 105]
return f
# Initial values
x11 = 0.08
y11 = 0
x22 = 0.35
y22 = 0
x33 = 0.6
y33 = 0
# Initial paramaters regrouped + time array
time = np.linspace(0.0, 5, 1001)
w0 = [x11,y11,x22,y22,x33,y33]
x1,y1,x2,y2,x3,y3 = sci.integrate.odeint(membre_droite, w0, time, args=(bump,)).T
plt.plot(time,x1,'b')
plt.plot(time,x2,'g')
plt.plot(time,x3,'r')
plt.plot(time,y2,'yellow')
plt.plot(time,y3,'black')
plt.xlabel('t')
plt.grid(True)
plt.legend((r'$x_1$', r'$x_2$', r'$x_3$', r'$y_2$', r'$y_3$'))
plt.show()
阅读上一个链接，它应该可以工作，但我收到另一个错误：
(-250000 * (x1 - x11 - bump(t)) + 25000 * (x2 - x1 - x22 + x11) + 2500 * (y2 - y1)) / 90,
TypeError: 'list' object is not callable
【问题讨论】：
-
我已经对您的代码进行了一些修改，例如：
from scipy.integrate import odeint和x1,y1,x2,y2,x3,y3 = odeint(membre_droite, w0, time, args=(bump,)).T，它工作正常，您使用的是什么版本的 python 和库。 -
我使用的是 3.6 版本。不过，你对图书馆意味着什么？ （我对这一切都比较陌生，所以我看不出我可以改变什么）
-
这是我的输出：i.imgur.com/u864MSx.png
-
好吧，实际上，当我像你说的那样切换到“from scipy.integrate import odeint”时，它起作用了。我不知道是什么原因造成的，但非常感谢。
-
我将把那个答案写成答案。