使用 For 迴圈計算水平抛射資料

作者：王一哲
日期：2018/3/28

在前一篇文章〈水平抛射〉中，我們寫了一個模擬水平抛射運動的程式，再手動修高度 h 計算飛行時間 t 及水平射程 R。但是手動修改是很麻煩的，如果能多寫一個 for 迴圈，就能讓程式幫我們自動代入不同的 h 並將資料存成文字檔。以下共有兩個程式：

• 程式5-4：水平抛射, 改變h, 記錄R　（GlowScript 網站動畫連結）
• 程式5-5：由程式5-4修改, 用for迴圈改變h, 記錄t、R

程式 5-4：水平抛射, 改變h, 記錄R （取得程式碼）

"""
 VPython教學: 5-4.水平抛射, 改變h, 記錄R
 Ver. 1: 2018/2/19
 Ver. 2: 2019/9/6
 作者: 王一哲
"""
from vpython import *

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1     # 小球半徑
v0 = 5       # 小球水平初速
h = 15       # 小球離地高度5, 10, 15, 20, 25, 30, 35, 40, 45, 50
L = 50       # 地板長度
g = 9.8      # 重力加速度 9.8 m/s^2
t = 0        # 時間
dt = 0.001   # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Projectile", width=800, height=600, x=0, y=0, center=vec(0, h/2, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)
ball = sphere(pos=vec(-L/2, h, 0), radius=size, texture=textures.wood, make_trail=True, v=vec(v0, 0, 0), a=vec(0, -g, 0))

"""
 3. 物體運動部分, 小球觸地時停止
"""
while ball.pos.y - floor.pos.y > size + 0.5*floor.height:
    rate(1000)
    ball.v += ball.a*dt
    ball.pos += ball.v*dt
    t += dt

print(t, ball.pos.x + L/2)

這個程式是程式 5-1 的簡化版，由於我只想算出飛行時間 t 及水平射程 R，小球碰到地板後就可以停止動畫，因此 while 迴圈中設定的條件為

ball.pos.y - floor.pos.y > size + 0.5*floor.height

手動修改 h，分別以 5、10、15、……50代入，再將計算得到的 t 和 R 手動輸入 SciDAVis 並作圖，這個作法雖然可行，但是相當費時。

數據處理原則

通常我們只用 XY 散布圖，而且只畫數據點，不會將數據點連線，畫完數據點之後再加上擬合直線式曲線。如果應變數與自變數的關係圖為曲線，還要想辦法化為直線，最後的成果一定要用線性擬合，還要解釋最接近直線的斜率、y軸截距的物理意義，要求相當嚴謹。使用 SciDAVis 畫 XY 散布圖的方法請參考另一篇文章〈SciDAVis繪圖：XY散布圖〉。

先畫 t - h 關係圖，由於圖形看起來像是二次曲線，就試著用多項式擬合，最高次方為二次，結果為

$$y=0.656+0.070x-0.00058x^2$$

但這畢竟是猜測的結果，不一定正確。比較好的作法是取

$\log t$ 、

$\log h$，畫

$\log t-\log h$ 關係圖，畫出來的圖為斜直線，最接近直線為

$$y=0.500x-0.345$$

由斜率可以猜測

$t\propto \sqrt{h}$。理論上

$$h=\frac{1}{2}g{t}^2\Rightarrow \log h=\log \left(\frac{1}{2}g{t}^2\right)=\log \left(\frac{1}{2}g\right)+\log t^2=2\log t+c$$

$$\log t=0.5\log h+c^{\prime }$$

為了驗證這個猜測，新增一欄數據計算

$\sqrt{h}$，畫

$t-\sqrt{h}$ 關係圖，畫出來的圖為斜直線，最接近直線為

$$y=0.45178x-0.00034$$

理論上

$$h=\frac{1}{2}g{t}^2\Rightarrow t=\sqrt{\frac{2h}{g}}\Rightarrow slope=\sqrt{\frac{2}{g}}\approx 0.45175$$

理論值與模擬得到的值很接近，可以確定程式應該沒有問題。

於 SciDAVis 手動輸入數據

t - h 關係圖

log(t) - log(h) 關係圖

t - h^1/2 關係圖

程式 5-5：水平抛射, 用for迴圈改變h, 記錄t、R （取得程式碼）

"""
 VPython教學: 5-5.水平抛射, 用for迴圈改變h, 記錄t、R
 Ver. 1: 2018/3/21
 Ver. 2: 2018/6/13  加上用 matplotlib.pyplot 繪圖的部分
 Ver. 3: 2019/10/10 改用 with 開啟檔案, 上課用的版本
 作者: 王一哲
"""
from vpython import *
import matplotlib.pyplot as plt

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1     # 小球半徑
v0 = 5       # 小球水平初速
L = 50       # 地板長度
g = 9.8      # 重力加速度 9.8 m/s^2
dt = 0.001   # 時間間隔
data_y, data_x = [], []     # 儲存繪圖資料的串列

"""
 2. 畫面設定
"""
scene = canvas(title="Projectile", width=800, height=600, x=0, y=0, center=vec(0, 5, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)

# 開啟檔案 5-5_data.csv, 屬性為寫入, 先寫入欄位的標題
with open("5-5_data.csv", "w", encoding="UTF-8") as file:
    file.write("h(m), t(s), R(m)\n")

"""
 3. 主程式
"""
def main(h):
    t = 0
    ball = sphere(pos=vec(-L/2, h, 0), radius=size, texture=textures.wood, make_trail=True,
                  v=vec(v0, 0, 0), a=vec(0, -g, 0))
    while ball.pos.y - floor.pos.y > size + 0.5*floor.height:
        rate(1000)
        ball.v += ball.a*dt
        ball.pos += ball.v*dt
        t += dt
    return t, ball.pos.x + L/2

"""
 4. 用for迴圈改變h, 計算t、R, 寫入檔案
"""
for h in range(5, 51, 1):
    t, r = main(h)
    with open("5-5_data.csv", "a", encoding="UTF-8") as file:
        file.write(str(h) + "," + str(t) + "," + str(r) + "\n")
    data_y.append(t)
    data_x.append(h)

# 用 matplotlib.pyplot 繪圖
plt.figure(figsize=(8, 6), dpi=100)                   # 設定圖片尺寸
plt.xlabel(r'$h ~\mathrm{(m)}$', fontsize=16)         # 設定坐標軸標籤
plt.ylabel(r'$t ~\mathrm{(s)}$', fontsize=16)
plt.xticks(fontsize=12)                               # 設定坐標軸數字格式
plt.yticks(fontsize=12)
plt.grid(color='red', linestyle='--', linewidth=1)    # 設定格線顏色、種類、寬度
plt.plot(data_x, data_y, marker='o', markerfacecolor='blue', markersize=8, color='skyblue', linewidth=3)   # 繪圖並設定線條顏色、寬度
plt.savefig('t-h_plot.svg')                           # 儲存圖片
plt.savefig('t-h_plot.png')
plt.show()                                            # 顯示圖片

程式 5-5 模擬結果畫面截圖

程式 5-5 當中運用了串列 (list) 及自訂函式 (def) 兩個新的物件，詳細的說明請參考另外兩篇講義〈串列〉、〈自訂函式〉 ，以下先說明與程式 5-4 的不同之處。

1. 為了將資料儲存到文字檔中，需要加上以下的程式碼：

   ​​​​with open("5-5_data.csv", "w", encoding="UTF-8") as file
   ​​​​    file.write("h(m), t(s), R(m)\n")
   

   第一行程式碼用來開啟名稱為 5-5_data.csv 的文字檔，如果檔案已存在則覆蓋檔案內容，文字編碼為 UTF-8。第二行程式碼會在 5-5_data.csv 中寫入文字 h(m), t(s), R(m) 並換行。

2. 第 34 ~ 42 行，用 def 自訂名為 main 的函式，輸入變數 h 數值，自動跑完小球落下的過程，回傳落下需要的時間 t 以及水平射程 ball.pos.x + L/2 。如果要讓動畫不要執行地太快，要在 while 迴圈當中加上 rate(1000)。

3. 為了讓程式自動改變 h 的數值並代入後續的程式碼中運算，需要用到 for 迴圈。

   ​​​​for h in range(5, 51, 1):
   ​​​​    t, r = main(h)
   ​​​​    with open("5-5_data.csv", "a", encoding="UTF-8") as file:
   ​​​​        file.write(str(h) + "," + str(t) + "," + str(r) + "\n")
   ​​​​    data_y.append(t)
   ​​​​    data_x.append(h)
   

   h 的數值由 5 開始代入，每次增加 1，直到數量為 50 為止。將每個小球的飛行時間及水平射程回傳為 t 及 r，再用

   ​​​​file.write(str(h) + "," + str(t) + "," + str(r) + "\n")
   

   將 h、t、r 轉為字串、用逗號分隔、寫入檔案。

4. 由於是用 with 開啟檔案，不需要加上 file.close() 關閉檔案。

5. 儲存到文字檔中的資料格式如下（取得檔案）：

h(m), t(s), R(m)
5,1.0099999999999996,5.049999999998995
6,1.105999999999989,5.5299999999989
7,1.1949999999999792,5.974999999998811
8,1.2769999999999702,6.38499999999873
9,1.3549999999999616,6.774999999998652
10,1.4279999999999535,7.1399999999985795

6. 第 55 ~ 64 行的程式碼是利用 matplotlib 套件繪圖，如果想要用其它軟體讀取 5-5_data.csv 的資料並繪圖，可以刪除這幾行。

7. 由於 GlowScript 網站上不能將資料儲存成文字檔，也不支援 matplotlib，可以使用以下的程式碼，將不支援的部分刪掉，並用 print 輸出需要的資料。程式碼第1行 Web VPython 3.2 是在網站上建立程式時自動産生的，會自動採用最新的版本，不需要手動輸入，也不需要再寫 from vpython import *，這是 GlowScript 動畫連結。

Web VPython 3.2

"""
 1. 參數設定, 設定變數及初始值
"""
size = 1     # 小球半徑
v0 = 5       # 小球水平初速
L = 50       # 地板長度
g = 9.8      # 重力加速度 9.8 m/s^2
dt = 0.001   # 時間間隔

"""
 2. 畫面設定
"""
scene = canvas(title="Projectile", width=800, height=600, x=0, y=0, center=vec(0, 5, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, -size, 0), size=vec(L, 0.01, 10), texture=textures.metal)

"""
 3. 主程式
"""
def main(h):
    t = 0
    ball = sphere(pos=vec(-L/2, h, 0), radius=size, texture=textures.wood, make_trail=True,
                  v=vec(v0, 0, 0), a=vec(0, -g, 0))
    while ball.pos.y - floor.pos.y > size + 0.5*floor.height:
        rate(1000)
        ball.v += ball.a*dt
        ball.pos += ball.v*dt
        t += dt
    return t, ball.pos.x + L/2

"""
 4. 用for迴圈改變h, 計算t、R, 寫入檔案
"""
for h in range(5, 51, 1):
    t, r = main(h)
    print("{}, {}, {}".format(h, t, r))

於 GlowScript 網站執行的結果

數據處理原則

仿照程式 5-4 的方法，先畫 t - h 關係圖，多項式擬合的結果為

$$y=0.728+0.074x-0.0005x^2$$

再畫

$\log t-\log h$ 關係圖，畫出來的圖為斜直線，最接近直線為

$$y=0.500x-0.345$$

由斜率可以猜測

$t\propto \sqrt{h}$。接著新增一欄數據計算

$\sqrt{h}$，畫

$t-\sqrt{h}$ 關係圖，畫出來的圖為斜直線，最接近直線為

$$y=0.45179x-0.00043$$

斜率的理論值為 0.45175，理論值與模擬得到的值很接近，可以確定程式應該沒有問題。

t - h 關係圖

log(t) - log(h) 關係圖

t - h^1/2 關係圖

使用 Matplotlib 繪製的 t - h 關係圖

結語

本文是以改變 h 為例，可以試著用相同的方法改變其它的變數值，例如改變空氣阻力的係數，找出飛行時間與空氣阻力係數的關係。而且用 for 迴圈自動改變數值，再將運算結果存入文字檔中，比手動改變數值所需要的時間少很多，這樣才能發揮寫程式的好處。

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. texture: http://www.glowscript.org/docs/VPythonDocs/textures.html

---

tags: VPython