<style>
.markdown-body table{
display: unset;
}
</style>
# 終端速度
> 作者：王一哲
> 第1版：2018年3月20日
> 第2版：2024年3月21日
<br />
這篇文章本來應該是前一篇文章〈[自由落下](https://hackmd.io/s/S1e8LxzGQ)〉最後一部分的內容，但由於這個程式的寫法比較不一樣，需要解釋的東西較多，所以獨立出來另外寫一篇。這個程式的目標：當小球從高空落下時，同時受到重力及空氣阻力的作用，試著找出小球的運動過程及終端速度，同時將得到的資料存成文字檔。
<br />
## 物理原理
假設空氣阻力
$$f = -bv$$
小球落下時所受合力（向下為正）
$$F = mg - bv = ma$$
小球剛開始運動時
$$ v = 0 ~~~~~ f = 0 ~~~~~ a = g $$
當小球速度 $v$ 增加時、$f$ 增加、$a$ 減小。當小球所受合力為零時，小球不會再加速，此時的速度稱為終端速度 (terminal velocity)，通常代號為 $v_t$，由上式可知
$$
v_t = \frac{mg}{b}
$$
以下圖片是以小球質量 $m = 0.1 ~\mathrm{kg}$、重力加速度 $g = 9.8 ~\mathrm{m/s^2}$、空氣阻力係數 $b = 0.1 ~\mathrm{kg/s}$ 模擬得到的結果，$v_t = 9.799001457836292 ~\mathrm{m/s}$，$理論值 = 9.8 ~\mathrm{m/s}$。
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/jZW99jv.png">
<div style="text-align:center">小球終端速度，b = 0.1， y-t 圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/XE0pzKw.png">
<div style="text-align:center">小球終端速度，b = 0.1， v-t 圖</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/6GkhYff.png">
<div style="text-align:center">小球終端速度，b = 0.1， a-t 圖</div>
<br />
## 程式 4-4：終端速度 （[取得程式碼](https://github.com/YiZheWangTw/VPythonTutorial/blob/master/04.%E8%87%AA%E7%94%B1%E8%90%BD%E4%B8%8B/4-4_airdrag_vt.py)）
```python=
"""
VPython教學: 4-4.自由落下, 有空氣阻力, 求終端速度
Ver. 1: 2018/3/8
Ver. 2: 2019/2/2
Ver. 3: 2019/9/6
Ver. 4: 2024/3/21 修改 while 迴圈的停止條件，修改參數減少程式運作時間
作者: 王一哲
"""
from vpython import *
"""
1. 參數設定, 設定變數及初始值
"""
r = 1 # 小球半徑
m = 0.1 # 小球質量
h = 0 # 小球離地高度
g = 9.8 # 重力加速度 9.8 m/s^2
b = 0.1 # 空氣阻力 f=-bv
c1, c2 = color.red, color.green
t = 0 # 時間
dt = 0.001 # 時間間隔
"""
2. 畫面設定
"""
scene = canvas(title="Terminal Velocity", width=400, height=400, x=0, y=0, center=vec(0, h/2, 0), background=vec(0, 0.6, 0.6))
# b1 with air drag, b2 without air drag
b1 = sphere(pos=vec(2*r, h, 0), radius=r, color=c1, v=vec(0, 0, 0), a=vec(0, -g, 0))
b2 = sphere(pos=vec(-2*r, h, 0), radius=r, color=c2, v=vec(0, 0, 0), a=vec(0, -g, 0))
gd = graph(title="<i>y</i> - <i>t</i> plot", width=600, height=450, x=0, y=400,
xtitle="<i>t</i> (s)", ytitle="<i>y</i> (m)", fast=False)
gd2 = graph(title="<i>v</i> - <i>t</i> plot", width=600, height=450, x=0, y=700,
xtitle="<i>t</i> (s)", ytitle="<i>v</i> (m/s)", fast=False)
gd3 = graph(title="<i>a</i> - <i>t</i> plot", width=600, height=450, x=0, y=1000,
xtitle="<i>t</i> (s)", ytitle="<i>a</i> (m/s<sup>2</sup>)", fast=False)
yt1 = gcurve(graph=gd, color=c1)
yt2 = gcurve(graph=gd, color=c2)
vt1 = gcurve(graph=gd2, color=c1)
vt2 = gcurve(graph=gd2, color=c2)
at1 = gcurve(graph=gd3, color=c1)
at2 = gcurve(graph=gd3, color=c2)
# 開啟檔案 data.csv, 屬性為寫入, 先寫入欄位的標題
file = open("data.csv", "w", encoding="UTF-8")
file.write("t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n")
tp = 0
# 設定計算終端速度用的變數
eps = 1E-6
v1 = 0
v2 = -1E6
"""
3. 物體運動部分, 小球速度變化大於 eps 時繼續運作
"""
while abs(v2 - v1) > eps:
rate(1000)
# 更新小球受力、加速度、速度、位置，畫 y-t 及 v-t 圖
f = -b*b1.v
b1.a = vec(0, -g, 0) + f/m
b1.v += b1.a*dt
b1.pos += b1.v*dt
b2.v += b2.a*dt
b2.pos += b2.v*dt
yt1.plot(pos=(t, b1.pos.y))
vt1.plot(pos=(t, b1.v.y))
at1.plot(pos=(t, b1.a.y))
yt2.plot(pos=(t, b2.pos.y))
vt2.plot(pos=(t, b2.v.y))
at2.plot(pos=(t, b2.a.y))
# 每隔 0.1 秒將資料轉成字串後寫入檔案
tc = t
if tc == 0 or tc - tp >= 0.1:
file.write(str(t) + "," + str(b1.pos.y) + "," + str(b2.pos.y) + \
"," + str(b1.v.y) + "," + str(b2.v.y) + "," + str(b1.a.y) + \
"," + str(b2.a.y) + "\n")
tp = tc
# 更新 v1 為 v2 目前的值，更新 v2 為 b1.v.y
v1, v2 = v2, b1.v.y
# 更新時間
t += dt
print("t =", t, "vt =", v2)
file.close() # 關閉檔案
```
<br />
以下只說明這個程式與程式 4-3 的不同之處。
1. 為了使用 計算小球的加速度，需要定義小球質量 m。為了計算空氣阻力 f = -bv ，需要定義係數 b。在第60、61行
```python
f = -b*b1.v
b1.a = vec(0, -g, 0) + f/m
```
計算空氣阻力 f 及加速度 b1.a。
2. 第38~43行，由於 y、v、a 的數值差異很大，為了使 v 和 a 的曲線不會被壓扁，將 y-t、v-t、a-t 分開作圖。
3. 為了將資料存成文字檔，需要增加以下的程式碼，分別位於第46 ～ 48、73 ～ 78行。第46行，先用 open 開啟檔名為 data.csv 的文字檔，**w** 代表寫入，若檔案不存在會新增檔案並寫入資料，若檔案已存在會覆寫檔案內容，最後的 **encoding = "UTF-8"** 是指定文字編碼格式為 UTF-8 (8-bit Unicode Transformation Format)。
```python
file = open("data.csv", "w", encoding="UTF-8")
```
第47行，用 **file.write** 將字串 "t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n" 寫入檔案，最後的 **\n 代表換行**。
```python
file.write("t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)\n")
```
第85行，最後一定要用 **file.close() 關閉檔案**，否則程式會執行錯誤。我們在之後的課程〈[使用 For 迴圈計算水平抛射資料](https://hackmd.io/@yizhewang/SJoo2fGfQ)〉會介紹另一個寫法，使用 with 開啟檔案，使用這個寫法就不需要關閉檔案。
```python
file.close() # 關閉檔案
```
4. 如果將所有的資料都寫入文字檔會使檔案太大，大約每隔 0.1 秒記錄一次資料應該就夠詳細了，因此在第48行定義了變數 tp、初始值為0，在第73行定義了變數 tc，將此時的 t 數值指定給 tc。接著第74 ~ 77行，當 tc 等於 0 或 tc - tp >= 0.1（經過 0.1 秒）時，用 **str** 將格式為浮點數的變數 t、b1.pos.y、……等數值轉成字串，再用 **+** 接成一個較長的字串，最後用 file.write 寫入檔案。第78行，將 tc 的數值指定給 tp。
```python
tc = t
if tc == 0 or tc - tp >= 0.1:
file.write(str(t) + "," + str(b1.pos.y) + "," + str(b2.pos.y) + \
"," + str(b1.v.y) + "," + str(b2.v.y) + "," + str(b1.a.y) + \
"," + str(b2.a.y) + "\n")
tp = tc
```
5. 如果只想要趕快算出終端速度，並不想看小球落下過程的動畫，可以刪除 while 迴圈當中的 rate(1000) ，不限制每秒執行幾次，程式能跑多快就跑多快。如果覺得刪除 rate(1000) 之後程式運算速度還是太慢，可以再刪除與繪製函數圖形有關的程式碼。
6. 為了計算終端速度，在第51 ~ 53行定義變數 eps = 1E-6、v1 = 0、v2 = -1E6。第57行，while 迴圈運作的條件為
```python
abs(v2 - v1) > eps
```
也就是速度變化仍大於我們所設定的精準度時，代表小球尚未達到終端速度，需要繼續運算。接下來在第80行，更新 v1 為 v2 目前的值，更新 v2 為 b1.v.y。
```python
v1, v2 = v2, b1.v.y
```
7. 轉存成文字檔的資料格式如下（[取得檔案](https://raw.githubusercontent.com/YiZheWangTw/VPythonTutorial/master/04.%E8%87%AA%E7%94%B1%E8%90%BD%E4%B8%8B/data.csv)）：
```csv=
t(s), y1(m), y2(m), v1(m/s), v2(m/s), a1(m/s^2), a2(m/s^2)
0,-9.800000000000001e-06,-9.800000000000001e-06,-0.009800000000000001,-0.009800000000000001,-9.8,-9.8
0.10000000000000007,-0.048837982593961625,-0.0504798,-0.9419039213273656,-0.9898000000000009,-8.86696304171435,-9.8
0.20000000000000015,-0.18632103459137758,-0.19894980000000026,-1.7852642296382608,-1.969800000000004,-8.02275852889063,-9.8
0.3000000000000002,-0.40401831622843054,-0.44541979999999975,-2.548330013785355,-2.9497999999999855,-7.258928915129776,-9.8
0.4000000000000003,-0.6942928029302021,-0.7898897999999974,-3.23874594301281,-3.9297999999999664,-6.567821878866057,-9.8
0.5000000000000003,-1.050234574833993,-1.2323597999999951,-3.8634288540200266,-4.909799999999988,-5.942513659639614,-9.8
0.6000000000000004,-1.4655915907108468,-1.7728297999999951,-4.428637046335492,-5.8898000000000135,-5.376739693357867,-9.8
0.7000000000000005,-1.9347070527533745,-2.411299799999998,-4.940032980226851,-6.869800000000039,-4.864831851624774,-9.8
0.8000000000000006,-2.452462734727637,-3.147769800000003,-5.402740005277641,-7.849800000000064,-4.40166165637874,-9.8
0.9000000000000007,-3.0142277057300206,-3.982239800000011,-5.821393687957935,-8.82980000000009,-3.9825889009430098,-9.8
...
```
可以將資料匯入 LibreOffice Calc、MicroSoft Excel 之類的試算表軟體進行後續處理。
## 結語
在這個程式當中我們學到了將資料存成文字檔的方法，匯出的資料可以再利用其它的軟體處理並作圖。雖然在 Python 當中有一個很有名的繪圖套件 matplotlib，有興趣的同學可以參考〈[Matplotlib 繪圖技巧：讀取數據及線性擬合](https://hackmd.io/s/ByN63TEC7)〉、〈[Matplotlib 繪圖技巧：在同一張圖上繪製兩個函數、兩張圖並列](https://hackmd.io/s/Bkd0EMyCm)〉。
<br />
## 2022年9月29日補充：使用 NumPy 以及 Matplotlib 繪圖
如果想要使用 Python 讀取文字檔的資料並繪圖，可以使用 NumPy 以及 Matplotlib，程式碼如下。
<br />
```python=
import numpy as np
import matplotlib.pyplot as plt
# 從 data.csv 讀取資料
t, y1, y2, v1, v2, a1, a2 = np.loadtxt('data.csv', skiprows=1, delimiter=',', unpack=True)
# 繪圖
plt.figure(figsize=(8, 5), dpi=100) # 設定圖片尺寸
plt.xlabel(r'$t ~\mathrm{(s)}$', fontsize=16) # 設定坐標軸標籤
plt.ylabel(r'$v ~\mathrm{(m/s)}$', fontsize=16)
plt.xticks(fontsize=12) # 設定坐標軸數字格式
plt.yticks(fontsize=12)
plt.grid(color='0.5', linestyle='--', linewidth=1) # 設定格線顏色、種類、寬度
# 繪圖並設定線條顏色、寬度、圖例
line1, = plt.plot(t, v1, color='red', linewidth=3, label=r'$v_1$')
line2, = plt.plot(t, v2, color='blue', linewidth=3, label=r'$v_2$')
plt.legend(handles=[line1, line2], loc='lower left', fontsize=16)
plt.savefig('v-t_plot.svg') # 儲存圖片
plt.savefig('v-t_plot.png')
plt.show() # 顯示圖片
```
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="60%" width="60%" src="https://imgur.com/bQeNGnE.png">
<div style="text-align:center">使用 Matplotlib 繪製的 v-t 圖</div>
<br /><br />
假設有以下的實驗數據，若推測 y 與 x 成反比，可以使用以下的程式碼作圖並計算最接近直線的斜率及截距。
<div style="text-align:center">
| x | 0.2 | 0.4 | 0.6 | 0.8 | 1.0 |
| - | --- | --- | --- | --- | --- |
| y | 10.1| 4.9 | 3.3 | 2.6 | 1.9 |
</div><br /><br />
```python=
import numpy as np
import matplotlib.pyplot as plt
x = np.asarray([0.2, 0.4, 0.6, 0.8, 1.0])
y = np.asarray([10.1, 4.9, 3.3, 2.6, 1.9])
z = 1/x
# 產生計算最接近直線需要用到的陣列 A ，計算直線的斜率 m 和截距 c
A = np.vstack([z, np.ones(len(z))]).T
m, c = np.linalg.lstsq(A, y, rcond = -1)[0]
print(m, c)
plt.figure(figsize=(8, 6), dpi=100) # 設定圖片尺寸
plt.xlabel(r'$1/x$', fontsize=16) # 設定坐標軸標籤
plt.ylabel(r'$y$', fontsize=16)
plt.xticks(fontsize=12) # 設定坐標軸數字格式
plt.yticks(fontsize=12)
plt.grid(color='0.5', linestyle='--', linewidth=1) # 設定格線顏色、種類、寬度
# 繪圖並設定線條顏色、寬度、圖例
data, = plt.plot(z, y, color='blue', marker='o', markersize=10, linestyle='', label="Data")
fitline, = plt.plot(z, m*z + c, color='orange', linewidth=2, label="Fitted Line")
plt.legend(handles=[data, fitline], loc='upper left', fontsize=16)
plt.savefig('fit_plot.svg') # 儲存圖片
plt.savefig('fit_plot.png')
plt.show() # 顯示圖片
```
<br />
最接近直線的斜率、縱軸截距分別為 $slope = 2.028088701161563$、$intercept = -0.07080253431890154$，以下是繪圖成果。
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/aQKfx36.png">
<div style="text-align:center">使用 Matplotlib 繪製資料點及最接近直線</div>
<br /><br />
## VPython官方說明書
1. **canvas**: http://www.glowscript.org/docs/VPythonDocs/canvas.html
2. **sphere**: http://www.glowscript.org/docs/VPythonDocs/sphere.html
3. **graph**: http://www.glowscript.org/docs/VPythonDocs/graph.html
---
###### tags: `VPython`
Sign in with Wallet
Connect another wallet