一個地球人在臺北: 使用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.070 x - 0.00058 x^2$$

但這畢竟是猜測的結果，不一定正確。比較好的作法是取 $\log t$ 、$\log h$，畫 $\log t - \log h$ 關係圖，畫出來的圖為斜直線，最接近直線為

$$ y = 0.500 x - 0.345 $$

由斜率可以猜測 $t \propto \sqrt h$。理論上

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

$$ \log t = 0.5 \log h + c' $$

為了驗證這個猜測，新增一欄數據計算 $\sqrt h$，畫 $t - \sqrt h$ 關係圖，畫出來的圖為斜直線，最接近直線為

$$ y = 0.45178 x - 0.00034 $$

理論上

$$ h = \frac{1}{2}gt^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/9/6
 作者: 王一哲
"""
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_h, 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)

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

"""
 3. 用for迴圈改變h, 計算t、R, 寫入檔案
"""
for h in range(5, 51, 1):
    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(500)
        ball.v += ball.a*dt
        ball.pos += ball.v*dt
        t += dt
    print(h, t, ball.pos.x + L/2)
    file.write(str(h) + "," + str(t) + "," + str(ball.pos.x + L/2) + "\n")
    data_h.append(h)
    data_x.append(ball.pos.x + L/2)

file.close() # 關閉檔案

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

程式 5-5 與 5-4 幾乎一模一樣，以下只說明不同之處。

為了將資料儲存到文字檔中，需要加上以下的程式碼，第一行程式碼用來開啟名稱為 data.csv 的文字檔，如果檔案已存在則覆蓋檔案內容，文字編碼為 UTF-8。第二行程式碼會在 data.csv 中寫入文字 h(m), t(s), R(m) 並換行。
```
file = open("data.csv", "w", encoding = "UTF-8")
file.write("h(m), t(s), R(m)\n")
```
為了讓程式自動改變 h 的數值並代入後續的程式碼中運算，需要用到 for 迴圈。h 的數值由 5 開始代入，每次增加 1，直到數量為 50 為止。為了計算每個小球的飛行時間，每次代入 h 的數值之後，將時間 t 重設為 0。

for h in range(5, 51, 1):
    t = 0
    ball = sphere(pos = vector(-L/2, h, 0), radius = size, texture = textures.wood, make_trail = True)

物體運動部分，也就是 while 迴圈的部分，必須整個放在 for 迴圈裡面。在 while 迴圈執行完之後，使用以下程式碼印出高度 h、飛行時間 t、水平射程 R = ball.pos.x + L/2，再將 h、t、R 轉為字串、用逗號分隔、寫入檔案。。
```
print(h, t, ball.pos.x + L/2)
file.write(str(h) + "," + str(t) + "," + str(ball.pos.x + L/2) + "\n")
```
最後一定要加上 file.close() 關閉檔案。

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

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

數據處理原則

仿照程式 5-4 的方法，先畫 t - h 關係圖，多項式擬合的結果為 $$ y = 0.728 + 0.074 x - 0.0005 x^2 $$ 再畫 $ \log t - \log h $ 關係圖，畫出來的圖為斜直線，最接近直線為 $$ y = 0.500 x - 0.345 $$ 由斜率可以猜測 $ t \propto \sqrt h $。接著新增一欄數據計算 $ \sqrt h $，畫 $ t - \sqrt h $ 關係圖，畫出來的圖為斜直線，最接近直線為 $$ y = 0.45179 x - 0.00043 $$ 斜率的理論值為 0.45175，理論值與模擬得到的值很接近，可以確定程式應該沒有問題。