双摆运动的GGB模拟

理论推导参考该作者：https://www.jianshu.com/p/844287eb0b34

最终运用分析力学的方法，推导出如下的动力学方程

首先，在GGB中输入初始参数，如下：

g=9.8

悬挂点A=(0,0)

摆长l_1=10

摆长l_2=10

质量m_1=5

质量m_2=5

终止时间t_f=100

球1的初始角度=pi/6

球2的初始角度=pi/2

球1的初始角速度=0

球2的初始角速度=0

动画播放速度speed=slider(0.1,1,0.1)

滑动条a=slider(0,1)

设置a的播放速度为speed

输入函数，如下：

xa(θ_{1},θ_{2},ω_{1})=l_{1} ω_{1}^(2) sin(θ_{1}-θ_{2})-g sin(θ_{2})

xb(θ_{1},θ_{2},ω_{2})=-(l_{2} m_{2} ω_{2}^(2) sin(θ_{1}-θ_{2})+(m_{1}+m_{2}) g sin(θ_{1}))

θ_{1}'(t,θ_{1},θ_{2},ω_{1},ω_{2})=ω_{1}

θ_{2}'(t,θ_{1},θ_{2},ω_{1},ω_{2})=ω_{2}

ω_{1}'(t,θ_{1},θ_{2},ω_{1},ω_{2})=((xb(θ_{1},θ_{2},ω_{2})-m_{2} cos(θ_{1}-θ_{2}) xa(θ_{1},θ_{2},ω_{1}))/(l_{1} (m_{1}+m_{2})-m_{2} l_{1} cos^(2)(θ_{1}-θ_{2})))

ω_{2}'(t,θ_{1},θ_{2},ω_{1},ω_{2})=((xa(θ_{1},θ_{2},ω_{1})-l_{1} cos(θ_{1}-θ_{2})*((xb(θ_{1},θ_{2},ω_{2})-m_{2} cos(θ_{1}-θ_{2}) xa(θ_{1},θ_{2},ω_{1}))/(l_{1} (m_{1}+m_{2})-m_{2} l_{1} cos^(2)(θ_{1}-θ_{2}))))/(l_{2}))

这里为了偷懒，我并没有对方程进行化简，有兴趣的可以自行化简，消去xa与xb以及合并同类项。

接着，数值解微分方程，如下：

解常微分方程组({θ_{1}',θ_{2}',ω_{1}',ω_{2}'},0,{θ_{01},θ_{02},ω_{01},ω_{02}},t_{f})

对解出的第一和第二个方程进行描点，例如：

B=描点(numericalIntegral1,a)

C=描点(numericalIntegral2,a)

时间t=x(B)

球P_{1}=(l_{1} sin(y(B)),-l_{1} cos(y(B)))

球P_{2}=P_{1}+(l_{2} sin(y(C)),-l_{2} cos(y(C)))

线f=线段(A,P_{1})

线h=线段(P_{1},P_{2})

大致模型制作完成，剩下的美观因素，自行优化即可。

成品如下：