執行計算值數組 - 的Python

Question

我試圖計算和繪製的時間將溶液的振幅amp運動的依賴t微分方程（見rhs6）作爲wd的力系數f的多個值的函數發現於f_array。

到目前爲止，我已經得到了代碼amp針對wd繪製的單個值f。其結果是一個共振峯：

針對wd用於f一個值（其產生上述圖像）繪製amp的代碼如下。

from pylab import * 
from scipy.integrate import odeint 
from numpy import * 
import math 

#Parameters. 
k = 2.0 
m = 1.0 
w0 = (k/m)**(1/2) 
alpha = 0.2 
l = alpha/(2*m) 
f = 1.0 
wd = w0 + 0.025 
beta = 0.2 
t_fixed = 2.0 

#Arrays. 
t = linspace(0.0, 400.0, 400.0) 
wd_array = linspace(w0-1.0, w0+1, 400.0) 
f_array = linspace(10.0, 100.0, 3.0) 

#Time step. 
init_x = 0.0 
init_v = 0.0 
dx = 15.0 
dv = 0.0 
init_cond = [init_x,init_v] 
init_cond2 = [init_x + dx,init_v + dv] 

def rhs6(c,t,wd): 
    c0dot = c[1] 
    c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t) 
    return [c0dot, c1dot] 

amp_array=[] 
for wd in wd_array: 
    res = odeint(rhs6, init_cond, t, args=(wd,)) 
    amp = max(res[:,0]) 
    amp_array.append(amp) 

plot(wd_array, amp_array) 
xlabel('Driving frequency, wd') 
ylabel('Ampltiude, amp')  
show() 

現在我希望能夠找到amp對wd爲f多個值。我試圖通過在f_array上做出for循環聲明來做到這一點。然而，我的方法不起作用，我得到的錯誤：

setting an element with a sequence。

因爲這是一個很好的嘗試，下面是我的。

from pylab import * 
from scipy.integrate import odeint 
from numpy import * 
import math 

#Parameters. 
k = 2.0 
m = 1.0 
w0 = (k/m)**(1/2) 
alpha = 0.2 
l = alpha/(2*m) 
f = 1.0 
wd = w0 + 0.025 
beta = 0.2 
t_fixed = 2.0 

#Arrays. 
t = linspace(0.0, 200.0, 200.0) 
wd_array = linspace(w0-1.0, w0+1, 200.0) 
f_array = linspace(10.0, 200.0, 3.0) 

#Time step. 
init_x = 0.0 
init_v = 0.0 
dx = 15.0 
dv = 0.0 
init_cond = [init_x,init_v] 
init_cond2 = [init_x + dx,init_v + dv] 

def rhs6(c,t,wd,f): 
    c0dot = c[1] 
    c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t) 
    return [c0dot, c1dot] 

full_array = zeros(len(f_array))  
for index,force in enumerate(f_array): 
    amp_list = [] 
    for wd in wd_array: 
     res = odeint(rhs6, init_cond, t, args=(wd,force)) 
     amp = max(res[:,0]) 
     amp_list.append(amp) 
     print(res) 
    amp_array = array(amp_list) 
    full_array[index] = amp_array 

for f in full_array: 
    plot(wd, amp) 

show() 

任何想法？

Answer 1

你的問題是，full_array是一個numpy數組，你試圖設置一個元素作爲列表，因此setting an element with a sequence。

爲了解決這個可以改爲創建一個二維numpy的陣列，然後設置每一行中是amp_array如下

from pylab import * 
from scipy.integrate import odeint 
from numpy import * 
import math 

#Parameters. 
k = 2.0 
m = 1.0 
w0 = (k/m)**(1/2) 
alpha = 0.2 
l = alpha/(2*m) 
f = 1.0 
wd = w0 + 0.025 
beta = 0.2 
t_fixed = 2.0 

#Arrays. 
t = linspace(0.0, 200.0, 200.0) 
wd_array = linspace(w0-1.0, w0+1, 200.0) 
f_array = linspace(10.0, 200.0, 3.0) 

#Time step. 
init_x = 0.0 
init_v = 0.0 
dx = 15.0 
dv = 0.0 
init_cond = [init_x,init_v] 
init_cond2 = [init_x + dx,init_v + dv] 

def rhs6(c,t,wd,f): 
    c0dot = c[1] 
    c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t) 
    return [c0dot, c1dot] 

full_array = zeros((len(f_array),len(wd_array))) 
for index,force in enumerate(f_array): 
    amp_list = [] 
    for wd in wd_array: 
     res = odeint(rhs6, init_cond, t, args=(wd,force)) 
     amp = max(res[:,0]) 
     amp_list.append(amp) 
     # print(res) 
    amp_array = array(amp_list) 
    full_array[index,:] = amp_array 

for f in full_array: 
    plot(wd_array, f) 

show() 

或者，也可以按照gboffi's的建議，使用任一np.vstack或np.hstack創建你的full_array。

您的for循環結束繪製數組也是不正確的，所以我編輯了它。該圖看起來像

Answer 2

我碰到你的代碼和回溯是

--------------------------------------------------------------------------- 
ValueError        Traceback (most recent call last) 
<ipython-input-25-84f8768cb39a> in <module>() 
     7     print(res) 
     8   amp_array = array(amp_list) 
----> 9   full_array[index] = amp_array 
    10 
ValueError: setting an array element with a sequence. 

你爲什麼不建立full_array增量，使用numpy.vstack或numpy.hstack？

看來您的大部分時間都花在解決ODE上，所以 構造響應數組中的效率問題應該沒有關係。