python - 对值数组进行计算 - Python

Question

我正在尝试计算并绘制与时间相关的t运动微分方程解的振幅amp(请参阅rhs6)作为 wd 的函数，用于在 f_array 中找到力系数 f 的多个值。

到目前为止，我已经获得了针对 wd 绘制 amp 和 f 单个值的代码。结果是共振峰:

下面给出了针对 f 的一个值(生成上图)针对 wd 绘制 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()

现在我希望针对 wd 查找 amp 来获取 f 的多个值。我尝试通过在 f_array 上创建一个 for 循环语句来实现此目的。但是，我的方法不起作用，并且出现错误:

使用序列设置元素。

由于展示尝试很好，下面是我的尝试。

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 数组，并且您尝试将其中的元素设置为列表，因此使用序列设置元素。

要解决此问题，您可以创建一个二维 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 循环也不正确，因此我对其进行了编辑。该图看起来像