# 簡諧運動

> 作者：王一哲
> 日期：2018/4/5

<br />

在水平光滑桌面上有一個木塊質量為 m，用一條彈性係數為 k 的彈簧連接到左側的牆壁上，若將木塊向右拉一段距離 R 再由靜止釋放，木塊所受合力與加速度的關係為

$$ F = -kx = ma ~\Rightarrow -kx = m \frac{d^2 x}{dt^2} $$

此時木塊的運動方式稱為簡諧運動(simple harmonic motion, S.H.M.)，由上式可以解出

$$ x(t) = R \cos (\omega t + \phi) $$

$$ v(t) = -\omega R \sin (\omega t + \phi) $$

$$ a(t) = -\omega^2 R\cos (\omega t + \phi) $$

上式中 $\omega$ 為角頻率

$$ \omega = \sqrt{\frac{k}{m}} $$

週期為

$$ T = 2\pi \sqrt{\frac{m}{k}} $$

理論上我們只要在 VPython 中設定好木塊所受的彈簧回復力與木塊離開平衡點位置的關係，應該就能畫出簡諧運動的過程與週期。以下共有2個程式：

1. 理想的簡諧運動。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/8-1SHM)）
2. 考慮阻尼(damping)的簡諧運動。　（[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/8-2SHMdamp)）

<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://i.imgur.com/3HwfNHG.png">
<div style="text-align:center">簡諧運動畫面截圖</div>
<br />

## 程式 8-1：簡諧運動 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/08.%E7%B0%A1%E8%AB%A7%E9%81%8B%E5%8B%95/8-1_SHM.py)）

"""
 VPython教學: 8-1.簡諧運動
 Ver. 1: 2018/2/25
 Ver. 2: 2019/9/8
 Ver. 3: 2022/4/21 加上週期理論值, 修改 print 的格式
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m = 4               # 木塊質量 4 kg
size = 1            # 木塊邊長 1 m
R = 5               # 振幅 5 m
k = 1               # 彈性常數 1 N/m
L0 = R + size       # 彈簧原長
i = 0               # 木塊回到初位置的次數
t = 0               # 時間
dt = 0.001          # 時間間隔
print("週期理論值為 {:f} s".format(2*pi*sqrt(m/k)))

"""
 2. 畫面設定
"""
# 產生動畫視窗、地板、木塊、彈簧
scene = canvas(title="Simple Harmonic Motion", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -(size+0.1)/2, 0), size=vec(2*L0, 0.1, R), texture=textures.metal)
wall = box(pos=vec(-L0, 0, 0), size=vec(0.1, size, R), texture=textures.metal)
block = box(pos=vec(R+size/2, 0, 0), size=vec(size, size, size), texture=textures.wood, v=vec(0, 0, 0))
spring = helix(pos=vec(-L0, 0, 0), radius=0.3*size, thickness=0.05*size, color=color.yellow)
spring.axis = block.pos - spring.pos - vec(size/2, 0, 0)
# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=block.pos + vec(0, size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_a = arrow(pos=block.pos + vec(0, 2*size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=400, xtitle="<i>t</i>(s)", 
           ytitle="blue: <i>x</i>(m), green: <i>v</i>(m/s), magenta: <i>a</i>(m/s<sup>2</sup>)")
xt = gcurve(graph=gd, color=color.blue)
vt = gcurve(graph=gd, color=color.green)
at = gcurve(graph=gd, color=color.magenta)

"""
 3. 物體運動部分, 重複3個週期
"""
while i < 3:
    rate(1000)
# 計算彈簧長度、伸長量、回復力
    spring.axis = block.pos - spring.pos - vec(0.5*size, 0, 0)
    F = -k * (spring.axis - vec(L0, 0, 0))
# 計算木塊加速度, 更新速度、位置
    block.a = F/m
    block.v += block.a*dt
    block.pos += block.v*dt
# 更新代表速度、加速度的箭頭位置、方向、長度
    arrow_v.pos = block.pos + vec(0, size, 0)
    arrow_a.pos = block.pos + vec(0, 2*size, 0)
    arrow_v.axis = block.v
    arrow_a.axis = block.a
# 畫出 x-t, v-t, a-t 圖
    xt.plot(pos=(t, block.pos.x - size/2))
    vt.plot(pos=(t, block.v.x))
    at.plot(pos=(t, block.a.x))
# 判斷木塊是否回到出發點, 計算回到出發點的次數
    if block.pos.x > R + size/2 and block.v.x > 0:
        i += 1
        print("第 {:d} 次回到右端點，經過時間為 {:f} s".format(i, t))
# 更新時間
    t += dt

<br />

### 參數設定

在此定義的變數有 m、size、R、k、L0、i、t、dt，用途都已經寫在該行的註解中。
<br />

### 畫面設定
1. 由於我希望木塊的中心位於 x 軸上，因此將地板的位置向下移動到 y = -(size+0.1)/2 處。
2. 為了讓木塊在 $-R \leq x \leq R$ 之間來回運動，因此將牆壁向左移動到 x = -L0 處。
3. 由於木塊的位置是以中心點為準，而且彈簧連接在木塊左，因此需要將木塊放在彈簧右端點的右側 size/2 處。
4. 為了畫出彈簧需要用到一個新的物件 helix（螺旋線）[[3](https://www.glowscript.org/docs/VPythonDocs/helix.html)]，可以調整的選項主要有
-- 位置(pos)，以螺旋線的端點為基準。
-- 半徑(radius)
-- 軸(axis)，以位置為起點，格式為向量。
-- 厚度(thickness)，線條本身的半徑，預設值為 radius/20。
-- 顏色(color)

<br />

### 物體運動

1. 先用 spring.axis = block.pos - spring.pos - vector(size/2, 0, 0) 計算並更新彈簧的軸方向及長度，再用 F = -k * (spring.axis - vector(L0, 0, 0)) 計算彈簧的伸長量、回復力。
2. 判斷木塊是否回到出發點的部分有很多種寫法，我是用 block.pos.x > R + size/2 且 木塊速度方向向右來判斷，若計算小球回到出發點時將次數 i 加 1 並印出時間。

<br />

### 模擬結果

由下圖可以看出 x、v、a 與時間的關係符合理論預測，的確是 cos、-sin、-cos 的型式。週期模擬值為 12.565，理論值為 12.566，相當接近。

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/TUw9Eum.png">
<div style="text-align:center">m = 4、k = 1、R = 5 的簡諧運動 x - t、v - t、a - t 關係圖</div>

<br />

## 程式 8-2：簡諧運動, 有阻尼 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/08.%E7%B0%A1%E8%AB%A7%E9%81%8B%E5%8B%95/8-2_SHM_damp.py)）

"""
 VPython教學: 8-2.簡諧運動, 有阻尼
 Ver. 1: 2018/2/25
 Ver. 2: 2019/9/8
 Ver. 3: 2022/4/21 加上週期理論值, 修改 print 的格式
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
m = 4               # 木塊質量 4 kg
size = 1            # 木塊邊長 1 m
R = 5               # 振幅 5 m
k = 1               # 彈性常數 1 N/m
L0 = R + size       # 彈簧原長
b = 0.1*sqrt(4*m*k) # 阻尼 f = -bv, overdamped: b^2 > 4mk, critical damping: b^2 = 4mk, underdamped: b^2 < 4mk
T = 2*pi*sqrt(m/k)  # 週期理論值
i = 0               # 木塊運動經過的週期次數
t = 0               # 時間
dt = 0.001          # 時間間隔
print("週期理論值為 {:f} s".format(T))

"""
 2. 畫面設定
"""
# 產生動畫視窗、地板、木塊、彈簧
scene = canvas(title="Simple Harmonic Motion", width=800, height=400, x=0, y=0, background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -(size+0.1)/2, 0), size=vec(2*L0, 0.1, R), texture=textures.metal)
wall = box(pos=vec(-L0, 0, 0), size=vec(0.1, size, R), texture=textures.metal)
block = box(pos=vec(R+size/2, 0, 0), size=vec(size, size, size), texture=textures.wood, v=vec(0, 0, 0))
spring = helix(pos=vec(-L0, 0, 0), radius=0.3*size, thickness=0.05*size, color=color.yellow)
spring.axis = block.pos - spring.pos - vec(size/2, 0, 0)
# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=block.pos + vec(0, size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_a = arrow(pos=block.pos + vec(0, 2*size, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=400, xtitle="<i>t</i>(s)",
           ytitle="blue: <i>x</i>(m), green: <i>v</i>(m/s), magenta: <i>a</i>(m/s<sup>2</sup>)")
xt = gcurve(graph=gd, color=color.blue)
vt = gcurve(graph=gd, color=color.green)
at = gcurve(graph=gd, color=color.magenta)

"""
 3. 物體運動部分, 重複5個週期
"""
vp = block.v.x
while i < 5 and t < 5*T:
    rate(1000)
# 計算彈簧長度、伸長量、回復力
    spring.axis = block.pos - spring.pos - vec(0.5*size, 0, 0)
    F = -k * (spring.axis - vec(L0, 0, 0)) - b * block.v
# 計算木塊加速度, 更新速度、位置
    block.a = F/m
    block.v += block.a*dt
    block.pos += block.v*dt
# 更新代表速度、加速度的箭頭位置、方向、長度
    arrow_v.pos = block.pos + vec(0, size, 0)
    arrow_a.pos = block.pos + vec(0, 2*size, 0)
    arrow_v.axis = block.v
    arrow_a.axis = block.a
# 畫出 x-t, v-t, a-t 圖
    xt.plot(pos = (t, block.pos.x - size/2))
    vt.plot(pos = (t, block.v.x))
    at.plot(pos = (t, block.a.x))
# 判斷木塊是否經過一個週期
    vc = block.v.x
    if vp > 0 and vc < 0:
        i += 1
        print("完成第 {:d} 次振盪，經過時間為 {:f} s".format(i, t))
    vp = vc
# 更新時間
    t += dt

<br />

### 理論計算 [[6](http://hyperphysics.phy-astr.gsu.edu/hbase/oscda.html)]

假設木塊所受阻尼 $f = -bv$，由牛頓第二定律可得

$$ ma + bv + kx = 0 $$

$$ m \frac{d^2 x}{dt^2} + b \frac{dx}{dt} + kx = 0 $$

共有 3 種情形：
1. $b^2 > 4mk ~~~~~$過阻尼 (overdamped)
2. $b^2 = 4mk ~~~~~$臨界阻尼 (critical damping)
3. $b^2 < 4mk ~~~~~$次阻尼 (underdamped)

<br />

### 程式設計部分

程式 8-2 與 8-1 很像，只有 2 個不同之處。

1. 在計算木塊所受合力時需要改成 F = -k * (spring.axis - vec(L0, 0, 0)) - b * block.v

2. 由於木塊不會回到出發點，在判斷木塊是否經過一個週期時，改成用速度來判斷，如果木塊的速度原來向右、後來向左，代表木塊經過一個週期。

    vc = block.v.x
    if vp > 0 and vc < 0:
        i += 1
        print(i, t)
    vp = vc

<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/GW9ktDC.png">
<div style="text-align:center">b = 6，過阻尼</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/StiTIYn.png">
<div style="text-align:center">b = 4，臨界阻尼</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/xqNpGAj.png">
<div style="text-align:center">b = 1，次阻尼</div>
<br />

<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/0OWv3Ti.png">
<div style="text-align:center">b = 0.3，次阻尼</div>
<br />

## 參考資料
1. **canvas**: http://www.glowscript.org/docs/VPythonDocs/canvas.html
2. **box**: http://www.glowscript.org/docs/VPythonDocs/box.html
3. **helix**: http://www.glowscript.org/docs/VPythonDocs/helix.html
4. **arrow**: http://www.glowscript.org/docs/VPythonDocs/arrow.html
5. **graph**: http://www.glowscript.org/docs/VPythonDocs/graph.html
6. **Damped Harmonic Oscillator**: http://hyperphysics.phy-astr.gsu.edu/hbase/oscda.html

---

###### tags:`VPython`