# 重力及簡諧
> 作者：王一哲
> 日期：2018/4/15
<br />
這是在高二下的課程當中一定會出現卻又很抽象的題目：
> 在外太空有兩個質量為 $M$ 的星球，星球質量均勻分布且位置固定，兩星球的球星距離為 $2d$，在中垂線上距離 $x$ 處有一個質量為 $m$ 的質點。若質點原為靜止，若只考慮重力的作用，當 $x \ll d$ 時，求 $m$ 運動的週期為何？
<img style="display: block; margin-left: auto; margin-right: auto" height="50%" width="50%"
src="https://i.imgur.com/9LzKJEj.png">
<br />
#### 解析：
先對 $m$ 畫力圖，則 $m$ 所受合力
$$ F_x = -\frac{2GMmx}{(d^2 + x^2)^{\frac{3}{2}}} \approx -\frac{2GMm}{d^3} x = -kx $$
$$ T = 2\pi \sqrt{\frac{m}{k}} = 2\pi \sqrt{\frac{d^3}{2GM}}$$
<img style="display: block; margin-left: auto; margin-right: auto" height="50%" width="50%"
src="https://i.imgur.com/bfaVxHw.png">
<div style="text-align:center"><i>m</i> 的力圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%"
src="https://imgur.com/CoEApAD.png">
<div style="text-align:center"><i>m</i> 所受合力與距離 <i>x</i> 的關係圖 (0 ≤ x ≤ 100)</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%"
src="https://imgur.com/b8C4qYQ.png">
<div style="text-align:center"><i>m</i> 所受合力與距離 <i>x</i> 的關係圖 (0 ≤ x ≤ 8)</div>
<br />
通常還會有加強版，將兩個質量為 $m$ 的星球換成一個質量均勻分布、總質量為 $m$ 、半徑為 $r$ 的圓環。若其餘的條件相同，則 $m$ 所受合力及運動週期：
$$ F_x = -\frac{GMmx}{(r^2 + x^2)^{\frac{3}{2}}} \approx -\frac{GMm}{r^3} x = -kx $$
$$ T = 2\pi \sqrt{\frac{m}{k}} = 2\pi \sqrt{\frac{r^3}{GM}}$$
<img style="display: block; margin-left: auto; margin-right: auto" height="50%" width="50%"
src="https://i.imgur.com/bwn61kS.png">
<div style="text-align:center"> <i>m</i> 與 <i>M</i> 圓環示意圖</div>
<br />
這次的目標就是要將這兩個題目畫出來。
<br />
## 程式 11-1：重力造成的簡諧運動, 初速為0, 從端點出發（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/11.%E9%87%8D%E5%8A%9B%E5%8F%8A%E7%B0%A1%E8%AB%A7/11-1_gravity_SHM.py)） （[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/11-1gravitySHM)）
```python=
"""
VPython教學: 11-1.重力造成的簡諧運動, 初速為0, 從端點出發
Ver. 1: 2018/2/23
Ver. 2: 2018/3/14
Ver. 3: 2019/9/8
作者: 王一哲
"""
from vpython import *
"""
1. 參數設定, 設定變數及初始值, 當 h << d 時週期理論值 3.84669 s
"""
size = 1 # 小球半徑
m = 1 # 小球質量
M = 2E13 # 星球質量
d = 10 # 星球之間的距離*0.5倍
h = 3 # 小球在星球連線的中垂線上的距離
G = 6.67E-11 # 重力常數
v0 = 0 # 小球初速
i = 0 # 小球回到初位置的次數
t = 0 # 時間
dt = 0.001 # 時間間隔
T = 2*pi*sqrt(d**3/(2*G*M)) # 週期理論值
print("週期理論值 = {:f} s".format(T))
"""
2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="Gravity and SHM", width=600, height=600, x=0, y=0, center=vec(0, 0, 0),
background=color.black, range=1.3*d)
# 產生 2 個質量為 M 的星球
s1 = sphere(pos=vec(-d, 0, 0), radius=size, color=color.blue)
s2 = sphere(pos=vec(d, 0, 0), radius=size, color=color.blue)
# 產生質量為 m 的星球並設定初速度
ball = sphere(pos=vec(0, h, 0), radius=0.4*size, v=vec(0, v0, 0), color=color.red)
# 畫星球間的連線、上端點、下端點
line = cylinder(pos=s1.pos, axis=s2.pos-s1.pos, radius=0.1*size, color=color.yellow)
top = cylinder(pos=vec(-2, h, 0), axis=vec(4, 0, 0), radius=0.1*size, color=color.white)
bottom = cylinder(pos=vec(-2, -h, 0), axis=vec(4, 0, 0), radius=0.1*size, color=color.white)
# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=ball.pos+vec(1, 0, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_a = arrow(pos=ball.pos+vec(2, 0, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)",
ytitle="blue: <i>y</i> (m), green: <i>v</i> (m/s), magenta: <i>a</i> (m/s<sup>2</sup>)")
yt = gcurve(graph=gd, color=color.blue)
vt = gcurve(graph=gd, color=color.green)
at = gcurve(graph=gd, color=color.magenta)
"""
3. 物體運動部分, 小球來回 5 次時停止
"""
while i < 5:
rate(1000)
# 計算小球所受重力並存到變數 F
r1 = ball.pos-s1.pos
r2 = ball.pos-s2.pos
F1 = -(G*M*m)/r1.mag2*r1.norm()
F2 = -(G*M*m)/r2.mag2*r2.norm()
F = F1+F2
# 計算運動中小球的加速度、速度、位置, 畫出代表速度、加速度的箭頭
ball.a = F/m
ball.v += ball.a*dt
ball.pos += ball.v*dt
arrow_v.pos = ball.pos+vec(1, 0, 0)
arrow_a.pos = ball.pos+vec(2, 0, 0)
arrow_v.axis = ball.v
arrow_a.axis = ball.a
# 畫出 y-t, v-t, a-t 圖
yt.plot(pos=(t, ball.pos.y))
vt.plot(pos=(t, ball.v.y))
at.plot(pos=(t, ball.a.y))
# 判斷小球是否回到出發點, 計算回到出發點的次數
if ball.pos.y > h:
i += 1
print("第 {:d} 次回到出發點, t = {:f} s".format(i, t))
# 更新時間
t += dt
```
<br />
### 參數設定
在此設定變數為 size、m、M、d、h、G、v0、i、t、dt，用途已寫在該行的註解當中。
<br />
### 畫面設定
1. 由於 VPython 預設的視角是從 +z 軸方向朝原點方向看去，畫面右方為 +x 軸方向、畫面上方為 +y 軸方向，原題目中 m 在 x 軸方向上運動，在動畫中改為在 y 軸方向上運動。
2. 產生星球 s1、s2 及小球 ball，並設定小球初位置為 (0, h, 0)、初速度為 0。如果之後想要試著讓小球從不同距離出發，只要修改參數設定當中的 h 即可。
3. 畫星球間的連線，也就是平衡點的位置；出發時的高度為上端點，因此下端點應該位於 y = -h 處。
4. 產生表示速度、加速度的箭頭，為了使箭頭不要疊在一起，分別畫在小球右側距離為 1、2 兩處。
5. 開啟繪圖視窗，畫小球的位置、速度、加速度與時間關係圖。
<br />
### 物體運動
1. 為了使動畫在小球來回 5 次時，將 while 迴圈當中的條件設定為 i < 5，並用 **if(ball.pos.y > h)** 判斷小球是否回到出發點，若回到出發點則印出經過的時間 t 及回到出發點的次數 i 。
2. 利用萬有引力定律計算小球所受合力 F，再由 F = ma 計算小球的加速度，更新速度及位置。
3. 最後畫 y - t、v - t、a - t 圖、更新時間。
<br />
### 模擬結果
若以程式當中設定數值為例，當 h 很小的時候，週期理論值為 3.84669；若 h = 0.00001 則週期為 3.8449999999996876；若 h = 3 則週期為 4.036999999999683；若 h = 9 則週期為 5.439000000000151。若從 y - t、v - t、a - t 圖來看，當 h = 3 時圖形相當接近 cos、-sin、-cos 的樣子，但是當 h = 9 時圖形就很奇怪了。
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/BuSmtSX.png">
<div style="text-align:center">程式 11-1 畫面截圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/JGl1CLR.png">
<div style="text-align:center">程式 11-1：h = 3 的 y - t、v - t、a - t 圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/4TjLGTH.png">
<div style="text-align:center">程式 11-1：h = 9 的 y - t、v - t、a - t 圖</div>
<br />
## 程式 11-2.重力造成的簡諧運動, 初速不為0, 從平衡點出發（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/11.%E9%87%8D%E5%8A%9B%E5%8F%8A%E7%B0%A1%E8%AB%A7/11-2_gravity_SHM_2.py)） （[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/11-2gravitySHM2)）
```python=
"""
VPython教學: 11-2.重力造成的簡諧運動, 初速不為0, 從平衡點出發
Ver. 1: 2018/2/23
Ver. 2: 2018/3/14
Ver. 3: 2019/9/8
作者: 王一哲
"""
from vpython import *
"""
1. 參數設定, 設定變數及初始值, 當 h << d 時週期理論值 3.84669 s
"""
size = 1 # 小球半徑
m = 1 # 小球質量
M = 2E13 # 星球質量
d = 10 # 星球之間的距離*0.5倍
h = 0 # 小球在星球連線的中垂線上的距離
G = 6.67E-11 # 重力常數
v0 = 5 # 小球初速
i = 0 # 小球回到初位置的次數
t = 0 # 時間
dt = 0.001 # 時間間隔
T = 2*pi*sqrt(d**3/(2*G*M)) # 週期理論值
print("週期理論值 = {:f} s".format(T))
"""
2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="Gravity and SHM", width=600, height=600, x=0, y=0,
center=vec(0, 0, 0), range=1.6*d, background=color.black)
# 產生 2 個質量為 M 的星球
s1 = sphere(pos=vec(-d, 0, 0), radius=size, color=color.blue)
s2 = sphere(pos=vec(d, 0, 0), radius=size, color=color.blue)
# 產生質量為 m 的星球並設定初速度
ball = sphere(pos=vec(0, h, 0), radius=0.4*size, color=color.red, v=vec(0, v0, 0))
# 畫星球間的連線
line = cylinder(pos=s1.pos, axis=s2.pos-s1.pos, radius=0.1*size, color=color.yellow)
# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=ball.pos+vec(1, 0, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.green)
arrow_a = arrow(pos=ball.pos+vec(2, 0, 0), axis=vec(0, 0, 0), shaftwidth=0.3*size, color=color.magenta)
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)",
ytitle="blue: <i>y</i> (m), green: <i>v</i> (m/s), magenta: <i>a</i> (m/s<sup>2</sup>)")
yt = gcurve(graph=gd, color=color.blue)
vt = gcurve(graph=gd, color=color.green)
at = gcurve(graph=gd, color=color.magenta)
"""
3. 物體運動部分, 小球來回 5 次時停止
"""
while i < 6:
rate(1000)
# 計算小球所受重力並存到變數 F
r1 = ball.pos-s1.pos
r2 = ball.pos-s2.pos
F1 = -(G*M*m) / r1.mag2 * r1.norm()
F2 = -(G*M*m) / r2.mag2 * r2.norm()
F = F1+F2
# 計算運動中小球的加速度、速度、位置, 畫出代表速度、加速度的箭頭
ball.a = F/m
ball.v += ball.a*dt
ball.pos += ball.v*dt
arrow_v.pos = ball.pos+vec(1, 0, 0)
arrow_a.pos = ball.pos+vec(2, 0, 0)
arrow_v.axis = ball.v
arrow_a.axis = ball.a
# 畫出 y-t, v-t, a-t 圖
yt.plot(pos=(t, ball.pos.y))
vt.plot(pos=(t, ball.v.y))
at.plot(pos=(t, ball.a.y))
# 判斷小球是否回到出發點, 計算回到出發點的次數, 出發時會多算1次
if ball.pos.y >= 0 and ball.v.y >= v0:
if i >= 1: print("第 {:d} 次回到出發點, t = {:f} s".format(i, t))
i += 1
# 更新時間
t += dt
```
<br />
這個程式與 11-1 幾乎一樣，差別在於設定參數時 h = 0、v0 = 5，由於無法猜出上、下端點在何處，所以沒有先畫出這兩個位置。由於設定的 v0 不大，小球只在連心線附近運動，週期為 4.060999999999691，與理論值 3.84669 相當接近， y - t、v - t、a - t 圖也相當接近 cos、-sin、-cos 的樣子。
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/d3S3ysj.png">
<div style="text-align:center">程式 11-2 畫面截圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/bKhMaFI.png">
<div style="text-align:center">程式 11-2：v0 = 5 的 y - t、v - t、a - t 圖</div>
<br />
## 程式 11-3.重力造成的簡諧運動, 圓環, 初速為0, 從端點出發（[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/11.%E9%87%8D%E5%8A%9B%E5%8F%8A%E7%B0%A1%E8%AB%A7/11-3_gravity_SHM_ring.py)） （[GlowScript 網站動畫連結](http://www.glowscript.org/#/user/yizhe/folder/Public/program/11-3gravitySHMring)）
```python=
"""
VPython教學: 11-3.重力造成的簡諧運動, 圓環, 初速為0, 從端點出發
Ver. 1: 2018/2/25
Ver. 2: 2018/3/14
Ver. 3: 2019/9/8
作者: 王一哲
"""
from vpython import *
"""
1. 參數設定, 設定變數及初始值, 當 h << r 時週期理論值 3.44058 s
"""
size = 0.4 # 小球半徑
m = 1 # 小球質量
M = 5E13 # 圓環的質量
r = 10 # 圓環的半徑
N = 100 # 圓環分割成小球的數量
dM = M/N # 圓環分割後每個小球的質量
h = 9 # 小球在圓環中垂線上的距離
G = 6.67E-11 # 重力常數
v0 = 0 # 小球初速
j = 0 # 小球回到初位置的次數
t = 0 # 時間
dt = 0.001 # 時間間隔
T = 2*pi*sqrt(r**3/(G*M)) # 週期理論值
print("週期理論值 = {:f} s".format(T))
"""
2. 畫面設定
"""
# 產生動畫視窗
scene = canvas(title="Gravity and SHM", width=600, height=600, x=0, y=0,
center=vec(0, 0.3*h, 0), range=1.6*r, background=color.black)
# 產生質量為 m 的星球並設定初速度
ball = sphere(pos=vec(0, h, 0), radius=size, color=color.red, v=vec(0, v0, 0))
# 畫圓環間的連線、上端點、下端點
line1 = cylinder(pos=vec(-r, 0, 0), axis=vec(2*r, 0, 0), radius=0.1*size, color=color.yellow)
line2 = cylinder(pos=vec(0, 0, -r), axis=vec(0, 0, 2*r), radius=0.1*size, color=color.yellow)
top = cylinder(pos=vec(-2, h, 0), axis=vec(4, 0, 0), radius=0.1*size, color=color.white)
bottom = cylinder(pos=vec(-2, -h, 0), axis=vec(4, 0, 0), radius=0.1*size, color=color.white)
# 產生表示速度、加速度的箭頭
arrow_v = arrow(pos=ball.pos+vec(1, 0, 0), axis=vec(0, 0, 0), shaftwidth=0.8*size, color=color.green)
arrow_a = arrow(pos=ball.pos+vec(2, 0, 0), axis=vec(0, 0, 0), shaftwidth=0.8*size, color=color.magenta)
# 產生質量為 M 的圓環及空白 list, 用 for 迴圈產生圓環分割後的小球並填入空白 list 中
stars_ring = ring(pos=vec(0, 0, 0), axis=vec(0, 1, 0), radius=r, thickness=0.4*size, color=color.green)
stars = []
for i in range(0, N):
stars.append(sphere(pos=vec(r*cos(i*2*pi/N), 0, r*sin(i*2*pi/N)), radius=0.5*size, color=color.blue))
# 繪圖部分
gd = graph(title="plot", width=600, height=450, x=0, y=600, xtitle="<i>t</i> (s)",
ytitle="blue: <i>y</i> (m), green: <i>v</i> (m/s), magenta: <i>a</i> (m/s<sup>2</sup>)")
yt = gcurve(graph=gd, color=color.blue)
vt = gcurve(graph=gd, color=color.green)
at = gcurve(graph=gd, color=color.magenta)
"""
3. 物體運動部分, 小球來回 5 次時停止
"""
while j < 5:
rate(1000)
# 用 for 迴圈將 stars 當中的小球依序讀取出來指定給變數 star, 計算小球的重力並存到變數 F
F = vec(0, 0, 0)
for star in stars:
r = ball.pos-star.pos
f = -(G*dM*m)/r.mag2*r.norm()
F += f
# 計算運動中小球的加速度、速度、位置, 畫出代表速度、加速度的箭頭
ball.a = F/m
ball.v += ball.a*dt
ball.pos += ball.v*dt
arrow_v.pos = ball.pos+vec(1, 0, 0)
arrow_a.pos = ball.pos+vec(2, 0, 0)
arrow_v.axis = ball.v
arrow_a.axis = ball.a
# 畫出 y-t, v-t, a-t 圖
yt.plot(pos=(t, ball.pos.y))
vt.plot(pos=(t, ball.v.y))
at.plot(pos=(t, ball.a.y))
# 判斷小球是否回到出發點, 計算回到出發點的次數
if ball.pos.y >= h and ball.v.y >= 0:
j += 1
print("第 {:d} 次回到出發點, t = {:f} s".format(j, t))
# 更新時間
t += dt
```
<br />
### 程式設計部分
這是以程式 11-1 為基礎修改而成的，不同之處共有以下幾點：
1. 在畫面設定部分，將 s1、s2 改為圓環，並將圓環分割為 N 個小球，每個小球的質量為 M/N，利用 for 迴圈將 N 個小球平均地放在圓環上，將資料儲存在名為 stars 的串列當中。
2. 在物體運動部分，利用 for 迴圈將 stars 中的小球一次取一個出來並命名為 star，計star 與質量為 m 的小球 ball 之間的萬有引力 f，再由 F += f 計算 ball 合力 F。
<br />
### 模擬結果
若以程式當中設定數值為例，當 h 很小的時候，週期理論值為 3.44058；若 h = 0.00001 則週期為 3.4389999999997323；若 h = 5 則週期為 3.9069999999996807；若 h = 9 則週期為 4.864999999999959。若從 y - t、v - t、a - t 圖來看，當 h = 5 時圖形相當接近 cos、-sin、-cos 的樣子，但是當 h = 9 時圖形就很奇怪了。
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://i.imgur.com/VkPbAEY.png">
<div style="text-align:center">程式 11-3 畫面截圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/c6v2Eht.png">
<div style="text-align:center">程式 11-3：h = 5 的 y - t、v - t、a - t 圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/TO51ZiV.png">
<div style="text-align:center">程式 11-3：h = 9 的 y - t、v - t、a - t 圖</div>
<br />
## 結語
同學們可以試著調整參數，模擬各種比較極端的情境，例如h很小或是h很大，看看模擬的結果是否符合課本所寫的內容。
<br />
## VPython官方說明書
1. **canvas**: http://www.glowscript.org/docs/VPythonDocs/canvas.html
2. **sphere**: http://www.glowscript.org/docs/VPythonDocs/sphere.html
3. **arrow**: http://www.glowscript.org/docs/VPythonDocs/arrow.html
4. **ring**: http://www.glowscript.org/docs/VPythonDocs/ring.html
5. **graph**: http://www.glowscript.org/docs/VPythonDocs/graph.html
---
###### tags: `VPython`
Sign in with Wallet
Connect another wallet