使用scipy.integrate.odeint

Question

当我试图解决的颂歌基本上类似于这一个，但多一个弹簧和阻尼器==>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 

Answer 1

更改为：

from scipy.integrate import odeint 

和

x1,y1,x2,y2,x3,y3 = odeint(membre_droite, w0, time, args=(bump,)).T 

完整代码：

from scipy.integrate import odeint 
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 = 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()