紙蜻蜓的受風面積與紙蜻蜓落地時間的關係 #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當作是急遽上升好了，v-t圖怎麼可能跟x-t圖一樣？也就是說，不是v-t圖的程式碼錯了，就是x-t圖和v-t圖的程式碼錯了。為了找到bug，我找到半夜兩點，但還是找不到，只好先睡一覺，明天思緒清楚了，再繼續找。

我又從早上9點找到12點，總算找到了。

為了能讓讀者體驗一下我的艱辛，下面是我有問題的程式碼，來找找看，是哪裡出錯了。

  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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
import matplotlib.pyplot as plt

all_position = []
all_speed = []

all_position_2 = []
all_speed_2 = []

all_position_3 = []
all_speed_3 = []

#F=(1/2)CρSV^2
#x = vot + 1/2at^2

min_time = 0.001
data_count = 3000

#初始化
speed = 0
x = 0
#可調整
mass = 10
affected_area = 0



def drag_force(current_speed, affected_area, drag_coefficient = 0.4, air_density = 1.225):
    #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):
    #print("gravity_force:", mass * gravity)
    return mass * gravity

def final_acceleration(gravity_force, drag_force, mass):
    #print("final_acceleration:", (gravity_force - drag_force) / mass)
    return (gravity_force - drag_force) / mass

def delta_time_move(acceleration, min_time, current_speed):
    #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):
    #print("delta_time_speed:", current_speed + acceleration * min_time)
    return current_speed + acceleration * min_time

for i in range(data_count):
    x += delta_time_move(final_acceleration(gravity_force(mass), drag_force(speed, affected_area), mass), min_time, speed)
    speed += delta_time_speed(final_acceleration(gravity_force(mass), drag_force(speed, affected_area), mass), min_time, speed)
    #print(delta_time_speed(final_acceleration(gravity_force(mass), drag_force(speed, affected_area), mass), min_time, speed))
    #print()
    #print("--------------------------------------------")
    #print(">>> now position <<<", x)
    #print(">>> now speed <<<", speed)
    #print("--------------------------------------------")
    #print()
    all_position.append(x)
    all_speed.append(speed)

#print(all_position)
#print(all_speed)

#初始化
speed2 = 0
x2 = 0
#可調整
mass2 = 200
affected_area2 = 100


for i in range(data_count):
    x2 += delta_time_move(final_acceleration(gravity_force(mass2), drag_force(speed2, affected_area2), mass2), min_time, speed2)
    speed2 += delta_time_speed(final_acceleration(gravity_force(mass2), drag_force(speed2, affected_area2), mass2), min_time, speed2)
    #print()
    #print("--------------------------------------------")
    #print(">>> now position <<<", x2)
    #print(">>> now speed <<<", speed2)
    #print("--------------------------------------------")
    #print()
    all_position_2.append(x2)
    all_speed_2.append(speed2)

#初始化
speed3 = 0
x3 = 0
#可調整
mass3 = 200
affected_area3 = 200


for i in range(data_count):
    x3 += delta_time_move(final_acceleration(gravity_force(mass3), drag_force(speed3, affected_area3), mass3), min_time, speed3)
    speed3 += delta_time_speed(final_acceleration(gravity_force(mass3), drag_force(speed3, affected_area3), mass3), min_time, speed3)
    #print()
    #print("--------------------------------------------")
    #print(">>> now position <<<", x3)
    #print(">>> now speed <<<", speed3)
    #print("--------------------------------------------")
    #print()
    all_position_3.append(x3)
    all_speed_3.append(speed3)

plt.subplot(1, 2, 1)
posLine1 = plt.plot(all_position, label="S = " + str(affected_area))
posLine2 = plt.plot(all_position_2, label="S = " + str(affected_area2))
posLine3 = plt.plot(all_position_3, label="S = " + str(affected_area3))
plt.legend(loc='best')
plt.title("x-t graph", fontsize = 24)
plt.xlabel("Time (1 = " + str(min_time) + "sec)", fontsize = 16)
plt.ylabel("x", fontsize = 16)
plt.tick_params(axis='both', labelsize = 12, color='red')

plt.subplot(1, 2, 2)
plt.title("v-t graph", fontsize = 24)
plt.xlabel("Time (1 = " + str(min_time) + "sec)", fontsize = 16)
plt.ylabel("v", fontsize = 16)
posLine1 = plt.plot(all_speed, label="S = " + str(affected_area))
posLine2 = plt.plot(all_speed_2, label="S = " + str(affected_area2))
posLine3 = plt.plot(all_speed_3, label="S = " + str(affected_area3))
plt.legend(loc='best')
plt.tick_params(axis='both', labelsize = 12, color='red')

plt.show()

找的到嗎？我想大家一定連找都不想找，直接看到這裡了吧。那我只好順讀者的意，公布答案。

答案是52, 76, 96的 '+=' 要改成 '='

說明一下 x=1 的意思是把1這個數字存入x變數裡, 會把原本x存的值給拋棄掉。

至於 x+=1 的意思是把原本存在x變數裡的值加上1，最後再存進x裡。

我們用等加速度第二公式求出來的是位移，也就是x的變化量。所以使用 "+=" 將原來的x加上位移(Δx)就是我們要的末x。

至於等加速度第一公式求出來的速度是末速度，不是速度變化量，可以直接把原來的值拋棄，直接存入算出的末速度，也就是使用 "="。

四、模擬數據:

以下是我寫的程式碼：

  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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
import matplotlib.pyplot as plt

all_position = []
all_speed = []

all_position_2 = []
all_speed_2 = []

all_position_3 = []
all_speed_3 = []

#F=(1/2)CρSV^2
#x = vot + 1/2at^2

min_time = 0.001
data_count = 3000

#初始化
speed = 0
x = 0
#可調整
mass = 10
affected_area = 0



def drag_force(current_speed, affected_area, drag_coefficient = 0.4, air_density = 1.225):
    #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):
    #print("gravity_force:", mass * gravity)
    return mass * gravity

def final_acceleration(gravity_force, drag_force, mass):
    #print("final_acceleration:", (gravity_force - drag_force) / mass)
    return (gravity_force - drag_force) / mass

def delta_time_move(acceleration, min_time, current_speed):
    #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):
    #print("delta_time_speed:", current_speed + acceleration * min_time)
    return current_speed + acceleration * min_time

for i in range(data_count):
    x += delta_time_move(final_acceleration(gravity_force(mass), drag_force(speed, affected_area), mass), min_time, speed)
    speed = delta_time_speed(final_acceleration(gravity_force(mass), drag_force(speed, affected_area), mass), min_time, speed)
    #print(delta_time_speed(final_acceleration(gravity_force(mass), drag_force(speed, affected_area), mass), min_time, speed))
    #print()
    #print("--------------------------------------------")
    #print(">>> now position <<<", x)
    #print(">>> now speed <<<", speed)
    #print("--------------------------------------------")
    #print()
    all_position.append(x)
    all_speed.append(speed)

#print(all_position)
#print(all_speed)

#初始化
speed2 = 0
x2 = 0
#可調整
mass2 = 200
affected_area2 = 100


for i in range(data_count):
    x2 += delta_time_move(final_acceleration(gravity_force(mass2), drag_force(speed2, affected_area2), mass2), min_time, speed2)
    speed2 = delta_time_speed(final_acceleration(gravity_force(mass2), drag_force(speed2, affected_area2), mass2), min_time, speed2)
    #print()
    #print("--------------------------------------------")
    #print(">>> now position <<<", x2)
    #print(">>> now speed <<<", speed2)
    #print("--------------------------------------------")
    #print()
    all_position_2.append(x2)
    all_speed_2.append(speed2)

#初始化
speed3 = 0
x3 = 0
#可調整
mass3 = 200
affected_area3 = 200


for i in range(data_count):
    x3 += delta_time_move(final_acceleration(gravity_force(mass3), drag_force(speed3, affected_area3), mass3), min_time, speed3)
    speed3 = delta_time_speed(final_acceleration(gravity_force(mass3), drag_force(speed3, affected_area3), mass3), min_time, speed3)
    #print()
    #print("--------------------------------------------")
    #print(">>> now position <<<", x3)
    #print(">>> now speed <<<", speed3)
    #print("--------------------------------------------")
    #print()
    all_position_3.append(x3)
    all_speed_3.append(speed3)

plt.subplot(1, 2, 1)
posLine1 = plt.plot(all_position, label="S = " + str(affected_area))
posLine2 = plt.plot(all_position_2, label="S = " + str(affected_area2))
posLine3 = plt.plot(all_position_3, label="S = " + str(affected_area3))
plt.legend(loc='best')
plt.title("x-t graph", fontsize = 24)
plt.xlabel("Time (1 = " + str(min_time) + "sec)", fontsize = 16)
plt.ylabel("x", fontsize = 16)
plt.tick_params(axis='both', labelsize = 12, color='red')

plt.subplot(1, 2, 2)
plt.title("v-t graph", fontsize = 24)
plt.xlabel("Time (1 = " + str(min_time) + "sec)", fontsize = 16)
plt.ylabel("v", fontsize = 16)
posLine1 = plt.plot(all_speed, label="S = " + str(affected_area))
posLine2 = plt.plot(all_speed_2, label="S = " + str(affected_area2))
posLine3 = plt.plot(all_speed_3, label="S = " + str(affected_area3))
plt.legend(loc='best')
plt.tick_params(axis='both', labelsize = 12, color='red')

plt.show()

以下是執行後的結果 (將3秒切成3000份) (空氣阻力係數=0.4) (空氣密度=1.225):






五、數據分析

從上圖可以看到，在一開始時，因沒有速度，導致一開始每個物體的空氣阻力都接近0，和自由落體差不多，但隨著物體漸漸開始有了速度，將受風面積的差距以空氣阻力的方式呈現出來，使受風面積越大者離自由落體的曲線越遠。最後因空氣阻力，導致速度變慢，最終以趨近等速度運動落下。

因此，我們可以說：受風面積越大，同高度落下時間越長，但不成正比，只是正相關。
有受到空氣阻力者最終以趨近等速度運動落下。

六、省思

原本我是希望可以藉由查詢的方式，找到我想要的資訊。但我搜尋了一段時間，都找不到我想要的，因此只好自行撰寫程式，了解紙蜻蜓的受風面積與紙蜻蜓落地時間的關係。做完了此活動，以下是我的省思。

(一)知識的學習是為了學習更多的知識


    前面提到，我一年級下學期時接觸了Python，還有我二年級時有學到順時速度與順時加速度，包含了lim和一點微積分的觀念。因為有這些知識，使我現在可以利用這些知識，完成此研究。學而不思則罔，如果只是了解其表面，不與其他的知識做連結，終究只是一個點，而不是個兼顧而實用的網子。

(二)程式設計是一個累人的事情


    我國小時曾經立志要當一個程式設計師，因為感覺就像駭客一樣，看著其他人完全看不懂的資料，將電腦掌握在手中。但經過了幾年的自學，有些問題即使稍微用頭腦轉一下就可以大略知道大略的解決方法，但實際實作時，沒兩三下就遇到錯誤，接下來就是幾小時的與Google互動的時間和忍受一事無成的感覺。但是，一件事最重要的是過程。結果是給人看的，過程是自己體會和學習的漣漪。

七、證明


為了證明這是我做的，我在製作的過程有用GitHub紀錄:https://github.com/DCtime/air_drag


八、下一步

對同學們來說，圖表雖然可以使數據的關係變得一目了然，但還是太枯燥了。為了解決這個問題，我進一步的利用VPython，不只使其可讀化，更使其可視化，詳情請看紙蜻蜓的受風面積與紙蜻蜓落地時間的關係 #2 [可視化使用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 ( "平年" )
閱讀完整內容

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，...
閱讀完整內容