双摆问题
所谓双摆,就是两个连在一起的摆。
设小球 m 1 m_1 m1和 m 2 m_2 m2的坐标分别为 ( x 1 , y 2 ) , ( x 2 , y 2 ) (x_1,y_2),(x_2,y_2) (x1,y2),(x2,y2),则
x 1 = l 1 sin θ 1 y 1 = − l 1 cos θ 1 x 2 = x 1 + l 2 sin θ 2 y 2 = y 1 − l 2 cos θ 2 x_1=l_1\sin\theta_1\\ y_1=-l_1\cos\theta_1\\ x_2=x_1+l_2\sin\theta_2\\ y_2=y_1-l_2\cos\theta_2 x1=l1sinθ1y1=−l1cosθ1x2=x1+l2sinθ2y2=y1−l2cosθ2
如图所示,设 O O O为坐标原点,如果这两个连在一起的摆在同一个平面内,则 m 1 m_1 m1和 m 2 m_2 m2点的坐标为
x 1 = l 1 sin θ 1 x 2 = x 1 + l 2 sin θ 2 y 1 = − l 1 cos θ 1 y 2 = y 1 − l 2 cos θ 2 \begin{aligned} &x_1=l_1\sin\theta_1&x_2=x_1+l_2\sin\theta_2\\ &y_1=-l_1\cos\theta_1&y_2=y_1-l_2\cos\theta_2 \end{aligned} x1=l1sinθ1y1=−l1cosθ1x2=x1+l2sinθ2y2=y1−l2cosθ2
接下来本来是要推公式的,考虑考虑到大家可能会有公式恐惧症,同时又喜欢看图,所以把公式挪到后面。
所以,只需知道角速度的微分方程,就可写出对应的代码,其方程如下
ω ˙ 1 = ( − m g l 1 sin θ 1 + μ λ ω 2 2 sin δ ) − μ cos δ ( − g l 1 sin θ 2 − ω 1 2 sin δ ) m − μ cos 2 δ ω ˙ 2 = − cos δ ( − m g l 1 sin θ 1 + μ λ ω 2 2 sin δ ) + m ( − g l 1 sin θ 2 − ω 1 2 sin δ ) m λ − μ λ cos 2 δ \dot\omega_1=\frac{(-m\frac{g}{l_1}\sin\theta_1+\mu\lambda\omega_2^2\sin\delta)-\mu\cos\delta(-\frac{g}{l_1}\sin\theta_2-\omega_1^2\sin\delta)}{m-\mu\cos^2\delta}\\ \dot\omega_2=\frac{-\cos\delta(-m\frac{g}{l_1}\sin\theta_1+\mu\lambda\omega_2^2\sin\delta)+m(-\frac{g}{l_1}\sin\theta_2-\omega_1^2\sin\delta)}{m\lambda-\mu\lambda\cos^2\delta} ω˙1=m−μcos2δ(−ml1gsinθ1+μλω22sinδ)−μcosδ(−l1gsinθ2−ω12sinδ)ω˙2=mλ−μλcos2δ−cosδ(−ml1gsinθ1+μλω22sinδ)+m(−l1gsinθ2−ω12sinδ)
从而转为代码得到
# 其中,lam,mu,G_L1,M为全局变量
def derivs(state, t):
dydx = np.zeros_like(state)
th1,om1,th2,om2 = state
dydx[0] = state[1]
delta = state[2] - state[0]
cDelta, sDelta = np.cos(delta), np.sin(delta)
sTh1,_,sTh2,_ = np.sin(state)
den1 = M - mu*cDelta**2
dydx[1] = (mu * om1**2 * sDelta * cDelta
+ mu * G_L1 * sTh2 * cDelta
+ mu * lam * om2**2 * sDelta
- M * G_L1 * sTh1)/ den1
dydx[2] = state[3]
den2 = lam * den1
dydx[3] = (- mu * lam * om2**2 * sDelta * cDelta
+ M * G_L1 * sTh1 * cDelta
- M * om1**2 * sDelta
- M * G_L1 * sTh2)/ den2
return dydx
接下来根据微分方程的解,便可进行绘图。
# 这段代码用于设置初值,并调用integrate求解微分方程组
import numpy as np
import scipy.integrate as integrate
G = 9.8
L1,L2 = 1.0, 1.0
G_L1 = G/L1
lam = L2/L1 #杆长度比L2/L1
mu = 1.0 #质量比M2/M1
M = 1+mu
# 生成时间
dt = 0.01
t = np.arange(0, 20, dt)
th1,th2 = 120.0, -10.0 #初始角度
om1,om2 = 0.0, 0.00 #初始角速度
state = np.radians([th1, om1, th2, om2])
# 微分方程组数值解
y = integrate.odeint(derivs, state, t)
# 真实坐标
x1 = L1*sin(y[:, 0])
y1 = -L1*cos(y[:, 0])
x2 = L2*sin(y[:, 2]) + x1
y2 = -L2*cos(y[:, 2]) + y1
至此,就得到了所有位置处的坐标,从而可以观察到双摆的轨迹如图所示
绘图代码为
import matplotlib.pyplot as plt
plt.scatter(x1,y1,marker='.')
plt.scatter(x2,y2,marker='.')
plt.show()
若将时间设置得长一点,然后在画图的时候更改一下颜色,就会看到双摆的运动区间,可见自然界还是挺有情怀的
其绘图代码为
plt.plot(x1,y1,marker='.',alpha=0.2, linewidth=0.2)
plt.plot(x2,y2,marker='.',alpha=0.2, linewidth=2, c='r')
plt.axis('off')
plt.show()
当然,也可以将其运动轨迹以一种三维的形式绘制出来
ax = plt.gca(projection='3d')
ax.plot3D(t,x1,y1,linewidth=1)
plt.show()
额……好吧,看来并没有什么情怀。
但是,如果把这两个小球分别当作两个星球,而我们又在一颗星球上,那么所观测到的另一颗星球的运动大致如下,不出意外是个圆,毕竟圆形二者之间的距离是恒定的。
绘图代码为
ax = plt.gca(projection='3d')
ax.plot3D(t,x2-x1,y2-y1,linewidth=0.5)
plt.show()
如果更改一下初值,则图形将有如下变化
初值设为
th1,th2 = 0, 0 #初始角度
om1,om2 = 120.0, 108.00 #初始角速度
运动过程
最后,还是传统技能,绘制一下双摆的运动过程如下
代码为
import matplotlib.animation as animation
# 下面为绘图过程
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect('equal')
ax.grid()
line, = ax.plot([], [], 'o-', lw=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
# 初始化图形
def init():
line.set_data([], [])
time_text.set_text('')
return line, time_text
def animate(i):
thisx = [0, x1[i], x2[i]]
thisy = [0, y1[i], y2[i]]
line.set_data(thisx, thisy)
time_text.set_text(time_template % (i*dt))
return line, time_text
ani = animation.FuncAnimation(fig, animate, range(1, len(y)),
interval=dt*1000, blit=True, init_func=init)
ani.save("dua_1.gif",writer='imagemagick')
plt.show()
附 公式推导过程
双摆的动能和势能分别为
T = ∑ i ∈ { 1 , 2 } m i 2 ( x ˙ i 2 + y ˙ i 2 ) V = − ∑ i ∈ { 1 , 2 } m i g y i \begin{aligned} T=&\sum_{i\in\{1,2\}}\frac{m_i}{2}(\dot x^2_i+\dot y^2_i)\\ V=&-\sum_{i\in\{1,2\}}m_igy_i \end{aligned} T=V=i∈{ 1,2}∑2mi(x˙i2+y˙i2)−i∈{ 1,2}∑migyi
系统的拉格朗日量为 L = T − V L=T-V L=T−V,带入角度坐标可得
L = m 1 + m 2 2 l 1 2 θ ˙ 1 2 + m 2 2 l 2 2 θ ˙ 2 2 + m 2 l 1 l 2 θ ˙ 1 θ ˙ 2 cos ( θ 1 − θ 2 ) + ( m 1 + m 2 ) l 1 g cos θ 1 + m 2 l 2 g cos θ 2 L=\frac{m_1+m_2}{2}l_1^2\dot\theta_1^2+\frac{m_2}{2}l^2_2\dot\theta_2^2+m_2l_1l_2\dot\theta_1\dot\theta_2\cos(\theta_1-\theta_2)\\+(m_1+m_2)l_1g\cos\theta_1+m_2l_2g\cos\theta_2 L=2m1+m2l12θ˙12+2m2l22θ˙22+m2l1l2θ˙1θ˙2cos(θ1−θ2)+(m1+m2)l1gcosθ1+m2l2gcosθ2
根据拉格朗日方程
d d t ∂ L ∂ θ ˙ i − ∂ L ∂ θ i = 0 \frac{\text d}{\text dt}\frac{\partial L}{\partial\dot\theta_i}-\frac{\partial L}{\partial\theta_i}=0 dtd∂θ˙i∂L−∂θi∂L=0
则有
d d t ( ( m 1 + m 2 ) l 1 2 θ ˙ 1 + m 2 l 1 l 2 θ ˙ 2 cos ( θ 1 − θ 2 ) ) = − m 2 l 1 l 2 θ ˙ 1 θ ˙ 2 sin ( θ 1 − θ 2 ) − ( m 1 + m 2 ) l 1 g sin θ 1 d d t ( m 2 l 2 2 θ ˙ 2 − m 2 l 1 l 2 θ ˙ 1 cos ( θ 1 − θ 2 ) ) = m 2 l 1 l 2 θ ˙ 1 θ ˙ 2 sin ( θ 1 − θ 2 ) − m 2 l 2 g sin θ 2 \frac{\text d}{\text dt}((m_1+m_2)l_1^2\dot\theta_1+m_2l_1l_2\dot\theta_2\cos(\theta_1-\theta_2))=-m_2l_1l_2\dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2)-(m_1+m_2)l_1g\sin\theta_1\\ \frac{\text d}{\text dt}(m_2l_2^2\dot\theta_2-m_2l_1l_2\dot\theta_1\cos(\theta_1-\theta_2))=m_2l_1l_2\dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2)-m_2l_2g\sin\theta_2\\ dtd((m1+m2)l12θ˙1+m2l1l2θ˙2cos(θ1−θ2))=−m2l1l2θ˙1θ˙2sin(θ1−θ2)−(m1+m2)l1gsinθ1dtd(m2l22θ˙2−m2l1l2θ˙1cos(θ1−θ2))=m2l1l2θ˙1θ˙2sin(θ1−θ2)−m2l2gsinθ2
记 μ = m 2 m 1 , λ = l 2 l 1 , δ = θ 1 − θ 2 \mu=\frac{m_2}{m_1},\lambda=\frac{l_2}{l_1}, \delta=\theta_1-\theta_2 μ=m1m2,λ=l1l2,δ=θ1−θ2,可简化为
( 1 + μ ) θ ¨ 1 + μ λ [ θ ¨ 2 cos δ − θ ˙ 2 2 sin δ ] = − ( 1 + μ ) g l 1 sin θ 1 λ θ ¨ 2 + θ ¨ 1 cos δ + θ ˙ 1 2 sin δ = − g l 1 sin θ 2 (1+\mu)\ddot\theta_1+\mu\lambda[\ddot\theta_2\cos\delta-\dot\theta_2^2\sin\delta] =-(1+\mu)\frac{g}{l_1}\sin\theta_1\\ \lambda\ddot\theta_2+\ddot\theta_1\cos\delta+\dot\theta_1^2\sin\delta=-\frac{g}{l_1}\sin\theta_2\\ (1+μ)θ¨1+μλ[θ¨2cosδ−θ˙22sinδ]=−(1+μ)l1gsinθ1λθ¨2+θ¨1cosδ+θ˙12sinδ=−l1gsinθ2
令 ω 1 = θ ˙ 1 , ω 2 = θ ˙ 2 \omega_1=\dot\theta_1,\omega_2=\dot\theta_2 ω1=θ˙1,ω2=θ˙2,记 m = 1 + μ m=1+\mu m=1+μ,上式可写为矩阵形式
[ m μ λ cos δ cos δ λ ] [ ω ˙ 1 ω ˙ 2 ] = [ − m g l 1 sin θ 1 + μ λ ω 2 2 sin δ − g l 1 sin θ 2 − ω 1 2 sin δ ] \begin{bmatrix}m&\mu\lambda\cos\delta\\ \cos\delta&\lambda \end{bmatrix}\begin{bmatrix}\dot\omega_1\\ \dot\omega_2 \end{bmatrix}=\begin{bmatrix} -m\frac{g}{l_1}\sin\theta_1+\mu\lambda\omega_2^2\sin\delta\\ -\frac{g}{l_1}\sin\theta_2-\omega_1^2\sin\delta \end{bmatrix} [mcosδμλcosδλ][ω˙1ω˙2]=[−ml1gsinθ1+μλω22sinδ−l1gsinθ2−ω12sinδ]
记 ω = [ ω 1 , ω 2 ] T \omega=[\omega_1,\omega_2]^T ω=[ω1,ω2]T, M = [ m μ λ cos δ cos δ λ ] M=\begin{bmatrix}m&\mu\lambda\cos\delta\\\cos\delta&\lambda\end{bmatrix} M=[mcosδμλcosδλ], F = [ − m g l 1 sin θ 1 + μ λ ω 2 2 sin δ − g l 1 sin θ 2 − ω 1 2 sin δ ] F=\begin{bmatrix}-m\frac{g}{l_1}\sin\theta_1+\mu\lambda\omega_2^2\sin\delta\\ -\frac{g}{l_1}\sin\theta_2-\omega_1^2\sin\delta\end{bmatrix} F=[−ml1gsinθ1+μλω22sinδ−l1gsinθ2−ω12sinδ],则
M ω ˙ = F → ω ˙ = M − 1 F M\dot\omega=F\to\dot\omega=M^{-1}F Mω˙=F→ω˙=M−1F
其中,
M − 1 = 1 m λ − μ λ cos 2 δ [ λ − μ λ cos δ − cos δ m ] M^{-1}=\frac{1}{m\lambda-\mu\lambda\cos^2\delta}\begin{bmatrix}\lambda&-\mu\lambda\cos\delta\\-\cos\delta&m\end{bmatrix} M−1=mλ−μλcos2δ1[λ−cosδ−μλcosδm]
展开可得则
ω ˙ 1 = ( − m g l 1 sin θ 1 + μ λ ω 2 2 sin δ ) − μ cos δ ( − g l 1 sin θ 2 − ω 1 2 sin δ ) m − μ cos 2 δ ω ˙ 2 = − cos δ ( − m g l 1 sin θ 1 + μ λ ω 2 2 sin δ ) + m ( − g l 1 sin θ 2 − ω 1 2 sin δ ) m λ − μ λ cos 2 δ \dot\omega_1=\frac{(-m\frac{g}{l_1}\sin\theta_1+\mu\lambda\omega_2^2\sin\delta)-\mu\cos\delta(-\frac{g}{l_1}\sin\theta_2-\omega_1^2\sin\delta)}{m-\mu\cos^2\delta}\\ \dot\omega_2=\frac{-\cos\delta(-m\frac{g}{l_1}\sin\theta_1+\mu\lambda\omega_2^2\sin\delta)+m(-\frac{g}{l_1}\sin\theta_2-\omega_1^2\sin\delta)}{m\lambda-\mu\lambda\cos^2\delta} ω˙1=m−μcos2δ(−ml1gsinθ1+μλω22sinδ)−μcosδ(−l1gsinθ2−ω12sinδ)ω˙2=mλ−μλcos2δ−cosδ(−ml1gsinθ1+μλω22sinδ)+m(−l1gsinθ2−ω12sinδ)