利用Python繪畫雙擺操作分享

更新時間：2022年04月01日 08:36:11 作者：微小冷

這篇文章主要介紹了利用Python畫雙擺，繪畫雙擺的過程主要包括以下步驟，雙擺問題、運(yùn)動過程及公式推導(dǎo)過程，下文詳細(xì)介紹，需要的小伙伴可以參考一下

1.雙擺問題

所謂雙擺，就是兩個連在一起的擺。

接下來本來是要推公式的，考慮考慮到大家可能會有公式恐懼癥，同時又喜歡看圖，所以把公式挪到后面。

所以，只需知道角速度的微分方程，就可寫出對應(yīng)的代碼，其方程如下:

從而轉(zhuǎn)為代碼得到：

# 其中，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

接下來根據(jù)微分方程的解，便可進(jìn)行繪圖。

# 這段代碼用于設(shè)置初值，并調(diào)用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 ? ? ?#質(zhì)量比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])

# 微分方程組數(shù)值解
y = integrate.odeint(derivs, state, t)

# 真實(shí)坐標(biāo)
x1 = L1*sin(y[:, 0])
y1 = -L1*cos(y[:, 0])
x2 = L2*sin(y[:, 2]) + x1
y2 = -L2*cos(y[:, 2]) + y1

至此，就得到了所有位置處的坐標(biāo)，從而可以觀察到雙擺的軌跡如圖所示

繪圖代碼為：

import matplotlib.pyplot as plt
plt.scatter(x1,y1,marker='.')
plt.scatter(x2,y2,marker='.')
plt.show()

若將時間設(shè)置得長一點(diǎn)，然后在畫圖的時候更改一下顏色，就會看到雙擺的運(yùn)動區(qū)間，可見自然界還是挺有情懷的

其繪圖代碼為：

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()

當(dāng)然，也可以將其運(yùn)動軌跡以一種三維的形式繪制出來

ax = plt.gca(projection='3d')
ax.plot3D(t,x1,y1,linewidth=1)
plt.show()

額……好吧，看來并沒有什么情懷。

但是，如果把這兩個小球分別當(dāng)作兩個星球，而我們又在一顆星球上，那么所觀測到的另一顆星球的運(yùn)動大致如下，不出意外是個圓，畢竟圓形二者之間的距離是恒定的。

繪圖代碼為：

ax = plt.gca(projection='3d')
ax.plot3D(t,x2-x1,y2-y1,linewidth=0.5)
plt.show()

如果更改一下初值，則圖形將有如下變化

初值設(shè)為：

th1,th2 = 0, 0 ?#初始角度
om1,om2 = 120.0, 108.00 ? ? ? #初始角速度

2.運(yùn)動過程

最后，還是傳統(tǒng)技能，繪制一下雙擺的運(yùn)動過程如下：

代碼為：

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()