VPython進階教學：蛇擺 Pendulum Wave

作者：王一哲
日期：2018/7/14




最近看到建國中學曾靖夫老師的 VPython 講義〈Lecture 7 多物件控制（list） - 波動現象模擬〉，其中一項作業就是要學生做出蛇擺的動畫。蛇擺的英文為 pendulum wave，如果照字面翻譯應該譯為擺波，由一組數個的單擺構成的波，假設整組蛇擺的週期為 $T$，其中週期最長的單擺在這段時間內共擺動 $N$ 次，則週期次長的單擺擺動 $N-1$ 次，依此類推。我們之前已經做過〈單擺〉的動畫，可以以此為基礎改造一下程式碼，應該很容易就能做出蛇擺的動畫。關於蛇擺的物理原理介紹請參考中央大學科學教育中心的網頁：〈蛇擺〉、〈大型蛇擺〉。

這是我用 VPython 在 Glowscript 網站製作的動畫，連結為 https://www.glowscript.org/#/user/yizhe/folder/Public/program/PendulumWave。


這個程式是進階教材，我們要善用類別 (class)的特性，使程式在使用上更加方便。程式的目標：

1. 輸入蛇擺的週期、週期最長的單擺在蛇擺的一個週期內擺動的次數、單擺個數，自動計算所有單擺的週期及擺長。
2. 建立一個類別 (class)，輸入單擺的週期、懸掛位置，自動產生擺錘、繩子，掛於指定的位置，將擺錘放置於最大擺角。計算單擺對應的轉動慣量，設定起始的角速度、角加速度為0。
3. 在類別中建立一個方法 (mehtod)，輸入經過的時間，計算擺錘所受力矩、角加速度、角速度，更新擺錘的位置及繩子的方向。

VPython教學進階內容: 蛇擺

"""
 VPython教學進階內容: 蛇擺 pendulum wave, 使用 class
 Ver. 1: 2018/7/12
 Ver. 2: 2018/7/14 加上計時器
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
g = 9.8             # 重力加速度
Tpw = 30            # 蛇擺的週期
N = 15              # 週期最長的單擺在蛇擺的1個週期內擺動的次數
Tmax = Tpw / N      # 週期最長的單擺擺動週期
Lmax = Tmax**2 * g / (4 * pi**2)   # 週期最長的單擺擺長 
num = 20            # 單擺的個數
width = Lmax        # 將畫面邊長設定為最長的單擺擺長
m = 1               # 小球質量
size = width / (2 * num)           # 小球半徑
theta0 = radians(30)# 起始擺角, 用 radians 將單位換成 rad
i = 0               # 小球經過週期次數
t = 0               # 時間
dt = 0.001          # 時間間隔

"""
 2. 畫面設定
"""
# 產生動畫視窗、天花板
scene = canvas(title="Pendulum Wave", width=600, height=600, x=0, y=0, background=color.black)
scene.camera.pos = vec(-1.5*width, 0.5*width, width)
scene.camera.axis = vec(1.5*width, -0.5*width, -width)
roof = box(pos=vec(0, (width + size)/2, 0), size=vec(width, size, 0.5*width), color=color.cyan)
timer = label(pos=vec(0.9*width, 0.9*width, 0), text="t =  s", space=50, height=24,
              border=4, font="monospace")

# 新增類別 Pendulum, 輸入週期T, 懸掛位置loc, 編號idx, 自動產生小球及對應的繩子
# 設定方法 update, 輸入經過的時間dt, 更新單擺狀態
class Pendulum:
    def __init__(self, T, loc, idx):
        self.T = T
        self.loc = loc
        self.idx = idx
        self.L = self.T**2 * g / (4 * pi**2)
        self.I = m * self.L**2
        self.alpha = 0
        self.omega = 0
        self.theta = theta0
        self.ball = sphere(pos=vec(self.loc, width/2 - self.L*cos(theta0), self.L*sin(theta0)), 
                           radius=size, color=vec(1 - self.idx/num, 0, self.idx/num))
        self.rope = cylinder(pos=vec(self.loc, width/2, 0), axis=self.ball.pos - vec(self.loc, width/2, 0),
                             radius=0.1*size, color=color.yellow)
    def update(self, dt):
        self.dt = dt
        self.alpha = -m*g*self.ball.pos.z/self.I
        self.omega += self.alpha*self.dt
        self.theta += self.omega*self.dt
        self.ball.pos = vec(self.loc, 0.5*width - self.L*cos(self.theta), self.L*sin(self.theta))
        self.rope.axis = self.ball.pos - vec(self.loc, 0.5*width, 0)

# 利用自訂類別 Pendulum 產生 num 個單擺
pendulums = []
for i in range(num):
    T = Tpw / (N + i)
    loc = width*(-0.5+(i/(num-1)))
    pendulum = Pendulum(T, loc, i)
    pendulums.append(pendulum)

"""
 3. 物體運動部分
"""
while True:
    rate(500)
    for pendulum in pendulums:
        pendulum.update(dt)
    timer.text = "t = {:.1f} s".format(t)
    t += dt

程式設計部分

1. 首先設定一些參數，主要會調整的是蛇擺的週期Tpw，週期最長的單擺在蛇擺的1個週期裡擺動的次數N，單擺的數量num。
2. 產生動畫視窗及天花板，調整動畫視窗的視角。
3. 用 label 物件於畫面右上角加上計時器，方便使用者看出蛇擺的週期。
4. 自訂類別 Pendulum。
   
   1. 定義初始化方法 __init__，輸入的參數為 self, T, loc, idx，其中 self 是建立類別時一定要有的，但在使用類別建立物件時不用輸入。T是單擺週期，loc是懸掛位置，idx是單擺編號。利用類別建立物件時會執行 __init__ ，程式會計算對應的擺長，自動產生擺錘、繩子，掛於指定的位置，將擺錘放置於最大擺角。計算單擺對應的轉動慣量，設定起始的角速度、角加速度為0。
   2. 定義方法 update，輸入的參數為經過的時間 dt，使用時程式會計算擺錘所受力矩、角加速度、角速度，更新擺錘的位置及繩子的方向。
5. 使用自訂類別 Pendulum 以及 for 迴圈產生 num 個單擺，存入串列 pendulums。


6. 物體運動部分，由於我們希望動畫持續執行，因此使用 while True，讓 while 迴圈不會結束。再用 for 迴圈將串列 pendulums 當中的單擺一次取出一個，利用自訂類別中的方法 update 更新狀態，同時更新計時器中顯示的時間。由於我們把大部分的功能都寫在自訂類別當中，所以物體運動部分的程式碼非常地少。


7. 由於時間 t 是浮點數，如果顯示所有的位數時會使得計時器的位置一直左右跳動，為了避免這個現象我們使用以下的語法限制輸出的位數到小數點以下第1位為止，兩種語法皆可，但建議使用第1種語法。
   

   timer.text = "t = {:.1f} s".format(t)
   timer.text = "t = " + str('%.1f' % t) + " s"
   

   舉例來說，假設變數 a = 0.123456789，但如果寫成
   

   print("{:.6f}".format(a))
   print('%.6f' % a)
   

   python 會顯示 0.123457。

模擬結果

我將程式碼上傳到 GlowScript 網站，請按以下連結觀看動畫http://www.glowscript.org/#/user/yizhe/folder/Public/program/PendulumWave GlowScript 網站在指定浮點數顯示位數的部分，只支援第1種語法。

蛇擺模擬動畫（第1版，沒有計時器）

蛇擺模擬動畫（第2版，有計時器）

蛇擺模擬動畫（GlowScript版）

結語

我們在這個動畫裡利用自訂的 Class 產生單擺，在使用上比較方便但在程式設計上比較麻煩，如果只是想做出同樣的動畫，可以使用 list 儲存單擺的週期、擺長等資料，再用 for 迴圈讀取、執行。Python 裡有許多好用但有點難懂的工具，例如這次用到的Class，還有科學計算套件Numpy、繪圖套件Matplotlib、讀取資料套件Pandas……等等，如果之後的動畫中需要用到會再寫講義，敬請期待。

---

HackMD 版本連結：https://hackmd.io/@yizhewang/SyEWFOEmm