旋转圆盘上弹簧摆的Python仿真优化问询

问题背景 #

参考旋转圆盘摆的仿真场景，将悬挂摆线替换为弹簧，构建随圆盘旋转的弹簧-质量系统仿真。通过SymPy推导拉格朗日方程时，出现计算耗时过长的问题，同时不确定初始条件的合理性，寻求更高效的计算方式或问题简化思路。

当前实现代码 #

import sympy as sp
from IPython.display import display

R = sp.symbols('R')
omega = sp.symbols('omega')
t = sp.symbols('t')
phi = sp.Function('phi')(t)
theta = sp.Function('theta')(t)
s = sp.Function('s')(t)
L = sp.symbols('L')
m = sp.symbols('m')
k = sp.symbols('k')
g = sp.symbols('g')

x = R*sp.cos(omega*t)+(L+s)*(sp.sin(theta)*sp.cos(phi))
y = R*sp.sin(omega*t)+(L+s)*(sp.sin(theta)*sp.sin(phi))
z = -(L+s)*sp.cos(theta)

xs = sp.diff(x,t)
ys = sp.diff(y,t)
zs = sp.diff(z,t)

v = xs**2 + ys**2 + zs**2
Ekin = 0.5*m*v
Epot = g*(L+s)*sp.cos(theta)+0.5*k*s**2

# 注意：拉格朗日量应为动能减势能，此处代码存在错误
L = Ekin + Epot

#display(L)

ELTheta = sp.diff(sp.diff(L,sp.Derivative(theta,t)), t) + sp.diff(L,theta)
ELPhi = sp.diff(sp.diff(L,sp.Derivative(phi,t)), t) + sp.diff(L,phi)
ELs = sp.diff(sp.diff(L,sp.Derivative(s,t)), t) + sp.diff(L,s)

Eq1 = sp.Eq(ELTheta,0)
Eq2 = sp.Eq(ELPhi,0)
Eq3 = sp.Eq(ELs,0)

LGS = sp.solve((Eq1,Eq2,Eq3),(sp.Derivative(theta,t,2),sp.Derivative(phi,t,2),sp.Derivative(s,t,2)))

thetadd = sp.simplify(LGS[sp.Derivative(theta,t,2)])
phidd = sp.simplify(LGS[sp.Derivative(phi,t,2)])
sdd = sp.simplify(LGS[sp.Derivative(s,t,2)])

关键问题修正与优化方案 #

1. 拉格朗日量核心错误修正 #

拉格朗日量的定义为动能（T）减去势能（V），即 L = Ekin - Epot，而非代码中的Ekin + Epot。这个错误不仅会导致动力学方程完全失效，还会增加不必要的计算复杂度，修正后可大幅减少后续化简的计算量。

2. 手动预简化表达式，降低SymPy计算负担 #

直接让SymPy展开所有三角函数项会导致计算量爆炸，可先手动利用三角恒等式简化：

• 合并x、y方向速度项：xs² + ys²可利用cos²α + sin²α=1和sin(A-B)=sinAcosB - cosAsinB化简，避免冗余的三角函数展开；
• 对含三角函数的表达式先调用sp.trigsimp()化简后再求导，例如：

  v_simplified = sp.trigsimp(xs**2 + ys**2 + zs**2)
  Ekin = 0.5*m*v_simplified
  

3. 旋转坐标系简化问题 #

切换到随圆盘旋转的非惯性坐标系（旋转角速度为ω），此时圆盘静止，质量块的位置表达式可消除ωt项，只需引入离心力和科里奥利力作为惯性力，拉格朗日量的形式会更简洁，推导效率显著提升。

4. 避免直接求解方程组 #

SymPy的solve()处理复杂方程组时耗时较长，可手动整理欧拉-拉格朗日方程，将二阶导数单独移项求解：例如对θ的方程，将含θ̈的项留在左边，其余项移到右边，再直接除以系数得到θ̈，无需调用solve()求解整个方程组。

初始条件建议 #

• 平衡位置初始：先计算圆盘匀速旋转时的静态平衡状态，此时质量块的弹簧伸长量s0、摆角θ0满足合力为零，以此为初始条件（s=s0, θ=θ0, ṡ=0, θ̇=0），可验证仿真的稳定性；
• 小扰动初始：在平衡位置附近施加微小偏移（如s=s0+0.01, θ=θ0+0.01）或初始角速度，观察系统的振动响应；
• 零初始速度：初始时设置所有速度项为零，避免非物理的初始冲击，简化仿真的初始状态验证。

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