为了账号安全,请及时绑定邮箱和手机立即绑定

如何使用 scipy.integrate.odeint 正确实现具有力时间因变量的强制质量弹簧系统

如何使用 scipy.integrate.odeint 正确实现具有力时间因变量的强制质量弹簧系统

aluckdog 2022-09-13 19:47:53
对于这5种情况中的每一种,我都试图通过odeint函数(弹簧质量函数)进行数值积分,参数F随时间变化。但是,它显示的错误:第 244 行,in odeint int(布尔(第一))值错误:设置具有序列的数组元素。from scipy.integrate import odeint as odeimport matplotlib.pyplot as graphimport numpy as npdef spring(y,t,par):    m = par[0]    k = par[1]    c = par[2]    F = par[3]     ydot=[0,0]    ydot[0] = y[1]    F = np.array(F)    ydot[1] = F/m-c*y[1]/m-k*y[0]/m    return ydot#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,4m=10.0k=3553.0par = [m,k]h = 0.01cont = 0for case in cases:     x0 = case[2]    xdot0 = case[3]    y = [x0,xdot0]    par.append(case[4])    ti = case[0]    tf = case[1]    t = np.arange(ti, tf,h)    F = []    for time in t:        if cont == 3:            F.append(1000*np.sin(np.pi*time+np.pi/2))        elif (cont == 4) and (time >= 0.5):             F.append(1000)        else:            F.append(0)    cont = cont + 1    par.append(F) #F is a vector    Y = ode(spring,y,t,args=(par,))    Xn = Y[:,0]    Vn = Y[:,1]    graph.plot(t,Xn)    graph.show()    graph.plot(t,Vn)    graph.show()    graph.plot(Xn,Vn)    graph.show()
查看完整描述

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

并检查发生了什么。


查看完整回答
反对 回复 2022-09-13
?
MM们

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 个示例给出了混沌振荡的结果

//img1.sycdn.imooc.com//63206e2f0001d58606320207.jpg

查看完整回答
反对 回复 2022-09-13
?
慕丝7291255

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


查看完整回答
反对 回复 2022-09-13
  • 3 回答
  • 0 关注
  • 123 浏览
慕课专栏
更多

添加回答

举报

0/150
提交
取消
微信客服

购课补贴
联系客服咨询优惠详情

帮助反馈 APP下载

慕课网APP
您的移动学习伙伴

公众号

扫描二维码
关注慕课网微信公众号