紙蜻蜓的受風面積與紙蜻蜓落地時間的關係 #2 [可視化使用VPython]

-

作品的網址在這裡:(https://www.glowscript.org/#/user/DCtime/folder/MyPrograms/program/MyFirstVPython) 

一開始，我是使用 Python + Matplotlib 畫出有空氣阻力落體的x-t圖和v-t圖。( Python + Matplotlib 設計流程請看 https://dctimelearninghowtocode.blogspot.com/2021/01/python-matplotlib.html )老師看完我的專案後，提議可以使用VPython將得到的數據用3D模型"直接看到"東西在移動，不只以圖表呈現。因此我就開始製作這個專案。

一、學習VPython

學習國英數，課本、習作、講義、學校授課配合政策很輕易的成為我們學生的學習資源，但學校沒有自動提供VPython的學習資源，因此學習VPython的第一步，就是尋找可使用的學習資源。老師在推薦VPython，有提到可以去Google搜尋 "VPython石明豐教授" 關鍵字。沿著這個線索，找到了一個NTU做的網站(https://tcjd71.wixsite.com/vpython)，裡面有許多教學文章提供下載。再搭配一些Youtube的資源，和先前學的一些Python的知識，不到一天就把我需要的知識全部學完了。

VPython是一個數據可視化的一個工具，簡單來說，就是一個3D動畫製作工具，而非一個物理模擬工具。VPython只提供圖形渲染的工具，移動只能透過直接設定位置的方式，而不是製造環境去使物件動起來。

二、製作物理動畫

正式使用VPython前，要先設定VPython的開發環境。我試著用終端機下載VPython，但她總是跳出一些奇怪的錯誤訊息，可能是因為我的設備太老舊，導致有不相容的問題。為了解決這個問題，我將錯誤碼貼上Google，希望可以找到解決方法。就在此時，我找到了GlowScript，一個可在線上開發的VPython編輯器。雖然沒有像Pycharm那麼好用，但至少解決環境的問題。

(圖一：GlowScript 主頁面, 資料來源：螢幕截圖)

接下來，就是產生出物體有空氣阻力之落體之相關數據(速度與位移)。在上一篇的文章，我們已經完成這個步驟，Ctrl+C + Ctrl+V，這步驟就完成了。( Python + Matplotlib 設計流程請看 https://dctimelearninghowtocode.blogspot.com/2021/01/python-matplotlib.html )

有了數據，接著使用VPython將冷冰冰的數據轉化成3D動畫。在此模擬中，我做兩個長方體，一個受風面積大，另一個受風面積小。這樣就可以比較兩者的掉落情形。

(圖二：一個長方體掉落，資料來源：螢幕截圖)

但是只有兩個物體在掉落，只能看出誰掉得比較快，無法看出掉落的速度等等。因此，要為這兩個長方體畫出 x-t 圖和 v-t 圖，才能精確地看出掉落的過程。

另外，為了確定兩個物體的運動方式是正確的，在此再製作一個圓球，代表自由落體，照理來講，如果有空氣阻力的物體掉落的速度比自由落體快，那就代表此程式有問題，需要修正。

(圖三：完成的專案，資料來源：螢幕截圖)

經過多次的修正，以上就是最後完成的樣子。

從這個動畫可以看到，受風面積越大，最終掉落速度越小，掉落的距離越短。而且有受到空氣阻力的物體，最後會呈現一個趨近等速的狀態，自由落體則是以等加速度落下。


以下為此專案的程式碼：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
GlowScript 3.0 VPython

# ---animation settings---
time_for_a_loop = 1.0 / 30.0;
max_frame = 10000;
time = 0

# ---environment settings---
drag_coefficient = 0.4
air_density = 1.225
gravity = 9.8

# ---box1 status---
length_1 = 2
width_1 = 2
height_1 = 1

mass_1 = 10
speed_1 = 0
position_1 = vector(0, 0, 0)

# ---box2 status---
length_2 = 1
width_2 = 1
height_2 = 1

mass_2 = 10
speed_2 = 0
position_2 = vector(5, 0, 0)

# ---free-fall ball status---
mass_3 = 10
speed_3 = 0
position_3 =vector(10, 0, 0)

# --- small def---
def drag_force(current_speed, affected_area, drag_coefficient = 0.4, air_density = 1.225):
    '''cal the air drag force'''
    #print("drag_force:", drag_coefficient * air_density * affected_area * current_speed * current_speed / 2)
    try:
        return drag_coefficient * air_density * affected_area * current_speed * current_speed / 2
    except ZeroDivisionError:
        return 0
    
def gravity_force(mass, gravity = 9.8):
    '''cal the gravity affected on the item'''
    #print("gravity_force:", mass * gravity)
    return mass * gravity

def final_acceleration(gravity_force, drag_force, mass):
    '''cal the acceleration with g force and drag force and mass effected'''
    #print("final_acceleration:", (gravity_force - drag_force) / mass)
    return (gravity_force - drag_force) / mass

def delta_time_move(acceleration, min_time, current_speed):
    '''cal the delta time movement'''
    #print("delta_time_move:", current_speed * min_time + acceleration * min_time * min_time / 2)
    return current_speed * min_time + acceleration * min_time * min_time / 2

def delta_time_speed(acceleration, min_time, current_speed):
    '''cal the delta time speed (final speed)'''
    #print("delta_time_speed:", current_speed + acceleration * min_time)
    return current_speed + acceleration * min_time

# --- main def ---
def cal_speed(mass, gravity, current_speed, affected_area, drag_coefficient, air_density, min_time):
    return delta_time_speed(final_acceleration(gravity_force(mass, gravity), drag_force(current_speed, affected_area, drag_coefficient, air_density), mass), min_time, current_speed)

def cal_position(now_position, mass, gravity, current_speed, affected_area, drag_coefficient, air_density, min_time):
    change_vector = vector(0, -1 * delta_time_move(final_acceleration(gravity_force(mass, gravity), drag_force(current_speed, affected_area, drag_coefficient, air_density), mass), min_time, current_speed), 0)
    return now_position + change_vector
    
# make boxes
box1 = box(pos=position_1, length=length_1, width=width_1, height=height_1, make_trail=True, trail_type="points", color=color.red)
box2 = box(pos=position_2, length=length_2, width=width_2, height=height_2, make_trail=True, trail_type="points", color=color.blue)
sphere3 = sphere(pos=position_3, color=color.yellow, make_trail=True, trail_type="points", color=color.yellow)

# make speed graph window
speed_graph = graph(xtitle="Time", ytitle="speed")
speed_curve_1 = gcurve(color=color.red, label="Area = " + str(width_1 * length_1))
speed_curve_2 = gcurve(color=color.blue, label="Area = " + str(width_2 * length_2))
speed_curve_3 = gcurve(color=color.yellow, label="Free-fall")

# make position graph window
pos_graph = graph(xtitle="Time", ytitle="position")
pos_curve_1 = gcurve(color=color.red, label="Area = " + str(width_1 * length_1))
pos_curve_2 = gcurve(color=color.blue, label="Area = " + str(width_2 * length_2))
pos_curve_3 = gcurve(color=color.yellow, label="Free-fall")

for i in range(max_frame):
    rate(1.0 / time_for_a_loop)
    
    # cal time
    time = time + time_for_a_loop
    
    # cal speed
    speed_1 = cal_speed(mass_1, gravity, speed_1, length_1 * width_1, drag_coefficient, air_density, time_for_a_loop)
    speed_2 = cal_speed(mass_2, gravity, speed_2, length_2 * width_2, drag_coefficient, air_density, time_for_a_loop)
    speed_3 = cal_speed(mass_3, gravity, speed_3, 0, drag_coefficient, air_density, time_for_a_loop)
    speed_curve_1.plot(time, speed_1)
    speed_curve_2.plot(time, speed_2)
    speed_curve_3.plot(time, speed_3)
    
    # cal pos
    box1.pos = cal_position(box1.pos, mass_1, gravity, speed_1, length_1 * width_1, drag_coefficient, air_density, time_for_a_loop)
    box2.pos = cal_position(box2.pos, mass_2, gravity, speed_2, length_2 * width_2, drag_coefficient, air_density, time_for_a_loop)
    sphere3.pos = cal_position(sphere3.pos, mass_3, gravity, speed_3, 0, drag_coefficient, air_density, time_for_a_loop)
    pos_curve_1.plot(time, box1.pos.y * -1.0)
    pos_curve_2.plot(time, box2.pos.y * -1.0)
    pos_curve_3.plot(time, sphere3.pos.y * -1.0)


三、用途

回到最初做得那受風面積與落地時間之關係的實驗：我們預測受風面積越大，落地時間越久，並從數據得知紙蜻蜓的受風面積與落地時間呈現正相關，但不是正比或一次函數的關係，此數據符合我們之預測。但我們試著去使用公式去解釋為何他會有此行為時，因數學程度不足，所以使用此程式去模擬此公式的現象：用VPython去吸引同學目光，藉由程式模擬，在不使用高深數學的狀況下使同學信服。最後，因為一些原因，並沒有將此專案搬上講臺，有些可惜。

四、製作的過程，我學到...

（一）電腦是計算機

小時候接觸電腦時，我一直認為電腦只是個遊樂器材，一個無限遊戲的儲藏庫。但隨著我電腦知識的累積，知道電腦的歷史，開始了解電腦是一個高級計算機，一個可以用撰寫程式或圖形介面使其計算出我要的答案。玩遊戲也是藉由其快速的計算達成我們的目的（螢幕投射出來，使我們獲得成就感和使我們開心的畫面）。

（二）電腦可應用在科學研究裡

最近所上的探究與實作課程，使我們可以體驗和學習科學研究的方法。電腦可以在提出預測時使用，也可以用來分析數據，彙整成圖表，更可以在實驗完成，建立模型後，利用程式去跑，比對預測是否符合。


















留言

張貼留言

這個網誌中的熱門文章

Zerojudge 基礎題庫a004 文文的求婚 (Python)

-
圖片
西元年被4整除且不被100整除，或被400整除者即為閏年 => year % 4 == 0 and year % 100 != 0 or year % 400 == 0 如果是True就是閏年 False 則是平年 對!就這樣~ Python 程式碼: 1 2 3 4 5 6 7 8 9 10 while 1 : try : x = int ( input ()) except : break if x % 4 == 0 and x % 100 != 0 or x % 400 == 0 : print ( "閏年" ) else : print ( "平年" )
閱讀完整內容

紙蜻蜓的受風面積與紙蜻蜓落地時間的關係 #1 [實驗歷程與Python Matplotlib]

-
圖片
一、問題: 最近在探究與實作中，我和我同學要以紙蜻蜓為實驗研究主題，完成一個實驗報告。為了簡單且方便，我和我的組員研究紙蜻蜓的受風面積與紙蜻蜓落地時間。 甚麼是紙蜻蜓? 紙蜻蜓是一個古時候的玩具，可自行用紙摺出來，當從高處往下掉，紙蜻蜓會旋轉且緩慢地掉落。 最後我們做出的實驗數據如下圖: ▲圖1:實驗數據( 紙蜻蜓的受風面積與紙蜻蜓落地時間之關係圖) 從上圖看，當受風面積越大，落下時間越長。一切看起來合理，但是我們要以公式來論證時，出現了問題。先來看一下公式: 空氣阻力的公式：F=(1/2)CρSV^2 牛頓第二運動定律：F=ma 等加速度第一公式：V = Vo + at 等加速度第二公式： Δx = Vot + (1/2)at^2 假想有一個長方體，一次是直立起來自由落體，另一次是橫躺的自由落體。一開始因為初速都是0，所以利用空氣阻力公式，空氣阻力正比於受風面積，再利用牛頓第二運動定律，得到受風面積與加速度呈負相關(因為阻力的方向與速度方向一定相反，不論阻力多大，速度一定是向下，也就是與阻力方向相反)所以可以推測直立自由落體的長方體會比橫躺的掉下去比較多( Δx)，但隨著時間的推移，直立的因一開始的加速度比較大，速度較大。空氣阻力的公式告訴我們，速度越大，阻力越大，而且還是二次方，所以直立的阻力會比橫躺的阻力比較大。阻截一下：直立的長方體受風面積較小，但速度較大，無法知道誰落下時間較久(阻力影響速度，速度影響阻力)。 二、解決方法: 一年級下學期的自主學習時間，我學Python物件導向程式設計。因為物件導向是Python入門的其中一部份，沒有一本書只寫物件導向，所以我直接買一本Python的入門書。在這本書裡，有提到matplotlib模組，這模組是專門將資料轉成圖表的。 ▲圖2:我的Python入門書 既然腦袋想不出來，那就給電腦做吧。電腦一開始的目的也是用來計算一些軍事的東西，也就是物理的東西。 首先，因為他是變加速度運動，所以會用到微分..應該吧。我們可以把整個運動的x軸切成好多個部分，並把每一個小部分當成等加速度運動去處理。 三、遇到的困難: 我一開始做這個作品時，已經是晚上一點了。照理來講，如果受風面積是0，那應該是自由落體，也就是x-t圖是某種彎曲向上的直線，v-t圖則是通過原點的斜直線，但運行後的結果如下圖: 很明顯的，怎麼可能是鉛直線？即使x-t當...
閱讀完整內容

Zerojudge 基礎題庫a013 羅馬數字 (Python)

-
圖片
一、題目重述 這題真的需要好好講一下，因為題目敘述真的好長。 我重新敘述一下題目好了: 輸入: 兩個羅馬數字 輸出: 一個羅馬數字，值為輸入的兩值相減 二、解法大綱 首先，羅馬數字只是用來單純標示他的值，無法做計算 而且電腦只看得懂阿拉伯數字，所以在做計算時，須先把羅馬轉成阿拉伯 因為題目是輸入兩個羅馬數字，無法直接做相減的動作，所以第一步先把羅馬轉成阿拉伯。 轉成阿拉伯後，就可以相減了，記得加絕對值，因為題目說可能第一數比較小。 計算出答案後，最後把它變成羅馬數字，就可以了。 三、解法問題 有兩個問題，也就是題目的主要部分。 1.羅馬轉成阿拉伯 我們先來觀察一個羅馬數字: III 我想聰明的你，這一定是三，那我問: 你怎麼知道? 我先給一下轉換表: I 1  V5  X10  L50  C100  D500  M1000 好，回來正題，III可以看成是 I + I + I 轉換後也就是 1 + 1 + 1，答案是3 再一個羅馬數字 IV 如果你用上面的方法，1 + 5 = 6，可是答案是四耶~ 那是因為羅馬數字有一個特殊的規則， 數碼限制： 同一數碼最多只能連續出現三次，如40不可表示為XXXX，而要表示為XL (羅馬數字 - 維基百科) 如果我們使用IIII來表示4，感覺I太多了，所以用IV來表示，意思是 5 - 1 換一個方式來表達，也就是 -1 + 4 那  MCDXXXVII  呢?我們先把一個一個分開來，變成: M, C, D, X, X, X, V, I, I 接下來，一一轉成數字，變成: 1000, 100, 500, 10, 10, 10, 5, 1, 1 接著，我們使用一個小伎倆: 在較大的羅馬數字的左邊記上較小的羅馬數字，表示大數字減小數字 我們先找右邊數字比自己大的，找到就把它變成負的 1000, 100, 500, 10, 10, 5, 1, 1 1000右邊是100，比1000小，不變成負的，下一個 100右邊是500，比100大，把它變成負的(100 => -100) 500右邊是10，...
閱讀完整內容