小球鉛直簡諧運動與橫波

作者：王一哲
日期：2020/10/19

前言

當我需要繪製橫波的動畫時，通常是使用 GeoGebra 搭配以下的數學式 $$ f(x) = R \sin(kx - \omega t) = R \sin \left(\frac{2 \pi}{\lambda} \cdot x - \frac{2 \pi}{T} \cdot t \right) $$

使用 GeoGebra 繪製的橫波動畫

橫波上的每個質點在鉛直方向做簡諧運動，如果想要用 VPython 做出一樣的效果，理論上只要讓一排的小球每隔固定的時間差落到正下方的彈簧上，小球與彈簧結合後開始做鉛直簡諧運動即可。為了盡量避免使用 for 迴圈，我試著使用 numpy 陣列 (ndarray) 取代大部分的串列 (list)，以下是我測試成功的動畫及程式碼。

使用 VPython 畫出小球做鉛直簡諧運動模擬橫波的動畫

VPython 程式碼

"""
 VPython教學: 使用多個質點的鉛直簡諧運動組成橫波
 日期: 2020/10/19
 作者: 王一哲
"""
from vpython import *
import numpy as np

"""
 1. 參數設定, 設定變數及初始值
"""
m, r, h = 0.1, 0.3, 4   # 小球質量, 半徑, 起始高度
g = 9.8             # 重力加速度
k, L0 = 1.0, 8      # 彈簧彈性常數, 原長
T = 2*pi*sqrt(m/k)  # 簡諧運動週期
t, dt = 0, 0.0005   # 時間, 時間間隔
length, num = 8, 41 # 波長, 小球個數

"""
 2. 畫面設定
"""
# 動畫視窗
scene = canvas(title="Particles and Wave", width=600, height=600, x=0, y=0,
               background=color.black, range=1.7*length, center=vec(1.5*length, 0, 0))
# 地板
floor = box(pos=vec(1.5*length, -L0, 0), size=vec(3*length, r, length), color=color.blue)

xs = np.linspace(0, 3*length, num)   # x 坐標
ys = np.ones(num)*h                  # y 坐標
vs = np.zeros(num)                   # y 方向速度
ay = np.ones(num)*(-g)               # y 方向加速度
delays = np.linspace(0, 3*T, num)    # 時間延遲
states = np.full(num, False)         # 小球與彈簧是否接觸

# 小球
balls = [sphere(pos=vec(xs[i], ys[i], 0), radius=r, color=color.red) for i in range(num)]

# 彈簧
springs = [helix(pos=vec(xs[i], -L0, 0), radius=0.6*r, thickness=0.4*r, axis=vec(0, L0, 0), color=color.yellow) for i in range(num)]

"""
 3. 物體運動部分
"""
while True:
    rate(1000)
    states[ys <= 0] = True
    ay[np.nonzero(states == True)] = -g - k*ys[np.nonzero(states == True)]/m
    vs[delays <= t] += ay[delays <= t]*dt
    ys += vs*dt
    for i in range(num):
        balls[i].pos = vec(balls[i].pos.x, ys[i], 0)
        if states[i]:
            springs[i].axis = vec(0, ys[i]+L0, 0)
    t += dt

1. 第 28 ~ 33 行：這裡利用了幾個 numpy 陣列儲存資料，xs 是小球與彈簧的 x 坐標，ys 是小球的起始高度，vs 是小球的初速度，ay 是小球的 y 方向加速度，delays 是小球開始落下的時間，states 是小球與彈簧是否接觸的狀態。使用到的 numpy 陣列函數有

   np.linspace(起始值, 結束值, 等分數量) # 將起始值到結束值之間按照數量等分並回傳為陣列
   np.zeros(長度) # 依照指定的長度產生所有元素皆為0的陣列
   np.ones(長度) # 依照指定的長度產生所有元素皆為1的陣列
   np.full(長度, 值) # 依照指定的長度產生所有元素皆為輸入值的陣列
   

2. 第 36 ~ 39 行：利用單行的 for 迴圈將小球 balls 及彈簧 springs 物件儲存到串列中。
3. 第 46 行：為了簡化算式，彈簧頂端起始高度為 y = 0，當 ys 由大於 0 變為小於等於 0 時，代表小球與彈簧接觸，state 指定為 True。
4. 第 47 行：當小球與彈簧接觸後, 彈簧形變量 = ys。當 ys < 0 時，彈簧回復力方向向上，因此 k*ys 需要加上負號。再使用 np.nonzero(states == True) 列出小球與彈簧已接觸的元素索引值，將加速度由 -g 更新為 -g - k*dy/m。
5. 第 48 行：當時間 t 大於等於該小球的時間延遲值時，更新小球的速度。
6. 第 50 ~ 53 行：使用 for 迴圈更新小球的高度, 若小球與彈簧已接觸則更新彈簧的軸向量為 (0, ys[i]+L0, 0)。

不使用 NumPy 的寫法

"""
 VPython教學: 使用多個質點的鉛直簡諧運動組成橫波, 不使用 NumPy
 日期: 2020/12/5
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m, r, h = 0.1, 0.3, 4   # 小球質量, 半徑, 起始高度
g = 9.8             # 重力加速度
k, L0 = 1.0, 8      # 彈簧彈性常數, 原長
T = 2*pi*sqrt(m/k)  # 簡諧運動週期
t, dt = 0, 0.0005   # 時間, 時間間隔
length, num = 8, 41 # 波長, 小球個數

"""
 2. 畫面設定
"""
# 動畫視窗
scene = canvas(title="Particles and Wave", width=600, height=600, x=0, y=0,
               background=color.black, range=1.7*length, center=vec(1.5*length, 0, 0))
# 地板
floor = box(pos=vec(1.5*length, -L0, 0), size=vec(3*length, r, length), color=color.blue)

xs = [i*3*length/(num-1) for i in range(num)]
ys = [h for i in range(num)]
vs = [0 for i in range(num)]
ay = [-g for i in range(num)]
delays = [i*3*T/(num-1) for i in range(num)]
states = [False for i in range(num)]

# 小球
balls = [sphere(pos=vec(xs[i], ys[i], 0), radius=r, color=color.red) for i in range(num)]

# 彈簧
springs = [helix(pos=vec(xs[i], -L0, 0), radius=0.6*r, thickness=0.4*r, axis=vec(0, L0, 0), color=color.yellow) for i in range(num)]

"""
 3. 物體運動部分
"""
while True:
    rate(1000)
    for i in range(num):
        if ys[i] <= 0:
            states[i] = True
        if states[i]:
            ay[i] = -g - k*ys[i]/m
        if delays[i] <= t:
            vs[i] += ay[i]*dt
        ys[i] += vs[i]*dt
        balls[i].pos = vec(balls[i].pos.x, ys[i], 0)
        if states[i] == True:
            springs[i].axis = vec(0, ys[i]+L0, 0)
    t += dt

主要修改以下兩處，修改後的程式碼可以在 GlowScript 網站上執行，這是 GlowScript 網站動畫連結。

1. 第 27 ~ 32 行：將原來的 numpy 陣列全部改成 list，並且使用單行的 for 迴圈縮短程式碼。計算 xs 及 delays 時需要除以 num-1 而不是 num，因為 41 個小球之間只有 40 個間隔。
2. 第 45 ~ 55 行：使用 for 迴圈逐一更新小球及彈簧的狀態。

結語

由於 numpy 陣列有許多很神奇的函數及使用方法，如果能夠善加利用可以少寫很多個 for 迴圈，有效提升程式的運算速度。但是 VPython 的物件只能儲存到串列中，因此在更新整串的 VPython 物件時仍然需要使用 for 迴圈，這是無法避免的。

---

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