Try   HackMD

終端速度

作者：王一哲
第1版：2018年3月20日
第2版：2024年3月21日


這篇文章本來應該是前一篇文章〈自由落下〉最後一部分的內容，但由於這個程式的寫法比較不一樣，需要解釋的東西較多，所以獨立出來另外寫一篇。這個程式的目標：當小球從高空落下時，同時受到重力及空氣阻力的作用，試著找出小球的運動過程及終端速度，同時將得到的資料存成文字檔。

物理原理

假設空氣阻力

f=−bv

小球落下時所受合力（向下為正）

F=mg−bv=ma

小球剛開始運動時

v=0     f=0     a=g

當小球速度

v 增加時、
f 增加、
a 減小。當小球所受合力為零時，小球不會再加速，此時的速度稱為終端速度 (terminal velocity)，通常代號為
vt，由上式可知
vt=mgb

以下圖片是以小球質量

m=0.1 kg、重力加速度
g=9.8 m/s2、空氣阻力係數
b=0.1 kg/s 模擬得到的結果，
vt=9.799001457836292 m/s，
理論值=9.8 m/s。

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
小球終端速度，b = 0.1， y-t 圖

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
小球終端速度，b = 0.1， v-t 圖

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
小球終端速度，b = 0.1， a-t 圖

程式 4-4：終端速度 （取得程式碼）






















































































"""
 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() # 關閉檔案

以下只說明這個程式與程式 4-3 的不同之處。

  1. 為了使用 計算小球的加速度，需要定義小球質量 m。為了計算空氣阻力 f = -bv ，需要定義係數 b。在第60、61行

    ​​​​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)。

    ​​​​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 代表換行。

    ​​​​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 迴圈計算水平抛射資料〉會介紹另一個寫法，使用 with 開啟檔案，使用這個寫法就不需要關閉檔案。

    ​​​​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。

    ​​​​    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 迴圈運作的條件為

    ​​​​abs(v2 - v1) > eps

    也就是速度變化仍大於我們所設定的精準度時，代表小球尚未達到終端速度，需要繼續運算。接下來在第80行，更新 v1 為 v2 目前的值，更新 v2 為 b1.v.y。

    ​​​​v1, v2 = v2, b1.v.y

  7. 轉存成文字檔的資料格式如下（取得檔案）：













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 繪圖技巧：讀取數據及線性擬合〉、〈Matplotlib 繪圖技巧：在同一張圖上繪製兩個函數、兩張圖並列〉。

2025年3月20日補充：線上版 GlowScript 網站可以執行的版本







































































Web VPython 3.2

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

print("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:
        print(f"{t:f}, {b1.pos.y:f}, {b2.pos.y:f}, {b1.v.y:f}, {b2.v.y:f}, {b1.a.y:f}, {b2.a.y:f}")
        tp = tc
# 更新 v1 為 v2 目前的值，更新 v2 為 b1.v.y
    v1, v2 = v2, b1.v.y
# 更新時間
    t += dt

2022年9月29日補充：使用 NumPy 以及 Matplotlib 繪圖

如果想要使用 Python 讀取文字檔的資料並繪圖，可以使用 NumPy 以及 Matplotlib，程式碼如下。





















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()                                          # 顯示圖片

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
使用 Matplotlib 繪製的 v-t 圖


假設有以下的實驗數據，若推測 y 與 x 成反比，可以使用以下的程式碼作圖並計算最接近直線的斜率及截距。

x 0.2 0.4 0.6 0.8 1.0
y 10.1 4.9 3.3 2.6 1.9 




























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()                              # 顯示圖片

最接近直線的斜率、縱軸截距分別為

slope=2.028088701161563、
intercept=−0.07080253431890154，以下是繪圖成果。

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
使用 Matplotlib 繪製資料點及最接近直線


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
Last changed by 
 
Yi-Zhe Wang
0
8391

Read more

APCS實作題2025年6月第4題：分組遊戲

日期：2025年6月27日 作者：王一哲 ZeroJudge 題目連結 https://zerojudge.tw/ShowProblem?problemid=q839

Jun 27, 2025
APCS實作題2025年6月第3題：貪心闖關

日期：2025年6月26日 作者：王一哲 ZeroJudge 題目連結 https://zerojudge.tw/ShowProblem?problemid=q838

Jun 27, 2025
APCS實作題2025年6月第2題：轉盤得分

日期：2025年6月26日 作者：王一哲 ZeroJudge 題目連結 https://zerojudge.tw/ShowProblem?problemid=q836

Jun 26, 2025
APCS實作題2025年6月第1題：小心陷阱

日期：2025年6月26日 作者：王一哲 ZeroJudge 題目連結 https://zerojudge.tw/ShowProblem?problemid=q836

Jun 26, 2025
Read more from Yi-Zhe Wang
Published on HackMD