雙重簡諧運動

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

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

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

模擬程式畫面截圖

程式 14-1.水平彈簧與小球組成的雙重簡諧運動（取得程式碼）　（GlowScript 網站動畫連結）








































































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

from vpython import *

"""
 1. 參數設定, 設定變數及初始值, 變數組合
"""
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 到達地板右側邊緣時停止
"""
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()
```
若彈簧伸長時形變量為正值，此時右側小球被彈簧向左拉，上式中得到的力量與彈簧的軸方向相反，與實際情形相符。
計算小球的加速度，更新速度、位置，計算木塊的動能、系統彈性位能、總能量，畫速度 - 時間關係圖及能量 - 時間關係圖，更新時間。