Scipy odeint 输出耦合 ODE 的错误解

Question

我正在尝试求解三个常微分方程组。此时我不确定我的方程是否设置错误或者是否是代码问题。

我正在尝试分析带有滑动点质量的摆，如下图所示：

我的自由体图和方程：

我的代码：

import scipy as sp
import numpy as np
import matplotlib.pyplot as plt

#pendulum
b1=2
g=9.81

#mass 1
m1=5
k1=10
b2=0
R1=5

#mass 2
m2=4
k12=9
b3=0
R2=10

#Initial conditions
r1_0=5
v_r_1_0=0

r2_0=12
v_r_2_0=0

theta_0=0
omega_0=0

S_0=[r1_0,v_r_1_0,r2_0,v_r_2_0,theta_0,omega_0]

#DiffEq function
def dSdt(t,S):
    r1,r1dot,r2,r2dot,theta,thetadot=S
    return [r1dot,
            (-k1*(r1-R1)-b2*r1dot-k12*(r1-R1-(r2-R2))+m1*g*np.cos(theta))/m1,
            r2dot,
            (-k12*(r1-R1-(r2-R2))-b3*r2dot+m2*g*np.cos(theta))/m2,
            thetadot,
            (-thetadot*b1-m1*g*(r1)*np.sin(theta)-m2*g*(r2)*np.sin(theta))/(m1*(r1**2)+m2*(r2**2))
           ]
#Solve
t=np.linspace(0,10,20)
sol=sp.integrate.odeint(dSdt,y0=S_0,t=t,tfirst=True)

r1_vals=sol.T[0]
r2_vals=sol.T[2]
theta_vals=sol.T[4]

plt.plot(t,sol.T[0])

这是 R1(t) 的绘图图片：

我已经为此苦苦挣扎了几个小时。如何获得不发散到无穷大的解决方案？ odeint 不能处理非线性方程吗？

Answer 1

两质量摆的正确运动方程如下图所示：

求解 ODE 系统的代码：

import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import display, clear_output,HTML


#pendulum
b1=100
g=9.81

#mass 1
m1=5
k1=10
b2=1

#mass 2
m2=4
k12=4
b3=0.5

#Initial conditions
r1_0=((m1+m2)*g/k1)+0
v_r_1_0=0

r2_0=((m2*g/k12)+r1_0)+5
v_r_2_0=0

theta_0=np.pi/3
omega_0=0

S_0=[r1_0,v_r_1_0,r2_0,v_r_2_0,theta_0,omega_0]

#DiffEq function
def dSdt(t,S):
    r1,r1dot,r2,r2dot,theta,thetadot=S
    return [r1dot,
            (-k1*(r1)-b2*r1dot-k12*(r1-r2)+m1*g*np.cos(theta))/m1,
            r2dot,
            (k12*(r1-r2)-b3*r2dot+m2*g*np.cos(theta))/m2,
            thetadot,
            (-thetadot*b1-m1*g*(r1)*np.sin(theta)-m2*g*(r2)*np.sin(theta))/(m1*(r1**2)+m2*(r2**2))
           ]
#Solve
t=np.linspace(0,60,400)
sol=sp.integrate.odeint(dSdt,y0=S_0,t=t,tfirst=True)
sol2=sp.integrate.solve_ivp(dSdt,t_span=(0,max(t)),y0=S_0,t_eval=t)

fig, axs = plt.subplots(2)
fig.suptitle('Mass Positions')
axs[0].set_ylim(0,30)
axs[0].plot(t, sol2.y[0],label='r1 [m]')
axs[0].plot(t, sol2.y[2],label='r2 [m]')
axs[1].set_xlabel('Time [sec]')
axs[0].set_ylabel('Distance from pivot [m]')
axs[1].set_ylabel('Theta [deg]')
axs[0].invert_yaxis()
axs[0].legend()

axs[1].plot(t,sol2.y[4]*(180/np.pi),label='theta [rad]')

剧情：