Scipy odeint求解耦合常微分方程输出错误解，结果发散至无穷

问题概述 #

我正在求解一个包含三个耦合常微分方程（ODE）的滑动质点摆系统，目前无法确定是方程推导错误还是代码问题导致解发散至无穷。已完成受力图绘制与方程推导，运行代码后得到的r₁(t)曲线显示解持续发散，想解决以下问题：

1. 如何让系统解收敛？
2. odeint是否无法处理非线性方程？

运行代码 #

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

# 摆参数
b1=2
g=9.81

# 质点1参数
m1=5
k1=10
b2=0
R1=5

# 质点2参数
m2=4
k12=9
b3=0
R2=10

# 初始条件
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]

# 微分方程定义
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))
           ]
# 求解ODE
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])

故障原因与解决建议 #

核心结论 #

odeint完全支持非线性方程求解，发散是系统物理设定或方程推导的问题，而非求解器能力限制。以下是具体排查方向：

1. 初始状态与力平衡错误

你的初始条件下，系统处于非平衡状态，且无径向阻尼，导致持续加速：

• 初始时r1=5=R1（质点1弹簧原长位置），r2=12>R2=10（质点2偏离自身弹簧平衡位置），此时连接两质点的弹簧k12产生向外拉质点1的力，叠加重力径向分量（theta=0时重力沿r正方向），直接给质点1一个向外的初始加速度。
• 径向阻尼b2=0、b3=0，没有能量耗散机制，质点会持续加速远离悬挂点，最终导致解发散。

解决方式：

• 添加径向阻尼：将b2、b3设为非零值（如b2=2、b3=1），通过阻尼消耗能量，让系统收敛到平衡位置。
• 修正初始条件：计算theta=0时的径向平衡位置，以此作为初始值。平衡条件为：

  m1*g = k1*(r1_eq - R1) + k12*(r1_eq - r2_eq)
  m2*g = k12*(r2_eq - r1_eq)
  

  解得平衡位置后代入初始条件，避免初始状态的不平衡力。

2. 方程推导细节检查

重点核对径向力的符号与弹簧伸长量定义：

• 当前k12的力项为-k12*(r1-R1-(r2-R2))，展开后为-k12*((r1-r2)-(R1-R2))，若k12原长为R2-R1，则该表达式的符号逻辑正确，但需确保受力分析中弹簧力的方向与位移对应。
• 重力径向分量的符号需与坐标系定义一致：若r为悬挂点指向质点的距离，theta=0时重力沿r正方向，当前符号正确，但需结合平衡位置调整初始状态。

修正后测试代码示例 #

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

# 摆参数
b1=2
g=9.81

# 质点1参数
m1=5
k1=10
b2=2  # 添加径向阻尼
R1=5

# 质点2参数
m2=4
k12=9
b3=1  # 添加径向阻尼
R2=10

# 计算theta=0时的径向平衡位置
r2_eq_offset = (m2*g)/k12
r1_eq = R1 + (m1*g + m2*g)/k1
r2_eq = r1_eq + r2_eq_offset

# 初始条件（平衡位置+小扰动）
r1_0=r1_eq
v_r_1_0=0

r2_0=r2_eq
v_r_2_0=0

theta_0=0.1
omega_0=0

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

# 微分方程定义
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))
           ]
# 求解ODE（增加时间点数量使曲线平滑）
t=np.linspace(0,10,1000)
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.figure(figsize=(12,8))
plt.subplot(311)
plt.plot(t,r1_vals,label='r1(t)')
plt.ylabel('r1')
plt.legend()
plt.subplot(312)
plt.plot(t,r2_vals,label='r2(t)')
plt.ylabel('r2')
plt.legend()
plt.subplot(313)
plt.plot(t,theta_vals,label='theta(t)')
plt.ylabel('theta')
plt.xlabel('t')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者remusconnor