如何使用scipy.integrate.odeint使用受力时间相关变量正确实现受力弹簧系统ode解算器

问题描述 投票:1回答:3

我试图针对这5种情况中的每一种,通过odeint函数,弹簧质量函数进行数值积分,参数F随时间变化。但是,它显示错误:

第244行,在odeint中int(bool(tfirst)))

ValueError:设置具有序列的数组元素。

from scipy.integrate import odeint as ode
import matplotlib.pyplot as graph
import numpy as np

def 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,4

m=10.0
k=3553.0
par = [m,k]
h = 0.01
cont = 0

for 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()

python dictionary scipy ode runge-kutta
3个回答
0
投票

在解决该问题后,尝试运行时代码中还会出现其他问题。但是为了精确起见,我只会回答您的问题。

问题出在这行代码:

ydot[1] = F/m-c*y[1]/m-k*y[0]/m

此处的F当前作为列表提供。在数学中,如果将矢量除以标量,则假定执行了逐元素除法。 Python列表不会复制此行为,但numpy数组会复制。因此,只需将列表转换为numpy数组,如下所示:

F = np.array(F)
ydot[1] = F/m-c*y[1]/m-k*y[0]/m

0
投票

问题出在spring(y,t,par)的返回值上>

odeint期望从ydot[0]返回2个单独的ydot[1]spring()值,但在ydot[0]中得到一个值,然后在len(F)中得到一个长度为ydot[1]的数组。

更改,

ydot[1] = F / m - c * y[1] / m - k * y[0] / m

to

ydot[1] = F[0] / m - c * y[1] / m - k * y[0] / m

并检查会发生什么。。


0
投票

将值数组F传递给派生函数,而没有任何将其条目与函数求值的时间t相关联的方法,这是完全错误的。请记住,参数yt是求解器在数值方法中计算下一阶段所需的状态向量和(单,标量)时间。最简单的方法是将F作为函数传递,以便可以直接计算其值。


0
投票

将值数组F传递给派生函数,而没有任何将其条目与函数求值的时间t相关联的方法,这是完全错误的。请记住,参数yt是求解器在数值方法中计算下一阶段所需的状态向量和(单,标量)时间。最简单的方法是将F作为函数传递,以便可以直接计算其值。

© www.soinside.com 2019 - 2024. All rights reserved.