# Lecture 8 進階多物件控制（array）
>建國高中特色選修課程 - 物理現象的程式設計與模擬
>作者：曾靖夫、賴奕帆
>日期：2018/7/5
---
## 一、NumPy Array（陣列）簡介
在Python預設的程式語言中，有list指令來一次儲存眾多元素，但是並沒有array這種資料型態。「array陣列」是Python的另一個套件NumPy（Numerical Python）中所內含資料型態，不過我們平常用的Visual Python套件裡面就含有NumPy，所以不需要再進行額外的程式安裝。
雖然array同樣是多數據處理的指令，但與list指令相比，它們在電腦的記憶體中儲存的方式有很大的差異。同學應該還記得之前使用list指令時，其中每一個元素可以是不同的資料型態，因此這些資料在記憶體中的儲存位置是很難以預測的。而array規定每一個元素都必須有相同的資料型態，它們在記憶體上的儲存位置會完全排在一起，因此array的存取速度會比list快非常多！下圖可以具體比較這兩種資料型態其儲存方式的差異：
<img style="display: block; margin-left: auto; margin-right: auto" width="45%"
src="https://i.imgur.com/e28jewS.png">
因為其卓越的運算效能，故NumPy Array被非常廣泛的應用在科學計算中。
一般常和NumPy搭配使用的套件還有很多，這裡簡單列出兩個最常用的，同學可以自行安裝玩玩看：
1. Matplotlib：可以畫出精美的2維、3維的數據圖和函數圖，畫圖的效果非常類似於Matlab。但重點是Matlab超貴，而Matplotlib免費！
周莫烦 - [Matplotlib 数据可视化神器 Python](https://www.youtube.com/playlist?list=PLXO45tsB95cKiBRXYqNNCw8AUo6tYen3l)
2. Scipy：Scipy是「Scientific Python」的縮寫，內含有解微分、積分、微分方程、傅立葉轉換…等數學的套裝指令，科學運算功能相當強大！
---
## 二、NumPy Array的性能測試
以下程式我們使用三種方式來產生數列「$0^2, 1^2, 2^2, 3^2,..., 9999^2$」總共100次。在執行結果中，你會發使用array運算速度大約快了約100多倍（每台電腦可能不同，要看性能），如果是更複雜的情況，運算速度一定差更多！
### Example 1 : NumPy Array性能測試
```python=
from timeit import default_timer as timer
import numpy as np
#1 先建立空的list再一一擺入元素
start = timer() #計時開始
for x in range(100): #重複做100次
j = [] #產生一個空的list
for i in range(10000): j.append(i**2) #將平方數一個個擺入list
end = timer() #計時結束
print (end - start) #計算時間差
#2 以list comprehensions方式直接建立平方數list
start = timer() #計時開始
for x in range(100): #重複做100次
j=[i**2 for i in range(10000)] #直接建立平方數list
end = timer() #計時結束
print (end - start) #計算時間差
#3 直接建立Array
start = timer() #計時開始
for x in range(100): #重複做100次
i = np.arange(10000) #建立一個0~9999的Array
j = i**2 #將Array裡面的每個元素平方
end = timer() #計時結束
print (end - start) #計算時間差
```
---
## 三、Array的基本運算規則
Array在NumPy中是一個多維度的陣列（稱為ndarray，n-dimensional array），其數據型態擺出來的樣子很像矩陣。如果是一個1維的array，擺出來就是一串數列，看起來和list很像，但事實上它們是不同的資料型態。以下我們提供幾個常見的1維陣列的基本運算範例：
### Example 2 : NumPy Array基本運算
```python=
import numpy as np
# step 1：1維的等差陣列
a = np.arange(0,4.0,0.5) #類似python內建的range，只是輸出是array的資料型態
print (a) #輸出為：[ 0. 0.5 1. 1.5 2. 2.5 3. 3.5]
print (type(a)) #輸出為：<class 'numpy.ndarray'>
# step 2：修改陣列中的內容
a[0] = 5 #將array a中index為0的元素重新指定數值為5
a[-1] = 100 #將array b中的最後一個元素改成100
print ('array a = ', a) #輸出為：[ 5. 0.5 1. 1.5 2. 2.5 3. 100.]
#step 3：先建立list再將轉為array
b = np.array(range(10)) #range(10)為list的資料型態，利用array( )指令就可轉換為array
print (type(range(10))) #輸出為：<class 'range'>
print (type(b)) #輸出為：<class 'numpy.ndarray'>
print ('array b = ', b) #輸出為：[0 1 2 3 4 5 6 7 8 9]
#step 4：對array中的元素進行運算
c = a[1:-1]**2 #將array中index為1~-1取出來，然後每個元素平方
print (c) #輸出為：[ 0.25 1. 2.25 4. 6.25 9.]
d = a[5:]*0.5 #將array中index為5以後的元素取出來，然後每個元素乘0.5
print (d) #輸出為：[1.25 1.5 50.]
#step 5：兩個array之間的運算
e = a[:-1] + b[-7:] #將array a和array b的index 0~-1與-7~最後元素取出相加
print (e) #各取出a和b七個元素後相加，請注意取出的個數要一樣多
```
這裡我們接著介紹2維的array，或者你暫時把它理解為矩陣也沒什麼問題。以下我們給幾個簡單的範例：
### Example 3 : 2維的Array
```python=
import numpy as np
a = np.array([[1, 2], [3, 4]], dtype = float) #建立一個2*2的array，裡面元素型態為float
print (a) #印出array a
print (type(a)) #印出a的資料型態
print (a.shape) #印出a的形狀，輸出結果為(2,2)
```
從程式執行的結果會發現a是一個 $2 \times 2$ 的矩陣，在數學上的表示法為：$$a=\begin{pmatrix}
1 & 2 \\
3 & 4
\end{pmatrix}$$下指令a.shape可給出a形狀，得到它是一個2列（row）、2行（column）的矩陣。其內容為，第1個 row 為 $\begin{pmatrix}1 & 2\end{pmatrix}$、第2個 row 為 $\begin{pmatrix}3 & 4\end{pmatrix}$；第1個 column 為 $\begin{pmatrix} 1\\3 \end{pmatrix}$、第2個 column 為 $\begin{pmatrix} 2\\4 \end{pmatrix}$。
除此之外，如果要將這個陣列中的元素更改、取出或做任何運算，則必須知道這個 array 的 index。如果是2維的就會有兩組 index 來表示各個元素的位置。如下圖所示：紅色數字為兩組 index，axis 為 array 預設的軸方向
<img style="display: block; margin-left: auto; margin-right: auto" width="45%"
src="https://i.imgur.com/4E7Hwip.png">
### Example 4 : 2維Array的索引值
```python=
import numpy as np
a = np.array([[1, 2], [3, 4]], dtype = float) #建立一個2*2的array，裡面元素型態為float
print ('the element at 1st row and 1st column = ', a[0][0]) #第1列，第1行 = 1.0
print ('the element at 1st row and 2nd column = ', a[0][1]) #第1列，第2行 = 2.0
print ('the element at 2nd row and 1st column = ', a[1][0]) #第2列，第1行 = 3.0
print ('the element at 2nd row and 2nd column = ', a[1][1]) #第2列，第2行 = 4.0
```
另外在科學上也常常利用到求和的指令 sum()，就是將矩陣中的元素加總，這常用在求平均值、向量量值、質心或重心位置…等。以下我們作簡單練習：
### Example 5 : 2維Array內的元素求和
```python=
import numpy as np
a = np.array([[1, 2], [3, 4]], dtype = float) #建立一個2*2的array，裡面元素型態為float
b = np.sum(a, axis = 0) #將矩陣的元素沿著直的的方向加起來，並以1維array儲存結果
c = np.sum(a, axis = 1) #將矩陣的元素沿著橫的的方向加起來，並以1維array儲存結果
print (b, b.shape) #可得array b為[4. 6.]的1維array，b.shape = (2,)
print (c, c.shape) #可得array c為[3. 7.]的1維array，c.shape = (2,)
```
:::info
* 1維和2維用shape指令取出陣列形狀的差異：
同學們會發現在這段段程式中，如果有一個1維的array裡面含有N個元素，則該陣列的shape為(N,)，這裡的N指的是元素的個數。雖然它看起來和1列N行的矩陣很像，但是它們實際上差了一個維度，1列N行的矩陣是2維的！舉例來說：
1維5個元素的array為：a = np.array([1, 2, 3, 4, 5]) → a.shape = (5,)
2維1列N行的矩陣為：b = np.array([[1, 2, 3, 4, 5]]) → b.shape = (1,5)
同學該可以明顯看出其差異，1維的array只要用1層中括號，而2維的array需要用到兩層中括號。
:::
這裡還有一些常見NumPy內建可以快速製造array的指令們：
### Example 6 : Array其他常用內建指令
```python=
import numpy as np
print (np.zeros((5,3))) #建立shape為(5,3)，每個元素都是0的Array
print (np.ones((5,3))) #建立shape為(5,3)，每個元素都是1的Array
print (np.diag([1,2,3])) #建立主對角線元素依序為1,2,3的方陣，shape為(3,3)
print (np.random.rand(5,5)) #建立shape為(5,5)的Array，每個元素介於0~1之間均勻隨機產生
print (np.random.randn(5,5)) #建立shape為(5,5)的Array，所有元素的產生呈常態分佈，平均值為0、標準差為1
```
:::info
* 更多 NumPy Array 的內建指令請參考 NumPy 官方網站：
1. [Array creation routines](https://docs.scipy.org/doc/numpy/reference/routines.array-creation.html) : 初階的陣列產生與控制指令。
2. [Random sampling](https://docs.scipy.org/doc/numpy/reference/routines.random.html) : 以陣列產生隨機亂數的各種指令，包含眾多分佈函數。
3. [Linear algebra](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html) : 以陣列進行線性代數等矩陣相關的運算指令。
:::
最後我們介紹更進階的用法－array 的維度擴充，例如：將1維 array 擴充成2維、2維擴充為3維。這裡只用到一個簡單的指令叫做「newaxis」，就可以輕易達成這件事。以下我們看看範例程式：
### Example 7 : Array的維度擴充 - 1維到2維
```python=
import numpy as np
a = np.array([1, 2, 3, 4]) #建立1維的array
b = a[:, np.newaxis] #意義等同array([[1], [2], [3], [4]])
c = a[np.newaxis, :] #意義等同array([[1, 2, 3, 4]])
print ('array b = ', b, 'shape of b is', b.shape) #印出b和其形狀
print ('array c = ', c, 'shape of c is', c.shape) #印出c和其形狀
```
以上程式執行之後你會發現，增加一個維度事實上就是多了一層中括號，只是有兩種括的方式。[:, newaxis]是將第一層中括號內的++每個元素括起來++；而[newaxis, :]是將第一層中括號內的++所有元素整個括起來++。這就使得原來1維4個元素的array分別變成shape是(4,1)和(1,4)的2維array，示意圖如下：
<img style="display: block; margin-left: auto; margin-right: auto" width="55%"
src="https://i.imgur.com/MfnPmfX.png">
根據前面的原理，我們試試將array的維度由2維擴充到3維，請執行以下程式碼：
### Example 8 : Array的維度擴充 - 2維到3維
```python=
import numpy as np
d = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]) #建立shape為(5,3)的array
e = d[:, np.newaxis] #意義同array([[[1, 2, 3]], [[4, 5, 6]], [[7, 8, 9]], [[10, 11, 12]]])
f = d[np.newaxis, :] #意義同array([[[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]]])
print (d, d.shape) #印出d和其形狀
print (e, e.shape) #印出d[:, newaxis]和其形狀
print (f, f.shape) #印出d[newaxis, :]和其形狀
print (e[1][0][2]) #印出d[:, newaxis]的第2列、第1行、第3排之元素 = 6
print (f[0][1][1]) #印出d[newaxis, :]的第1列、第2行、第2排之元素 = 5
```
以下我們用圖示來表示3維array的樣子，同學對照立體圖形看會比較具體：
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/ktUiVsw.png">
這裡的規則我們從程式碼拆解來看，請切記中括號要一層一層拆下去，每個逗點隔開的就是一個元素，看完一層後再繼續往內。以下舉例說明：
<img style="display: block; margin-left: auto; margin-right: auto" width="76%"
src="https://i.imgur.com/tM111vw.png">
<font color = "red">紅色中括號[...]，稱為第一層，這裡面放的元素會往axis = 0的方向擺，因此成為第1列、第2列…</font>。<font color = "blue">藍色中括號[...]，稱為第二層，這裡面放的元素會往axis = 1的方向擺，故成為第1行、第2行…</font>。<font color = "green">綠色中括號[...]，稱為第三層，這裡面放的元素會往axis = 2的方向擺，故成為第1排、第2排…</font>。所以你可以看出來，d[:, newaxis]會是一個共4列、1行、3排的array。如果是以下情況，原理也相同：
<img style="display: block; margin-left: auto; margin-right: auto" width="70%"
src="https://i.imgur.com/s7B8jkK.png">
可看出d[newaxis, :]將會是一個共1列、4行、3排的array。同學再對照上一頁的立體圖形，概念上將會清楚很多！
如果繼續往更高維度探討，可能很容易超出人腦的想像空間，因為我們活在3維空間。但是電腦不會，因為他只要一層一層作中括號下去就好了！
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/j9yzDz6.jpg">
<font size=2> <center>圖片來源：https://www.trin.cam.ac.uk/library/</center></font>
不過我們還是可以再想一下，例如你在圖書館裡也會看到類似的結構，如上圖。接著下個維度可能是「第X區」，再來是「第X樓」，再來是「第X棟」，只是我們不可能一直命名下去，同學腦海中記得這種一層一層的觀念即可！
---
## 四、重繩的靜力平衡問題
模擬一條有質量的重繩，作法是將繩子切成許多帶有質量的小球，然後小球與小球之間用高彈性係數的彈簧連接，而小球的位置、彈簧的位置與軸方向這裡將用 array 來計算。在靜力平衡中有一個著名的懸鍊線問題，就是將繩子的兩端吊著，中間部份的繩子自然垂下，範例程式如下：
### Example 9 : 懸練線問題
```python=
from vpython import*
import numpy as np
N = 50 #質點個數
g = 9.8 #重力加速度
size, m, k, d = 0.016, 0.1/N, N*1000.0, 2.0/N #球大小、質點質量、各力常數、間距
t, dt = 0, 0.0001 #初始時間、時間間隔
scene = canvas(width=1200, height = 600, center = vector(d*N/2.0*0.9, 0, 0))
balls = [sphere(radius=size, color=color.yellow) for i in range(N)]
springs = [helix(radius = size/2.0, thickness = size/5.0) for i in range(N-1)]
ball_pos, ball_v, ball_g = np.zeros((N, 3)), np.zeros((N,3)), np.zeros((N,3))
#以array建立初位置、初速度、重力加速度，全為0的N列3行的陣列
for i in range(N):
ball_pos[i][0] = d*i*0.9 #將球沿x方向的初位置，沿axis=0擺入pos array，並使兩端點距離為0.9倍繩子全長
ball_g[i][1] = -g #將重力加速度陣列的第2行改為-9.8，表示每個球在y方向受到的重力場
while True:
rate(1/dt)
t += dt #計時器
spring_axis = ball_pos[1:] - ball_pos[:-1] #每個彈簧的軸方向
b = np.sum(spring_axis**2, axis = 1)**0.5 #每個彈簧的長度
spring_axis_unit = spring_axis / b[:, np.newaxis] #每個彈簧軸方向的單位向量
fs = - k * (spring_axis - d*spring_axis_unit) #每個彈簧的作用力
fs[b<=d] = 0 #彈簧長度小於原長設為零，表示繩子的鬆弛狀態
ball_v[1:-1] += (fs[:-1] - fs[1:])/m*dt + ball_g[1:-1]*dt - 5.0*ball_v[1:-1]*dt
#計算第二個~倒數第二個球的速度
ball_pos += ball_v *dt #計算球的位置
for i in range(N): #畫球
balls[i].pos = vector(ball_pos[i][0], ball_pos[i][1], ball_pos[i][2])
for i in range(N-1): #畫彈簧
springs[i].pos = balls[i].pos
springs[i].axis = balls[i+1].pos - balls[i].pos
```
[在GlowScript IDE中執行 請點我](https://www.glowscript.org/#/user/cftseng/folder/108-Examples/program/L8-Example-9)
<img style="display: block; margin-left: auto; margin-right: auto" width="85%"
src="https://i.imgur.com/GrP8vRk.png">
以下我們詳細說明 Example 9 懸練線程式的細節：
首先是，++行數#12++，先了建立描述初始狀態的array：
```python=12
ball_pos, ball_v, ball_g = zeros((N, 3)), zeros((N,3)), zeros((N,3))
```
球的初位置、初速度和重力加速度都用先用一個N列3行元素都是0的陣列來建立存取空間，示意圖如下：
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/y08wygT.png">
你會發現只要用一個array就可以將所有球的三個分量一次描述清楚。
++行數#15~17++，將球的初位置沿著 $x$ 方向擺，並將每顆球受到的重力場設定在 $-y$ 方向：
```python=15
for i in range(N):
ball_pos[i][0] = i * d * 0.9
ball_g[i][1] = -g
```
這裡的，數字0.9是將繩子最兩端的球距離設定成0.9倍的繩子原長，這樣在程式執行後比較看得出來繩子受重力下垂的樣子。
++行數#22++，這裡一口氣計算了每個彈簧的軸方向：
```python=22
spring_axis = ball_pos[1:] - ball_pos[:-1]
```
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/5BEtDdC.png">
這裡分兩部份取出前49與後49個球：
1. ball_pos[1:]：取出index = 1到最後一個球的位置，如上圖，即第2～50個質點，共49個質點。
2. ball_pos[:-1]：取出index = 0到最後一個的前一個球的位置，即第1～49個質點，共49個質點。
以下我們用圖形擺起來看看這兩個指令的實際功能：
<img style="display: block; margin-left: auto; margin-right: auto" width="58%"
src="https://i.imgur.com/NzknCjH.png">
ball_pos[1:]和ball_pos[:-1]都是一個(N-1)×3的array，相減後結果如下：
<img style="display: block; margin-left: auto; margin-right: auto" width="72%"
src="https://i.imgur.com/wEmm3LH.png">
看到這裡你會發現，這是array指令發揮的超高效率，我們只透過極為簡單的陣列減法，就一口氣算出所有彈簧的軸方向，當然幾乎等同於得到所有彈簧的彈力了。
++行數#23++，算出每個彈簧的長度（純量）：
```python=23
b = np.sum(spring_axis**2, axis = 1)**0.5
```
這裡的spring_axis**2是將spring_axis這個array裡的所有元素平方
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/G66MQxY.png">
故sum(spring_axis**2, axis = 1)**0.5之輸出結果為：b = array([L1, L2, L3, …, LN-1])，再次強調這是1維N-1個元素的陣列！
++行數#24++：求每個彈簧軸方向的單位向量：
```python=24
spring_axis_unit = spring_axis / b[:, np.newaxis]
```
這裡的b[:, np.newaxis] = array([[L1], [L2], [L3], …, [LN-1]])，變成一個(N-1)×1的array，才能和spring_axis這個(N-1)×3的array相除。運算細節如下
<img style="display: block; margin-left: auto; margin-right: auto" width="72%"
src="https://i.imgur.com/paSWBg9.png">
++行數#25++：計算每一個彈簧的作用力：
```python=25
fs = - k*(spring_axis - d*spring_axis_unit)
```
fs 陣列的每一列都是彈力係數 k 乘上彈簧的伸長向量，故可以同時給出每個彈力的大小和方向，該陣列的每i列代表的就是第i個彈簧作用在第 i+1 個球的彈力向量，負號表示第 i+1 個球受力方向和伸長向量方向相反。
++行數#26++：當彈簧的伸長量小於原長時設定彈力數值為零：
```python=26
fs[b<=d] = 0
```
這是比較沒彈性的繩子模擬（拉緊不易伸長，鬆弛就軟掉），故繩子拉緊的狀態彈力係數k很大，因此伸長量並不明顯。但是繩子鬆弛的時候，應該會軟趴趴的，故此時的繩張力幾乎為零！
為了更了解這個指令的應用，同學可執行下列程式練習：
```python
import numpy as np
a = np.array([3, 5, 2, 4, 1])
print ('a = ', a, ', the elements of a>2 = ', a>2)
b = np.array([[1, 2, 3],[4, 5, 6],[7, 8, 9],[10, 11, 12],[13, 14, 15]])
b[a>2] = 0 #若a>2中的i個元素為True，則以0取代b陣列的第i列
print ('b = ', b)
```
++行數#28++：每個球所受彈力產生的加速度：
```python=28
ball_v[1:-1] += (fs[:-1] - fs[1:])/m*dt
```
這裡只需要處理第2個到第49個球，因為繩子兩端固定不動，端點速度在初始條件就已經設為零。每個球受的彈力分析請看下圖：
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/6I51sUA.png">
每一條彈簧會對兩端的質點施力，其大小相等方向向反。前面的 fs 陣列計算出來的彈力指的是上圖中的紅色作用力，所以 -fs 就是彈簧作用在另一端質點的作用力，如上圖中的綠色作用力。
### 檢核作業 8-1
在剛才的懸練線程式中雖然我們用了 array 來運算，但你會發現運算速度並沒有顯著提昇，主要是因為畫圖的物件數還是太多。事實上我們不需要每一個迴圈都更新畫面，例如我們可以每執行50個迴圈畫一張圖，就可以讓圖形呈現十分順暢，而對於人眼來說時間的精細度也是夠的。最後請用 arrow 繪製懸練線的受力分析，對照正課學過的靜力平衡觀念，分析此模擬的正確性。
<img style="display: block; margin-left: auto; margin-right: auto" width="80%"
src="https://i.imgur.com/FJDgpgC.png">
:::info
* 提示：
1. 設定畫面更新頻率：如 Lecture 1 和 Lecture 2 說明過的畫圖技巧，你可以在迴圈中做如下設定：
```python
...
while True:
...
T = t % (50*dt) #用餘數除法定出要更新圖形的時間點，每50個迴圈會歸零一次
if T + dt >50*dt and T < 50*dt: #在T被歸零前的那瞬間更新畫面圖形
畫圖的指令...
...
...
```
2. 為測試array的性能，請將質點數目增加到100個。
3. 改變球的擺放位置：
```python
for i in range(N):
ball_pos[i][0] = ... #每個球的初始擺放x座標
ball_pos[i][1] = ... #每個球的初始擺放y座標
```
為了使兩端高度不同，這些質點一開始排起來是擺斜的！
4. 如果你覺得一堆球和彈簧會讓畫面看起來比較雜亂，可以刪掉全部畫球和畫彈簧的指令，以曲線物件「curve」來取代。在程式中你可以做如下設定：
```python
...
c = curve([vector(..., ..., ...) for i in range(N)], radius = r)
#畫曲線描點，並給定初始狀態每個點的初位置
while True:
...
for i in range(N):
c.modify(i, pos=vector(ball_pos[i][0],ball_pos[i][1],ball_pos[i][2]))
#在迴圈中用modify更改每個質點的位置，事實上就是將ball_pos的資料抓進來即可
```
5. 兩端點的繩張力：產生兩個 arrow 將第一根與最後一根彈簧的彈力指定為 arrow.axis ，並將 arrow.pos 設定在兩個端點的球上，即可描繪出懸練線端點的受力情形。
6. 繩子所受重力：利用重心公式，以 sum 指令計算出繩子的重心位置，此為重力的作用點，而重力的量值為所有質點的重量和。最後產生 arrow 將 arrow.pos 指定給重心，arrow.axis 指定給所有質點的重力和。
7. 畫輔助線：請將兩個繩張力與重力箭頭的軸方向延伸，以 cylinder 畫出細線，並觀察這這三個箭頭的延伸是否交於一點。這即是課堂上學過的「不平行三力平衡必共點」！
:::
### 檢核作業 8-2
將繩子沿水平方向拉直，然後將其中一端點釋放，並讓繩子上每個質點都受到空氣阻力，請觀察繩子的擺動情形。
<img style="display: block; margin-left: auto; margin-right: auto" width="70%"
src="https://i.imgur.com/7c6shgR.png">
:::info
* 提示：
這裡只要使最後一個質點不固定，讓它受到受到最後一個彈簧的彈力和該質點重力作用即可盪下！
:::
### 指定作業 8
同學可以欣賞以下影片，請從懸練線的程式出發，改寫成軟質的重彈簧系統。將彈簧在均勻重力場中吊起，並使它受到空氣阻力作用而慢慢呈現靜止狀態，然後突然將上端釋放，你將觀察到彈簧的下端幾乎保持靜止，請模擬此現象。
<img style="display: block; margin-left: auto; margin-right: auto" width="90%"
src="https://i.imgur.com/sIE25hY.png">
<font size=2> <center>圖片來源：[Awesome HD Slinky Slow-Mo, Veritasium] https://www.youtube.com/watch?v=uiyMuHuCFo4 </center></font>
:::info
提示：請由範例程式 Example 9 改寫
1. 留下畫彈簧的指令，設每個彈簧的圈數「coils」為1圈，然後刪掉所有畫球的指令，電腦仍有計算球的運動，只是不畫出球。
2. 請在最上面的參數設定，設定彈簧的線條的粗細 springs.thickness = thick = size/5，然後將球與球之間的初始間距（即每個小彈簧原長）設定成 d = thick，彈力係數設定軟一些，例如：k = N * 0.5。
3. 請讓每個質點都受到空氣阻力的作用，否則彈簧系統垂下後會一直振動停不下來。當彈簧保持靜止不動之後，再將最上端的第一個球釋放。
4. 使用條件判斷釋放第一個球，釋放後將空氣阻力關掉。請在迴圈中設定釋放的時間點T0，釋放後代表球1不再固定，將受到它下面彈簧1、重力與空氣阻力的作用，可參考以下寫法：
```python
if t > T0: #設定釋放時機點，t < T0前球1不動，之後球1也受力落下
damp = 0 #掉落後關掉空氣阻力，模擬彈簧在真空環境落下
ball[0] += ... #請寫入球1釋放後的受力，使球1開始運動
```
5. 釋放後請設定彈簧每一圈之間的碰撞為完全非彈性碰撞：
計算每一個圈圈的距離，其實就是球和球之間的距離，請在while迴圈中利用for迴圈算過每個球之間的距離，當球和球間距小於彈簧原長瞬間，使球發生完全非彈性碰撞（碰撞後兩者具有相同速度），寫法如下：
```python
for j in arange(N-1): #由第1根彈簧算到最個最後1根，個數比球少1個
if (ball_pos[j][1] - ball_pos[j+1][1]) <= d: #球和球間距小於彈簧原長
ball_v[j] = ball_v[j+1] = (ball_v[j] + ball_v[j+1])/2 #碰撞後速度
ball_pos[j+1][1] = ball_pos[j][1] - d #碰撞後球的間距
```
這段程式用到了完全非彈性碰撞的公式，假設今天有兩個球質量各為 $m_1$ 和 $m_2$，碰撞後完全黏在一起，稱之為完全非彈性碰撞，碰撞完後的速度為 $\vec{v}_{C}=\dfrac{m_1\vec{v}_{1}+m_2\vec{v}_{2}}{m_1+m_2}$。
6. 運動分析：請畫出整條彈簧的質心與彈簧末端（最後一個球）的v-t圖
質心位置公式：$\vec{r}_{C}=\dfrac{m_1\vec{r}_{1}+m_2\vec{r}_{2}+...+m_N\vec{r}_{N}}{m_1+m_2+...+m_N}$，請用此公式，在迴圈中算出每一時刻的質心位置。因為這裡每個質點質量都相同，程式中可寫為：
```python
ball_cm = np.sum(ball_pos, axis = 0)/N #這是個1維3個元素的array
```
這3個元素分別是質心的x、y、z三個方向的座標位置，質心速度同理。
7. 程式的正確性說明：如果整條彈簧在無空氣阻力的情況下釋放，則彈簧質心的v-t圖應為一條的斜直線，而且斜率的大小恰為重力加速度g。同學們可以參考如下執行畫面：
<img style="display: block; margin-left: auto; margin-right: auto" width="90%"
src="https://i.imgur.com/JQwxqM7.png">
>
8. 執行結果：質心運動為非常完美的等加速度運動（綠色直線），而彈簧的尾端在彈簧完全縮起來之前，幾乎是靜止的（紅色曲線）！另外，物件視窗繪圖請自行調整到適當比例大小，力學要照真實的機制，但繪圖可以不用照真實比例，比如說彈簧的粗細可以粗一點，以免在螢幕上彈簧細到幾乎快看不到。
:::
---
Sign in with Wallet
Connect another wallet