Python艺术|用Python画双摆动力系统|双摆|matplotlib|scipy|

文章目录

- 双摆问题
- 运动过程
- 附公式推导过程

双摆问题所谓双摆，就是两个连在一起的摆。

文章图片

设小球 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?=l1?sinθ1?y1?=?l1?cosθ1?x2?=x1?+l2?sinθ2?y2?=y1??l2?cosθ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?=l1?sinθ1?y1?=?l1?cosθ1??x2?=x1?+l2?sinθ2?y2?=y1??l2?cosθ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δ(?ml1?g?sinθ1?+μλω22?sinδ)?μcosδ(?l1?g?sinθ2??ω12?sinδ)?ω˙2?=mλ?μλcos2δ?cosδ(?ml1?g?sinθ1?+μλω22?sinδ)+m(?l1?g?sinθ2??ω12?sinδ)?
从而转为代码得到

# 其中，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 integrateG = 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_textdef 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_textani = 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}∑?mi?gyi??
系统的拉格朗日量为 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?+m2??l12?θ˙12?+2m2??l22?θ˙22?+m2?l1?l2?θ˙1?θ˙2?cos(θ1??θ2?)+(m1?+m2?)l1?gcosθ1?+m2?l2?gcosθ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?+m2?l1?l2?θ˙2?cos(θ1??θ2?))=?m2?l1?l2?θ˙1?θ˙2?sin(θ1??θ2?)?(m1?+m2?)l1?gsinθ1?dtd?(m2?l22?θ˙2??m2?l1?l2?θ˙1?cos(θ1??θ2?))=m2?l1?l2?θ˙1?θ˙2?sin(θ1??θ2?)?m2?l2?gsinθ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 μ=m1?m2??，λ=l1?l2??,δ=θ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?+μλ[θ¨2?cosδ?θ˙22?sinδ]=?(1+μ)l1?g?sinθ1?λθ¨2?+θ¨1?cosδ+θ˙12?sinδ=?l1?g?sinθ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??]=[?ml1?g?sinθ1?+μλω22?sinδ?l1?g?sinθ2??ω12?sinδ?]
记 ω = [ ω 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=[?ml1?g?sinθ1?+μλω22?sinδ?l1?g?sinθ2??ω12?sinδ?]，则
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?]
展开可得则
【Python艺术|用Python画双摆】 ω ˙ 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δ(?ml1?g?sinθ1?+μλω22?sinδ)?μcosδ(?l1?g?sinθ2??ω12?sinδ)?ω˙2?=mλ?μλcos2δ?cosδ(?ml1?g?sinθ1?+μλω22?sinδ)+m(?l1?g?sinθ2??ω12?sinδ)?