Lecture 7 多物件控制（list） - 波動現象模擬

建國高中特色選修課程 - 物理現象的程式設計與模擬
作者：曾靖夫
日期：2018/7/5

一、list（串列）指令簡介

在前幾堂課中，我們做了非常多的動力學分析，也學會了如何插入物件來模擬動態效果。不過之前插入物件的方式都是「一個一個」放進去，但今天如果我們要一次分析100個物件的動力學，難道需要把插入物件的指令寫100次嗎？這樣程式還沒有任何功能，就已經寫幾百行了，很顯然我們需更聰明一點的寫法。在python的內建指令中，「list」指令正好可以做到這件事。

list為python中可以一次儲存多個「數字」、「字串」和「物件」的資料型態，在程式中的寫法只要用一個中括號把我們想要的東西寫進去並用逗號把每一元素隔開即可，例如：[ a, b, c,… ]。以下我們來做一點簡單的練習：

Example 1 : list指令基本練習










#list中的元素是有順序的，了解list中index的編號規則
my_list = [1, 2, 3, 4, 5, 'aa', 'bb']    #在list中塞入7個元素，每個元素可以是不同資料型態
print 'my_list = ', my_list     #印出my_list這個list的內容
print my_list[0]                #印出 index = 0 後面的元素
print my_list[2:]               #印出 index = 2 後面的所有元素
print my_list[:2]               #印出 index = 2 前面的所有元素
print my_list[2:5]              #印出 index = 2~5 之間的元素
print my_list[-1]               #印出 index = -1 後面的元素（最後一個元素）
print my_list[:-2]              #印出 index = -2 前面的所有元素
print my_list[:]                #印出所有元素

如果要清楚知道如何取出或控制list中的元素，我們要知道list的內容當中的編號規則，亦稱為「索引」，英文稱為index，下圖顯示python list的索引值：

上圖中的粗體的數字0~9為list中的10個「元素」，紅字為索引值（index），在python中內建有兩種索引值，上排是順著編、下排是逆著編，有時候如果要取後面那幾個元素，反序的編號很好用。

二、常用的list控制指令

同學們不需要記憶眾多操作list的各種指令，有需要的時候查一下官網就好！

Python3 Data Structures : 5.1. More on Lists

三、好多好多的單擺，擺錘波？！ - Pendulum wave（蛇擺）

在youtube上常有一些科學趣味影片，這是我們非常感興趣的主題之一。現在我們就用剛剛學到的list指令，一次掛上十幾個擺錘，並觀察這一堆擺錘的群體行為。如果你能想通這裡面的物理學，你就可以在程式中模擬Pendulum wave，中文譯為「蛇擺」，它會在特定時間點產生分組現象，同學們可以先在youtube欣賞蛇擺實驗短片。

圖片來源：https://www.youtube.com/watch?v=yVkdfJ9PkRQ

這堂課我們只教你如何掛很多個擺錘，但是如何做出蛇擺，同學們請努力想一想吧！也許你會用到單擺週期的公式，以下列出擺長和週期的關係：

這裡我們給出小幅度單擺的週期，當單擺擺角約為 $5\ {\rm deg}$ 以內時，擺動週期可近似為 $T=2\pi\sqrt{\dfrac{L}{g}}$ ， $L$ 為單擺擺長， $g$ 為重力加速度。

多擺錘系統的範例程式如下：

Example 2 : 用 list 做多擺錘系統模擬







































from vpython import *
g = 9.8                  #重力加速度9.8m/s^2
theta = 10*pi/180        #初始擺角設定
k = 100000               #彈力常數
m = 1.0                  #擺錘的質量
n = 12                   #單擺個數
d = 0.1                  #每個擺錘之間的間隔為d公尺 
size = d/3.5             #擺錘圓球半徑的大小
L = [0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.05]      #每一根擺桿的長度

def SpringForce(r,L):    #擺錘所受的彈力
    return - k*(mag(r)-L) * r / mag(r)

scene = canvas(width=600, height=800,center = vector(0, -L[int(n/2)]/2-d, d*n/2+0.15), range=0.9)     
ceiling = box(pos=vector(0,0,(n-1)*d/2), length=0.03, height=0.001, width=(n-1)*d*1.01, color=vector(0.7,0.7,0.7))    
timer = label(pos=vector(5*d,2*d,d*n/2), box = False, height = 20, color=color.yellow)     #在物件視窗用label做計時器

ball = []      #產生空白的list，稍後從迴圈將擺錘一個一個放進去
string = []    #產生空白的list，稍後從迴圈將擺桿一根一根放進去
for i in range(0,n,1):   #用range指令產生list，內容為0~11，並一一取出每一個元素
    ball.append(sphere(pos=vector(L[i]*sin(theta), -L[i]*cos(theta), d*i), v=vector(0,0,0), radius=size, color=color.yellow))
    string.append(cylinder(pos=vector(0,0,d*i), color=vector(0.7,0.7,0.7), radius=0.001))
    #用append指令將物件一個一個擺入list中，i為range這個list的元素

dt = 0.001    #時間間隔
t = 0         #初始時間
while True:
    rate(1/dt)
    t = t+dt            #累計時間
    Ts = t % 60         #計算「秒 sec」
    Tm = int( t/60 )    #計算「分 min」
    timer.text = str('Timer : \n') + str('%1.0f'%Tm) + str(' min ') + str('%1.2f '%Ts) + str(' sec')

    a = []    #產生一個空白的list，稍後計算每一個擺錘的加速度，並放入此list中
    for j in range(0,n,1):
        string[j].axis = ball[j].pos - string[j].pos                    #計算每一根擺桿的軸方向
        a.append(vector(0,-g,0)+SpringForce(string[j].axis, L[j])/m)    #計算每一個擺錘的加速度
        ball[j].v += a[j]*dt            #計算每一個擺錘的速度
        ball[j].pos += ball[j].v*dt     #計算每一個擺錘的位置

range指令的應用：range(start, stop[, step])，建立一個整數序列的list。以下是幾個簡單的例子，執行後的結果如下：

range(10)           [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]    #start預設值為 index = 0
range(1, 11)        [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]   #step預設值為1
range(0, 30, 5)     [0, 5, 10, 15, 20, 25]
range(0, 10, 3)     [0, 3, 6, 9]
range(0, -10, -1)   [0, -1, -2, -3, -4, -5, -6, -7, -8, -9]

進階作業 7-1

在Pendulum wave的模擬中，各個單擺的擺長是關鍵！請將單擺公式寫入程式中，一次算出所有單擺的擺長，並以list存放這些數據。請根據你設定的初始條件，清楚說明你的Pendulum wave週期為何？何時會看到清楚的分組現象（分三組或分兩組）？

提示：

Pendulum wave的特色：多個擺錘釋放一段時間之後，一定會再次呈現所有擺錘連成一線的現象。這代表在這段時間內，擺長最長的擺錘擺 $N$ 次，第二長的擺 $N+1$ 次，第三長的擺 $N+2$ 次…以次類推。
假設最長擺錘週期為 $T_{0}$ ，假設釋放擺錘經過 $N\times T_{0}$ 時間後連成一直線，則第二長的擺週期為 $\dfrac{NT_{0}}{N+1}$ ，第三長的擺週期為 $\dfrac{NT_{0}}{N+2}$ …以此類推。
根據單擺週期公式 $T=2\pi\sqrt{\dfrac{L}{g}}$ ， $L$ ，則可算出每個單擺的擺長，而且能推算出該Pendulum wave的週期為 $N\times T_{0}$ ，即相鄰兩次連成一線的時間間隔。
舉例來說，你可以在程式中設定 $T_{0}=2$ 秒， $N=24$ 次，則Pendulum wave的週期為 $2 \times 24$ 秒，單擺的個數你可以隨意設定。
假設共有 $n$ 個擺錘，請利用單擺週期公式、for迴圈和range設定 $n$ 的範圍，計算出所有單擺的擺長，並以list的資料型態儲存各擺長數據。
L = […週期公式… for i in range(n)]
其中第 i 個擺錘的擺長為 $L_{i}$ ，但在程式中請注意range(n) list的索引起始值為0。

四、波動的圖像

現在我們模擬一個球掉到一根直立的彈簧上，並且黏在彈簧上端，隨後球受到重力及彈力的作用而持續作上下震盪，我們姑且稱之為震盪子。假設這樣的系統，如果一次排幾十個，然後每一個震盪子後面都落後前一個固定的時間差，將宛如棒球場常見的觀眾群波浪舞，這是我們印象中最直觀的波動現象。

圖片來源(左)：https://www.facebook.com/TheDhakaTimes/videos/624880914255369/?t=7 圖片來源(右)：https://www.youtube.com/watch?v=r0quzKa0yg8

但請注意這是假的波動，真正的力學波需要考慮介質之間的交互作用，這將在本章最末討論真正的波動如何形成！以下我們練習用list指令排出數十個震盪子，使其彼此間固定落後特定時間而形成的波浪現象，如下圖所示：

範例程式如下：

Example 3 : 波動的圖像 - 多個自由落體結合震盪利用時間差形成波動







































from vpython import *
g = 9.8         #重力加速度
k = 80.0        #彈力常數
L0 = 0.6        #彈簧原長
H = 1.5         #球的初始高度
m = 0.5         #球的質量
d = 0.15        #震盪子的水平間距
size = 0.06     #球的大小
n = 40          #球的個數
T = 0.05        #相鄰兩球落下的時間差

scene = canvas(width=800, height=400, center=vector(d*n/2-d/2,L0,0))
floor = box(length=d*n, height=0.005, width=0.3, color=color.yellow, pos=vector(d*n/2-d/2,0,0))    #畫地板

def SpringForce(r,L):    #彈力
    return -k*(mag(r)-L)*r/mag(r)

dt = 0.001    #時間間隔
t = 0.0       #初始時間

ball = sphere(radius = size, color=color.red, pos=vector(0,H,0), v=vector(0,0,0))     #畫球
spring = helix(radius=0.03, thickness=0.01, pos=vector(0,0,0), axis=vector(0,L0,0))   #畫彈簧

touched = False         #設定球是否接觸到彈簧的判斷，初始設定為沒碰到(False)

while True:

    rate(1/dt/3)        #以1/3的速度慢速模擬
    t = t + dt          #計時器
    a = vector(0,-g,0)  #球碰到彈簧前的加速度

    if ball.pos.y - spring.pos.y - spring.axis.y <= size:      #球的下緣接觸到彈簧的條件
        touched = True  #接觸到了，將原先預設False改為True
    if touched:
        spring.axis = ball.pos-spring.pos - vector(0,size,0)   #球接觸到彈簧之後，將彈簧的尾端鎖在球的下緣

    a = a + SpringForce(spring.axis,L0)/m     #球的加速度
    ball.v = ball.v + a*dt                    #球的速度
    ball.pos = ball.pos + ball.v*dt           #球的位置

課堂作業 7-1

以上為一個球作自由落體運動掉落黏在彈簧上作上下振動的模擬，請將此震盪系統用list指令排成40個，並使兩兩之間有固定的時間落後。

步驟分解流程，同學請參考以下提示：

如何擺入n個球和彈簧：

for i in range(0,n,1):
    ball.append(sphere(radius = …, color = …, pos=(i*d,H,0), v=…))
    spring.append(helix(radius=…, thickness=…, pos=(i*d,0,0), axis=…))

如何使下一個球的運動固定落後一個時間差：
```
if t < j*T: a.append(vector(0,0,0))
else: a.append(vector(0,-g,0))
```
j 的值在每一個迴圈之後會增加 1，第 j 個球會在時間為 j*T 秒時開始落下，在達到此時間點之前保持靜止。
如何一次判斷n個接觸點：
```
touched = False            #判斷一個接觸點
touched = [False] * n      #判斷n個接觸點
```
給同學的範例程式中只有一個球，所以只須判斷一個接觸點，所以不需要用list。如果是n個球，就必須用list擺入n個False作為初始設定，意思是：一開始的時候n個球都還沒接觸到彈簧。

小提醒：以上使用到for迴圈或if-else條件判斷，要注意縮排！

五、力學波的模擬

在牛頓力學建構完成後，幾乎所有的科學家都相信，波動現象一定要在介質中傳播，而波到底怎麼傳遞將取決於介質的特性，所以只需要把介質特性研究清楚，我們就可掌握力學波的所有行為。即使一直到今天，你或許已經知道光可以在真空中傳播，但電磁波也會進入介質中，傳遞的時候一樣要考慮到介質的特性。

但波動問題一直是個相當複雜的問題，因為多質點系統的分析一直以來都不是容易的事，而且在數學計算上即使做了很多的近似與簡化還是相當有難度，但程式可以幫助我們高效率處理多體問題的數值計算，我們只要用直覺擺入合理的系統，專心思考物理的問題就好了！以下我們使用程式直接研究物理學中，相當重要的群體行為—波動。

這裡考慮最簡單的彈性「繩波」，所有的繩子都可以看成一串多質點系統，而相鄰的質點之間是用彈簧串接起來的，這是我們腦海中對於繩子最直覺的拆解圖像。參考程式如下：

Example 4 : 力學波 - 縱波的模擬








































from vpython import *

A = 0.4            #振幅
N = 50             #介質個數
omega = 2*pi/1.0   #振動角頻率
size = 0.1         #介質的大小
m = 0.1            #介質的質量
k = 500.0          #每一小段彈簧的彈力常數
d = 0.4            #介質之間的初始間隔
     
scene = canvas(title='Spring Wave', width=1200, height=300, range = 0.4*50/6, center = vector((N-1)*d/2, 0, 0)) 
ball = [sphere(radius=size, color=color.yellow, pos=vector(i*d, 0, 0), v=vector(0,0,0)) for i in range(N)]
#以list comprehensions方式產生N個球
spring = [helix(radius = size/2.0, thickness = d/15.0, pos=vector(i*d, 0, 0), axis=vector(d,0,0)) for i in range(N-1)]
#以list comprehensions方式產生N-1條彈簧


def SpringForce(r):    #彈力
    return - k*(mag(r)-d)*r/mag(r)
      
t, dt = 0, 0.001       #初始時間，時間差
while True:
    rate(1000)
    t += dt

    ball[0].pos = vector(A*sin(omega*t), 0, 0)    ##強迫第一個球沿x方向振動

    #在球與球之間放入彈簧
    #from index 0 to 49, the elements are 0,1, 2, 3, 4,... , 48, totally 49 springs
    for i in range(N-1):
        spring[i].pos = ball[i].pos
        spring[i].axis = ball[i+1].pos - ball[i].pos

    #計算每一個球受相鄰兩條彈簧的彈力
    #from index 1 to 50, the elements are 1, 2, 3, 4,... , 49
    for i in range(1, N):      
        if i == N-1: ball[-1].v += SpringForce(spring[-1].axis)/m*dt                            #最後一個球
        else: ball[i].v += (-SpringForce(spring[i].axis) + SpringForce(spring[i-1].axis))/m*dt  #非最後且非第一個的球

        ball[i].pos += ball[i].v*dt

課堂作業 7-2

請將第一個球的強迫振動改為沿y方向，並鎖住最後一個球使之不動，觀察繩波的傳遞與反射現象。執行畫面如下圖所示：

進階作業 7-2

請將前面分析的球－彈簧系統擺成沿鉛直方向，並使最上端的球在水平面上做圓周運動形成一強迫振動的繩擺。執行結果如下：

你似乎看到繩子的上端波長較長，下端的波長較短，你有猜到原因嗎？

提示：

將原來沿著x方向擺放的球和彈簧，改成沿y方向。
假設初始狀態為第一個球（繩子懸吊點）在 $\left ( A,0,0 \right )$ ，其他球請沿鉛直方向對齊第一個球。
強迫第一個球（最上面）做半徑為A圓周運動，請以 sin 和 cos 函數表示其在x-z平面上的圓周運動（Lecture 2內容）。
你也可以用柱狀物件插入一個轉盤，並用rotate指令控制它的轉動。
請在每一個球上加上重力場的作用 vector(0,-g,0)，以及空氣阻力的作用 -b * ball[i].v / m，式中 b 為空氣阻力係數。請比較有加空氣阻力和沒加空氣阻力有什麼明顯差異。

Array寫法：未來學到 Lecture 8 你也可以使用 array 重做這道題目，程式運行速度應更為順暢。

六、Future more…

2009年時，史丹弗大學數學暨機械工程學系的Joseph B. Keller教授，對一個常見的現象提出了一個耐人尋味的問題：「為什麼綁馬尾的慢跑者，頭髮會像節拍器一樣左右搖擺，而不是隨著身體的律動做上下的晃動？」於是他把馬尾的搖擺視做物理學及工程學上常見的單擺問題，進而做了一堆數學公式推導。他發現當上下運動頻率之半整數倍接近「馬尾單擺」的自然頻率時，方程式的解就會變得不穩定，於是乎馬尾辮子只能雀躍地左右擺而不是上下晃[1][2]。

圖片來源：https://www.mamamia.com.au/why-do-ponytails-swing-side-to-side/

參考文獻：
[1] 可愛馬尾的物理方程式－2012年搞笑諾貝爾物理獎：http://pansci.asia/archives/70394
[2] 馬尾擺盪的秘密─探討影響馬尾擺盪的因素：http://activity.ntsec.gov.tw/activity/race-1/55/pdf/030105.pdf

作業繳交：

平台連結：…

此平台由南港高中-高慧君老師設計建置