# VPython 進階教材：小球鉛直簡諧運動與橫波
> 作者：王一哲
> 第1版：2020/10/19
> 第2版：2020/12/5 補充不使用 NumPy 的寫法
<br />
## 前言
當我需要繪製橫波的動畫時，通常是使用 GeoGebra 搭配以下的數學式
$$
f(x) = R \sin(kx - \omega t) = R \sin \left(\frac{2 \pi}{\lambda} \cdot x - \frac{2 \pi}{T} \cdot t \right)
$$
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://lh4.googleusercontent.com/f-3MZb7NX3Fhdz2rJYLwSbItnIh022PsZQRTHnQiw5ZqJrxzuP-RNEzWXsSJSR9Sowi3fmIpXbXv4Q2M_xvD67TGg5lmeXBZxt5kxQtlb01mlEihJySO55ReTpvXegXkQITKm6S5">
<div style="text-align:center">使用 GeoGebra 繪製的橫波動畫</div>
<br />
橫波上的每個質點在鉛直方向做簡諧運動，如果想要用 VPython 做出一樣的效果，理論上只要讓一排的小球每隔固定的時間差落到正下方的彈簧上，小球與彈簧結合後開始做鉛直簡諧運動即可。為了盡量避免使用 for 迴圈，我試著使用 numpy 陣列 (ndarray) 取代大部分的串列 (list)，以下是我測試成功的動畫及程式碼。
<div style="text-align:center"><iframe src="https://www.glowscript.org/#/user/yizhe/folder/Public/program/ParticleWave" width="600" height="600"></iframe></div>
<div style="text-align:center">使用 VPython 畫出小球做鉛直簡諧運動模擬橫波的動畫</div>
<br />
## VPython 程式碼
```python=
"""
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
```
<br />
1. 第 28 ~ 33 行：這裡利用了幾個 numpy 陣列儲存資料，xs 是小球與彈簧的 x 坐標，ys 是小球的起始高度，vs 是小球的初速度，ay 是小球的 y 方向加速度，delays 是小球開始落下的時間，states 是小球與彈簧是否接觸的狀態。使用到的 numpy 陣列函數有
```python
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。np.nonzero(states == True) 也可以簡化成 np.nonzero(states)。
5. 第 48 行：當時間 t 大於等於該小球的時間延遲值時，更新小球的速度。
6. 第 50 ~ 53 行：使用 for 迴圈更新小球的高度, 若小球與彈簧已接觸則更新彈簧的軸向量為 (0, ys[i]+L0, 0)。
<br />
## 不使用 NumPy 的寫法
```python=
"""
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]*num
vs = [0]*num
ay = [-g]*num
delays = [i*3*T/(num-1) for i in range(num)]
states = [False]*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]:
springs[i].axis = vec(0, ys[i]+L0, 0)
t += dt
```
<br />
主要修改以下兩處，修改後的程式碼可以在 GlowScript 網站上執行，這是 [GlowScript 網站動畫連結](https://www.glowscript.org/#/user/yizhe/folder/Public/program/ParticleWave)。
1. 第 27、31 行：將原來的 numpy 陣列全部改成 list，並且使用單行的 for 迴圈縮短程式碼。計算 xs 及 delays 時需要除以 num-1 而不是 num，因為 41 個小球之間只有 40 個間隔。
2. 第 45 ~ 55 行：使用 for 迴圈逐一更新小球及彈簧的狀態。
<br />
## 結語
由於 numpy 陣列有許多很神奇的函數及使用方法，如果能夠善加利用可以少寫很多個 for 迴圈，有效提升程式的運算速度。但是 VPython 的物件只能儲存到串列中，因此在更新整串的 VPython 物件時仍然需要使用 for 迴圈，這是無法避免的。
<br />
---
###### tags:`VPython`
Sign in with Wallet
Connect another wallet