重力及簡諧

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

這是在高二下的課程當中一定會出現卻又很抽象的題目：

在外太空有兩個質量為

$M$ 的星球，星球質量均勻分布且位置固定，兩星球的球星距離為

$2d$，在中垂線上距離

$x$ 處有一個質量為

$m$ 的質點。若質點原為靜止，若只考慮重力的作用，當

$x\ll d$ 時，求

$m$ 運動的週期為何？

解析：

先對

$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}}$$

m 的力圖

m 所受合力與距離 x 的關係圖 (0 ≤ x ≤ 100)

m 所受合力與距離 x 的關係圖 (0 ≤ x ≤ 8)

通常還會有加強版，將兩個質量為

$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}}$$

m 與 M 圓環示意圖

這次的目標就是要將這兩個題目畫出來。

程式 11-1：重力造成的簡諧運動, 初速為0, 從端點出發（取得程式碼）　（GlowScript 網站動畫連結）

"""
 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

參數設定

在此設定變數為 size、m、M、d、h、G、v0、i、t、dt，用途已寫在該行的註解當中。

畫面設定

1. 由於 VPython 預設的視角是從 +z 軸方向朝原點方向看去，畫面右方為 +x 軸方向、畫面上方為 +y 軸方向，原題目中 m 在 x 軸方向上運動，在動畫中改為在 y 軸方向上運動。

2. 產生星球 s1、s2 及小球 ball，並設定小球初位置為 (0, h, 0)、初速度為 0。如果之後想要試著讓小球從不同距離出發，只要修改參數設定當中的 h 即可。

3. 畫星球間的連線，也就是平衡點的位置；出發時的高度為上端點，因此下端點應該位於 y = -h 處。

4. 產生表示速度、加速度的箭頭，為了使箭頭不要疊在一起，分別畫在小球右側距離為 1、2 兩處。

5. 開啟繪圖視窗，畫小球的位置、速度、加速度與時間關係圖。

物體運動

1. 為了使動畫在小球來回 5 次時，將 while 迴圈當中的條件設定為 i < 5，並用 if(ball.pos.y > h) 判斷小球是否回到出發點，若回到出發點則印出經過的時間 t 及回到出發點的次數 i 。

2. 利用萬有引力定律計算小球所受合力 F，再由 F = ma 計算小球的加速度，更新速度及位置。

3. 最後畫 y - t、v - t、a - t 圖、更新時間。

模擬結果

若以程式當中設定數值為例，當 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 時圖形就很奇怪了。

程式 11-1 畫面截圖

程式 11-1：h = 3 的 y - t、v - t、a - t 圖

程式 11-1：h = 9 的 y - t、v - t、a - t 圖

程式 11-2.重力造成的簡諧運動, 初速不為0, 從平衡點出發（取得程式碼）　（GlowScript 網站動畫連結）

"""
 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

這個程式與 11-1 幾乎一樣，差別在於設定參數時 h = 0、v0 = 5，由於無法猜出上、下端點在何處，所以沒有先畫出這兩個位置。由於設定的 v0 不大，小球只在連心線附近運動，週期為 4.060999999999691，與理論值 3.84669 相當接近， y - t、v - t、a - t 圖也相當接近 cos、-sin、-cos 的樣子。

程式 11-2 畫面截圖

程式 11-2：v0 = 5 的 y - t、v - t、a - t 圖

程式 11-3.重力造成的簡諧運動, 圓環, 初速為0, 從端點出發（取得程式碼）　（GlowScript 網站動畫連結）

"""
 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

程式設計部分

這是以程式 11-1 為基礎修改而成的，不同之處共有以下幾點：

1. 在畫面設定部分，將 s1、s2 改為圓環，並將圓環分割為 N 個小球，每個小球的質量為 M/N，利用 for 迴圈將 N 個小球平均地放在圓環上，將資料儲存在名為 stars 的串列當中。

2. 在物體運動部分，利用 for 迴圈將 stars 中的小球一次取一個出來並命名為 star，計star 與質量為 m 的小球 ball 之間的萬有引力 f，再由 F += f 計算 ball 合力 F。

模擬結果

若以程式當中設定數值為例，當 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 時圖形就很奇怪了。

程式 11-3 畫面截圖

程式 11-3：h = 5 的 y - t、v - t、a - t 圖

程式 11-3：h = 9 的 y - t、v - t、a - t 圖

結語

同學們可以試著調整參數，模擬各種比較極端的情境，例如h很小或是h很大，看看模擬的結果是否符合課本所寫的內容。

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