3 回答

TA贡献1808条经验 获得超4个赞
问题在于来自spring(y,t,par)
odeint
期望从 中返回 2 个单值,而是在 中获取单个值,然后在 中获取一个长度数组。ydot[0]
ydot[1]
spring()
ydot[0]
len(F)
ydot[1]
改变
ydot[1] = F / m - c * y[1] / m - k * y[0] / m
自
ydot[1] = F[0] / m - c * y[1] / m - k * y[0] / m
并检查发生了什么。

TA贡献1886条经验 获得超2个赞
将值数组传递给导数函数,而没有任何方法将其条目与函数的计算时间相关联,这是完全错误的。请记住,参数 和 是状态向量和(单个标量)时间,求解器需要计算数值方法中的下一阶段。最简单的方法是作为函数传递,以便可以直接计算其值。FtytF
def spring(y,t,par):
m,k,c,F = par
return [ y[1], F(t)/m-c*y[1]/m-k*y[0]/m ]
#cases:[ti,tf,x0,xdot0,c]
A=[0.0,0.3,0.1,0.0,37.7]
B=[0.0,3.0,0.0,1.0,37.7]
C=[0.0,3.0,0.1,1.0,377]
D=[0.0,5.0,0.0,0.0,37.7]
E=[0.0,3.0,0.0,0.0,37.7]
cases = [A, B, C, D, E] #0,1,2,3,4
m=10.0
k=3553.0
h = 0.01
使用元组赋值来解压缩参数元组更具有诗意味,在一行中完成以前分布在多个参数上的内容,也使参数的顺序更加可见。
删除一些不必要的并发症,例如具有额外的计数机制,请使用该机制。此外,不要使用显式构造的参数列表,因为在可选参数中临时构造它要容易得多(如果有数百个系统构造的条目,这将有所不同)。enumerate(list)parargs
在常规设置中,可能是使用 中的方法生成的插值函数。在这里,将给定的方程转换为lambda表达式以定义相应的标量函数就足够了。Fscipy.interpolate
for cont,case in enumerate(cases):
ti,tf,x0,xdot0,c = case
y = [x0,xdot0]
t = np.arange(ti, tf, h)
F = lambda t: 0
if cont == 3:
F = lambda t: 1000*np.sin(np.pi*t+np.pi/2)
elif (cont == 4):
F = lambda t: 1000 if t >= 0.5 else 0
Y = ode(spring,y,t,args=([m,k,c,F],))
Xn,Vn = Y.T
graph.figure(figsize=(9,3))
graph.subplot(1,3,1); graph.plot(t,Xn)
graph.subplot(1,3,2); graph.plot(t,Vn)
graph.subplot(1,3,3); graph.plot(Xn,Vn)
graph.tight_layout(); graph.show()
其他一切都像以前一样进行,使用子图将不同的图形组合在一张图片中。例如,第 4 个示例给出了混沌振荡的结果

TA贡献1859条经验 获得超6个赞
在修复问题后,当您尝试运行时,代码中还会出现其他问题。但为了准确起见,我将回答你的问题。
问题出在以下代码行中:
ydot[1] = F/m-c*y[1]/m-k*y[0]/m
此处的 F 当前作为列表提供。在数学中,如果将向量除以标量,则假定您执行元素除法。Python 列表不会复制这种行为,但麻木数组会复制。因此,只需将您的列表转换为numpy数组,如下所示:
F = np.array(F)
ydot[1] = F/m-c*y[1]/m-k*y[0]/m
添加回答
举报