我在 python 中使用 odeint 来求解微分方程，为什么我的质量位置发散到无穷大？

Question

我已经获得了涉及三个弹簧和两个质量的某个问题的运动方程。每侧的两个弹簧是非线性的，而中间的弹簧是线性的。它们也没有质量。 See the picture for clarity

我也把它放在一个代码中，以便能够求解二阶微分方程，我的代码如下：

   import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
initial=[-5,0,5,0] # [x,xdot,x2,x2dot]
t = np.linspace(0,5,10000) # Creating a vector which will represent time. 10 seconds divided in to 10000 intervals

def func(initials,t): # Defining function which is used in the odeint for solving diff.eq
    m1=5                    #Mass of M_1
    m2=5                    #Mass of M_2
    k12=10                  #Spring constant of spring connected and inbetween M_1 and M_2
    k1=10                   #Spring constant of spring connected to M_1
    k2=10                   #Spring constant of spring connected to M_2
    c1=2                    #constant of C for eq related to M_1
    c2=2                    #constant of C for eq related to M_2
    L1=5                    #Length of spring connected to M_1
    L2=5                    #Length of spring connected to M_2
    x1=initials[0]         #Initial values as chosen in row 5
    x2=initials[2]           #Initial values as chosen in row 5
    x1dotdot=(-k12*(x1-x2)/m1)-(k1*x1/m1)+(c1*2*np.pi/(L1*m1)*np.sin(2*np.pi*x1/L1))
    x2dotdot=(k12/m2*(x1-x2))-(k2*x2/m2)+(c2*2*np.pi/(L2*m2)*np.sin(2*np.pi*x2/L2))
    return(initials[1],x1dotdot,initials[3],x2dotdot)
output = odeint(func,initial,t)


plt.plot(t,output[:,0],'g:',linewidth = 2, label = 'M_1')
plt.plot(t,output[:,2],'y:',linewidth = 2, label = 'M_2')

plt.legend()
plt.xlabel('Time')
plt.ylabel('Velocities of the masses M_1 and M_2')
plt.show()

这里的问题是，当我选择通过返回首字母 (1) 和首字母 (3) 来绘制 velocities 时，它们以初始值 -5 和 5 而不是我给它们的 0 开头从头开始。当我绘制位置时，我期望 -5 和 5。

另一个问题也是职位。当我选择通过返回首字母（0）和首字母（2）来绘制它们时，它们会偏离infinity。

我不知道该怎么理解。

Answer 1

我认为您唯一的问题是混淆从odeint 返回的output 中的预期内容。它将按照定义的顺序返回状态向量（此处为：X = [x1, x1dot, x2, x2dot]），而不是状态向量的导数。

我在下面添加了一些额外的标签/cmets，希望能够澄清，并返回以下

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# State vector: X = [x1, x1dot, x2, x2dot]
initial=[-5,0,5,0] # [x,xdot,x2,x2dot]
t = np.linspace(0,5,10000) # Creating a vector which will represent time. 10 seconds divided in to 10000 intervals

def func(X,t): # Defining function which is used in the odeint for solving diff.eq
    m1=5                    #Mass of M_1
    m2=5                    #Mass of M_2
    k12=10                  #Spring constant of spring connected and inbetween M_1 and M_2
    k1=10                   #Spring constant of spring connected to M_1
    k2=10                   #Spring constant of spring connected to M_2
    c1=2                    #constant of C for eq related to M_1
    c2=2                    #constant of C for eq related to M_2
    L1=5                    #Length of spring connected to M_1
    L2=5                    #Length of spring connected to M_2
    x1=X[0]                 #Current value of x1
    x2=X[2]                 #Current value of x2
    x1dotdot=(-k12*(x1-x2)/m1)-(k1*x1/m1)+(c1*2*np.pi/(L1*m1)*np.sin(2*np.pi*x1/L1))
    x2dotdot=(k12/m2*(x1-x2))-(k2*x2/m2)+(c2*2*np.pi/(L2*m2)*np.sin(2*np.pi*x2/L2))

    # Need to return derivate of State vector X,
    # X    = [x1,    x1dot,    x2,    x2dot]
    # Xdot = [x1dot, x1dotdot, x2dot, x2dotdot]
    return(X[1], x1dotdot, X[3], x2dotdot)

output = odeint(func,initial,t)

fig, (ax1, ax2) = plt.subplots(2 ,1)

# Positions are in columns 0 and 2, X = [x1, x1dot, x2, x2dot]
ax1.plot(t,output[:,0],'g:',linewidth = 2, label = 'M_1')
ax1.plot(t,output[:,2],'y:',linewidth = 2, label = 'M_2')
ax1.set_xlabel('Time')
ax1.set_ylabel('Positions of the masses')
ax1.legend()

# Velocities are in columns 1 and 3, X = [x1, x1dot, x2, x2dot]
ax2.plot(t,output[:,1],'g:',linewidth = 2, label = 'M_1')
ax2.plot(t,output[:,3],'y:',linewidth = 2, label = 'M_2')
ax2.set_xlabel('Time')
ax2.set_ylabel('Velocities of the masses')
ax2.legend()

plt.show()