# 雙重簡諧運動
> 作者：王一哲
> 日期：2018/4/26
<br />
假設在光滑水平桌面上有兩個小球，質量分別為 m1 及 m2，用一條彈力常數為 k 的理想彈簧連結。若施力敲擊其中一個小球，使小球獲得動量，則整個系統會以類似毛毛蟲爬行的方式前進。我以前看到的動畫是用 Mathematica 製作的，這次我們改用 VPython 達成同樣的效果。
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/EPqYpPq.png">
<div style="text-align:center">模擬程式畫面截圖</div>
<br />
## 程式 14-1.水平彈簧與小球組成的雙重簡諧運動 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/14.%E9%9B%99%E9%87%8D%E7%B0%A1%E8%AB%A7%E9%81%8B%E5%8B%95/14-1_double_SHM.py)） （[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/14-1doubleSHM)）
```python=
"""
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
```
<br />
### 參數設定
在此設定變數為小球的半徑、質量、初速、顏色，彈簧的原長、彈性常數，x軸的範圍、時間、時間間隔，其中時間間隔 dt 設定為 0.00005，這是因為設定為 0.001 時計算系統能量的誤差較大，故選擇較小的數值。
<br />
### 畫面設定
產生動畫視窗、地板、小球、繪圖視窗的程式碼在之前動畫當中已經出現很多次，這裡就不再贅述。之前比較少用到的物件是 <span style="color:red; font-weight:bold">helix</span>，可以調整的選項主要有
1. 位置(pos)，以螺旋線的端點為基準。
2. 半徑(radius)
3. 軸(axis)，以位置為起點，格式為向量。
4. 厚度(thickness)，線條本身的半徑，預設值為 radius/20。
5. 顏色(color)
<br />
### 物體運動
1. 為了使動畫重覆執行，直到 b2 到達地板右側邊緣時停止，while 迴圈中的條件設定為
```python
b2.pos.x <= xmax
```
2. 更新彈簧的起點位置，由於彈簧左端連接在 b1 球心處，因此
```python
spring.pos = b1.pos
```
更新彈簧的長度及方向，也就是更新 axis，由於彈簧連接在兩個小球之間，因此
```python
spring.axis = b2.pos - b1.pos
```
3. 由虎克定律計算彈簧的回復力。彈簧的形變量為 **spring.axis.mag - L0**，若彈簧壓縮形變量為負值，為了使右側小球被彈簧向右推出，計算時要加上負號並乘以彈簧軸方向的單位向量，程式碼為
```python
force = -k*(spring.axis.mag - L0) * spring.axis.norm()
```
若彈簧伸長時形變量為正值，此時右側小球被彈簧向左拉，上式中得到的力量與彈簧的軸方向相反，與實際情形相符。
4. 計算小球的加速度，更新速度、位置，計算木塊的動能、系統彈性位能、總能量，畫速度 - 時間關係圖及能量 - 時間關係圖，更新時間。
<br />
### 模擬結果
我測試了 4 種不同的條件：
1. m1 = 0.1, v1 = 4.0, m2 = 0.2, v2 = 0.0
2. m1 = 0.1, v1 = 0.0, m2 = 0.2, v2 = 4.0
3. m1 = 0.2, v1 = 4.0, m2 = 0.1, v2 = 0.0
4. m1 = 0.2, v1 = 0.0, m2 = 0.1, v2 = 4.0
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/BPqzGGp.png">
<div style="text-align:center">條件1的模擬結果：速度 - 時間關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/9xhHnQG.png">
<div style="text-align:center">條件1的模擬結果：能量 - 時間關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/LqXJKX3.png">
<div style="text-align:center">條件2的模擬結果：速度 - 時間關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/6aSMuvv.png">
<div style="text-align:center">條件2的模擬結果：能量 - 時間關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/JofCc6a.png">
<div style="text-align:center">條件3的模擬結果：速度 - 時間關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/NsSPlB9.png">
<div style="text-align:center">條件3的模擬結果：能量 - 時間關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/xtNSgeH.png">
<div style="text-align:center">條件4的模擬結果：速度 - 時間關係圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/ZBxfRH8.png">
<div style="text-align:center">條件4的模擬結果：能量 - 時間關係圖</div>
<br />
## VPython官方說明書
1. **canvas**: http://www.glowscript.org/docs/VPythonDocs/canvas.html
2. **box**: http://www.glowscript.org/docs/VPythonDocs/box.html
3. **sphere**: http://www.glowscript.org/docs/VPythonDocs/sphere.html
4. **helix**: http://www.glowscript.org/docs/VPythonDocs/helix.html
5. **graph**: http://www.glowscript.org/docs/VPythonDocs/graph.html
---
###### tags"`VPython`
Sign in with Wallet
Connect another wallet