二阶微分方程组 - 编程之家

如何解决二阶微分方程组

我需要可视化连接到振荡器的摆锤。所以我有二阶微分方程组：

(m1 + m2)*x1'' + m2*l*(x2''*cos(x2) - x2'*x2'*sin(x2)) + c*q1 = 0

x1''*cos(x2) - x1'*x2'*sin(x2) + l*x2'' + x1'*x2'*x2'*sin(x2) + g*x2'*sin(x2) = 0

m1 - mass of oscillator
m2 - mass of pendulum

l - length of pendulum
c - spring rate
g = 9.8

x1 - start position of oscillator
x2 - angle of pendulum

我试图为python找到库来求解微分方程。我发现臭。因此，我需要制作y1=x1'和y2 = x2'，然后为该库表示y1'和y2'。

from math import cos,sin,pi

def vectorfield(w,p):
    q1,y1,q2,y2 = w
    m1,m2,c,l,g = p

    f = [y1,(-m2*l*(c*q1*cos(q2)/(m1+m2) - m2*l*y2*y2*sin(q2)*cos(q2)/(m1+m2) + y1*y2*sin(q2)
          - y1*y2*y2*sin(q2) - g*y2*sin(q2)) / (l - m2*l*cos(q2)*cos(q2)/(m1+m2))*cos(q2)
          + m2*l*y2*y2*sin(q2) - c*q1) / (m1 + m2),y2,(c*q1*cos(q2)/(m1+m2) - m2*l*y2*y2*sin(q2)*cos(q2)/(m1+m2) + y1*y2*sin(q2)
          - y1*y2*y2*sin(q2) - g*y2*sin(q2)) / (l - m2*l*cos(q2)*cos(q2)/(m1+m2))
         ]
    return f

def customCos(val):
    return cos(pi*val/180)

def customSin(val):
    return sin(pi*val/180)

然后将此函数放入odeint

m1 = 1.0
m2 = 0.5
# Spring constants
c = 200.0
l = 0.1
g = 9.8

# Initial conditions
# x1 and x2 are the initial displacements; y1 and y2 are the initial veLocities
q1 = 1.0
y1 = 0.0
q2 = 0.5
y2 = 0.0

stoptime = 10.0
numpoints = 20

t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

# Pack up the parameters and initial conditions:
p = [m1,g]
w0 = [q1,y2]

# Call the ODE solver.
wsol = odeint(vectorfield,w0,t,args=(p,))

但是结果看起来确实是错误的。那么我的计算是错误的还是我以错误的方式使用了这个库？