一個地球人在臺北: 雙重簡諧運動

作者：王一哲
日期：2018/4/26

假設在光滑水平桌面上有兩個小球，質量分別為 m1 及 m2，用一條彈力常數為 k 的理想彈簧連結。若施力敲擊其中一個小球，使小球獲得動量，則整個系統會以類似毛毛蟲爬行的方式前進。我以前看到的動畫是用 Mathematica 製作的，這次我們改用 VPython 達成同樣的效果。

模擬程式畫面截圖

程式 14-1.水平彈簧與小球組成的雙重簡諧運動

取得程式碼
GlowScript 網站動畫連結

"""
 VPython教學: 14-1.水平彈簧與小球組成的雙重簡諧運動
 Ver. 6: 2019/9/13
 作者: 王一哲
"""

from vpython import *

"""
 1. 參數設定, 設定變數及初始值, 變數組合
    (1) m1=0.1, m2=0.2, v1=4, v2=0, k=5
    (2) m1=0.1, m2=0.2, v1=0, v2=4, k=5
    (3) m1=0.2, m2=0.1, v1=4, v2=0, k=5
    (4) m1=0.2, m2=0.1, v1=0, v2=4, k=5
"""
r1, m1, v1, c1 = 0.1, 0.1, 4, color.red      # 球1的半徑、質量、初速、顏色
r2, m2, v2, c2 = 0.1, 0.2, 0, color.green    # 球2的半徑、質量、初速、顏色
L0, k = 1.0, 5.0                             # 彈簧的原長 L0=1 m, 彈性常數 k=5.0 N/m
xmax, xmin = 2.0, 0.0                        # x 軸範圍
t, dt = 0, 0.00005                           # 時間, 時間間隔，原為0.001但能量不夠準確, 故改為0.00005

"""
 2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="Double Simple Harmonic Motion", width=800, height=300, center=vec(0, 0.4, 0),
               background=vec(0, 0.6, 0.6))
# 產生地板
floor = box(pos=vec(0, -1.2*r1, 0), size=vec(2.0*(xmax - xmin), 0.05, 0.8), color=color.blue)
# 產生左側小球 b1, 右側小球 b2 並設定初速度
b1 = sphere(pos=vec(-L0, 0, 0), radius=r1, color=c1, v=vec(v1, 0, 0))
b2 = sphere(pos=vec(0, 0, 0), radius=r2, color=c2, v=vec(v2, 0, 0))
# 產生彈簧, 起點為 b1.pos, 方向為 b2.pos - b1.pos
spring = helix(pos=b1.pos, axis=b2.pos - b1.pos, radius=0.05, thickness=0.03)
# 繪圖部分
gd = graph(title="<i>E</i> - <i>t</i> plot", x=0, y=300, width=600, height=450, xtitle="<i>t</i> (s)",
           ytitle="red: <i>K</i><sub>1</sub>, green: <i>K</i><sub>2</sub>, orange: <i>U</i>, cyan: <i>E</i> (J)")
kt1 = gcurve(graph=gd, color=c1)
kt2 = gcurve(graph=gd, color=c2)
ut = gcurve(graph=gd, color=color.orange)
et = gcurve(graph=gd, color=color.cyan)
gd2 = graph(title="<i>v</i> - <i>t</i> plot", x=0, y=750, width=600, height=450, xtitle="<i>t</i> (s)",
            ytitle="red: <i>v</i><sub>1</sub>, green: <i>v</i><sub>2</sub> (m/s)")
vt1 = gcurve(graph=gd2, color=c1)
vt2 = gcurve(graph=gd2, color=c2)

"""
 3. 物體運動部分, b2 到達地板右側邊緣時停止
    由虎克定律求彈簧的回復力 force=-k*delta_x
    為了使 force 變為向量, 需要乘以 spring.axis 的單位向量
    mag(a)=a.mag => 計算向量 a 的量值
    a / mag(a)=a.norm()=norm(a) => 回傳向量 a 的單位向量
"""
while(b2.pos.x <= xmax):
    rate(1000)
# 更新彈簧的起點位置、長度、方向
    spring.pos = b1.pos
    spring.axis = b2.pos - b1.pos
# 計算彈簧回復力，更新小球的加速度、速度、位置
    force = -k*(spring.axis.mag - L0) * spring.axis.norm()
    b1.a = -force/m1     
    b2.a = force/m2
    b1.v += b1.a * dt
    b2.v += b2.a * dt
    b1.pos += b1.v * dt
    b2.pos += b2.v * dt
# 計算小球的動能、系統彈性位能、力學能並畫圖
    k1 = 0.5 * m1 * b1.v.mag2
    k2 = 0.5 * m2 * b2.v.mag2
    u = 0.5 * k * (spring.axis.mag - L0)**2
    e = k1 + k2 + u
    kt1.plot(pos=(t, k1))
    kt2.plot(pos=(t, k2))
    ut.plot(pos=(t, u))
    et.plot(pos=(t, e))
# 畫 v-t 圖
    vt1.plot(pos=(t, b1.v.x))
    vt2.plot(pos=(t, b2.v.x))
# 更新時間
    t += dt

參數設定

在此設定變數為小球的半徑、質量、初速、顏色，彈簧的原長、彈性常數，x軸的範圍、時間、時間間隔，其中時間間隔 dt 設定為 0.00005，這是因為設定為 0.001 時計算系統能量的誤差較大，故選擇較小的數值。

畫面設定

產生動畫視窗、地板、小球、繪圖視窗的程式碼在之前動畫當中已經出現很多次，這裡就不再贅述。之前比較少用到的物件是 helix，可以調整的選項主要有

位置(pos)，以螺旋線的端點為基準。
半徑(radius)
軸(axis)，以位置為起點，格式為向量。
厚度(thickness)，線條本身的半徑，預設值為 radius/20。
顏色(color)

物體運動

為了使動畫重覆執行，直到 b2 到達地板右側邊緣時停止，while 迴圈中的條件設定為
```
b2.pos.x <= xmax
```
更新彈簧的起點位置，由於彈簧左端連接在 b1 球心處，因此
```
spring.pos = b1.pos
```
更新彈簧的長度及方向，也就是更新 axis，由於彈簧連接在兩個小球之間，因此
```
spring.axis = b2.pos - b1.pos
```
由虎克定律計算彈簧的回復力。彈簧的形變量為 spring.axis.mag - L0，若彈簧壓縮形變量為負值，為了使右側小球被彈簧向右推出，計算時要加上負號並乘以彈簧軸方向的單位向量，程式碼為
```
force = -k*(spring.axis.mag - L0) * spring.axis.norm()
```
若彈簧伸長時形變量為正值，此時右側小球被彈簧向左拉，上式中得到的力量與彈簧的軸方向相反，與實際情形相符。
計算小球的加速度，更新速度、位置，計算木塊的動能、系統彈性位能、總能量，畫速度 - 時間關係圖及能量 - 時間關係圖，更新時間。