# Matplotlib 繪圖技巧：讀取數據及線性擬合
> 作者：王一哲 日期：2018/11/23
<br></br>
## 前言
對於物理實驗而言，我們在處理數據時通常會繪製 **xy 散佈圖**，將應變變因的數值畫在 y 軸、操縱變因的數值畫在 x 軸，如果數據點分布在一條斜直線上，就可以利用**最小平方法**找出最接近直線的斜率和截距，再利用物理定律解釋斜率和截距的成因。但如果數據點不是分布在一條斜直線上，就要先想辦法處理數據，圖形化為直線。
以雙狹縫干涉的實驗為例，假設雷射光波長為$\lambda$、狹縫間距為$d$、狹縫與屏幕間的距離為$L$，則干涉的寬度
$$\Delta y = \frac{\lambda L}{d}$$
如果我們固定$\lambda$和$d$，改變$L$並測量$\Delta y$，再畫出$\Delta y - L$關係圖，這些數據點應該會分布在一條斜直線上，而且斜率 = $\frac{\lambda}{d}$、截距 = 0。
如果我們固定$\lambda$和$L$，改變$d$並測量$\Delta y$，再畫出$\Delta y - d$關係圖，這些數據點會分布在一條曲線上。如果我們猜測$\Delta y$和$d$成反比，應該要計算$1/d$，再畫出$\Delta y - 1 / d$關係圖，如果數據點分布在一條斜直線上，才能證明$\Delta y$和$d$成反比。由理論式可知，直線的斜率 = $\lambda L$、截距 = 0。
<img height="100%" width="100%" src="https://i.imgur.com/cCh0T9Q.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">左：雙狹縫干涉 Δy - d 關係圖；右：雙狹縫干涉 Δy - 1/d 關係圖</div>
<br></br>
## 產生數據資料
我們需要有數據資料才能繪圖，我們利用以下的程式經由理論計算再加上一點隨機誤差，產生$\Delta y$和$d$的數據並儲存成csv檔。如果已經有實驗數據的同學可以直接跳到繪圖的部分。如果想要產生數據並直接繪圖的同學，可以將以下的兩個程式合併。
假設我們使用氦氖雷射進行實驗，則$\lambda = 633 ~\mathrm{nm}$，狹縫與屏幕間的距離設定為 $L = 1 ~\mathrm{m}$，狹縫間距為$d = 0.1, 0.2, ... , 0.5 ~\mathrm{mm}$。為了使數據大一點，以下的長度單位皆為 mm。
<br></br>
```python=
"""
雙狹縫干涉: 產生數據用的程式
日期: 2018/11/22
作者: 王一哲
"""
import numpy as np
# y = \lambda L / d, 由於 lambda 是保留字, 改用 length
# 產生 d = 0.1, 0.2, ... , 0.5 mm, 計算出 y 值 + 隨機誤差
length, L = 633E-6, 1E3
d = np.arange(0.1, 0.6, 0.1)
error = (np.random.random(size = len(d)) - 0.5) * 0.1
y = length * L / d * (1 + error)
data = np.vstack((d, y)).T
np.savetxt("interference_data.csv", data, delimiter = ',', newline = '\n', fmt = '%.2f')
```
<br></br>
1. 由於 lambda 是 Python 的保留字，改用 length 作為波長的變數名稱。
2. 使用 **np.arange(0.1, 0.6, 0.1)** 產生陣列 [0.1, 0.2, ... , 0.5] 並指定給 d。語法為
```python
np.arange[起始值, 結束值, 間距]
```
但是產生的數據不包含結束值。
3. 使用 **np.random.random(size = len(d)) - 0.5** 產生長度與 d 相等、數值介於 -0.5 到 +0.5 之間的浮點數。
4. 使用 **np.vstack((d, y)).T** 將陣列 d 和 y 合併後再轉置，由於這個部分比較不容易理解，詳細的說明在下一段。
5. 使用 **np.savetxt("interference_data.csv", data, delimiter = ',', newline = '\n', fmt = '%.2f')** 將 data 的內容儲存為 interference_data.csv，數據之間以',' 分隔，每一列的資料最末端加上換行符號 \n，數據的格式為 %.2f，取到小數點以下2位的浮點數。以這個實驗而言，能夠測量到 0.01 mm 已經很精準了，所以我們只取到小數點以下2位的數值。如何用電腦分析光的干涉、繞射照片請參考〈[ImageJ 教學：分析光的干涉、繞射照片](https://docs.google.com/document/d/e/2PACX-1vQhw-TqFuEHlZwgtudKii5LyYYDmV0l0o4_GYCFcZQ_mgEgtUE3NmRvDAwD-NYFJomYs8jBI99JdYDG/pub)〉。產生的資料檔如下
```
0.10,6.24
0.20,3.03
0.30,2.15
0.40,1.54
0.50,1.31
```
<br></br>
### numpy.vstack
首先引入函式庫 **import numpy as np**，再用 np.arange 產生陣列 A、B，我們可以看到 np.vstack([A, B]) 的效果是將 A、B 以垂直方向疊起來，C.T 則是將陣列 C 當中的元素行、列對調。
```python
A = np.arange(1, 4, 1)
B = np.arange(4, 7, 1)
C = np.vstack([A, B])
D = C.T
```
這些變數的內容如下
```python
A: array([1, 2, 3])
B: array([4, 5, 6])
C: array([[1, 2, 3],
[4, 5, 6]])
D: array([[1, 4],
[2, 5],
[3, 6]])
```
在上方的程式當中，如果使用 np.vstack((d, y)) 產生的陣列會是
```python
[[d1, d2, d3, d4, d5],
[y1, y2, y3, y4, y5]]
```
但是我們需要的格式是
```python
[[d1, y1],
[d2, y2],
[d3, y3],
[d4, y4],
[d5, y5]]
```
所以需要使用 np.vstack((d, y)).T 將陣列轉置。
<br></br>
## 讀取數據檔並繪圖
我們使用以下的程式讀取將剛才產生的數據檔，將數據處理過後作圖、計算最接近直線的斜率和截距。
```python=
"""
雙狹縫干涉: 讀取數據檔並繪圖的程式
日期: 2018/11/22
作者: 王一哲
"""
import numpy as np
import matplotlib.pyplot as plt
# 從 interference_data.csv 讀取資料, 資料以 ',' 分隔, 將第0、1欄的資料分別存入 d, y，再計算 1 / d 並存入 dx
d, y = np.loadtxt("interference_data.csv", delimiter = ',', usecols = (0, 1), unpack = True)
dx = 1 / d
# 產生計算最接近直線需要用到的陣列 A ，計算直線的斜率 m 和截距 c
A = np.vstack([dx, np.ones(len(dx))]).T
m, c = np.linalg.lstsq(A, y, rcond = -1)[0]
print(m, c)
# 建立繪圖物件 fig, 大小為 12 * 4.5, 內有 1 列 2 欄的小圖, 兩圖共用 x 軸和 y 軸
fig, (ax1, ax2) = plt.subplots(1, 2, sharex = False, sharey = True, figsize = (12, 4.5))
# 設定小圖 ax1 的坐標軸標籤, 格線顏色、種類、寬度, y軸繪圖範圍, 最後用 plot 繪圖
ax1.set_xlabel('d (mm)', fontsize = 14)
ax1.set_ylabel(r'$\Delta y \mathrm{(mm)}$', fontsize = 14)
ax1.grid(color = 'red', linestyle = '--', linewidth = 1)
ax1.set_xlim(0, np.max(d) * 1.1)
ax1.set_ylim(0, np.max(y) * 1.1)
ax1.plot(d, y, color = 'blue', marker = 'o', markersize = 10, linestyle = '')
# 設定小圖 ax2 的坐標軸標籤, 格線顏色、種類、寬度, 最後用 plot 繪圖
ax2.set_xlabel(r'$1 / d \mathrm{(mm)^{-1}}$', fontsize = 14)
ax2.set_ylabel(r'$\Delta y \mathrm{(mm)}$', fontsize = 14)
ax2.grid(color = 'red', linestyle = '--', linewidth = 1)
ax2.set_xlim(0, np.max(dx) * 1.1)
ax2.set_ylim(0, np.max(y) * 1.1)
ax2.plot(dx, y, color = 'blue', marker = 'o', markersize = 10, linestyle = '', label = "Original Data")
ax2.plot(dx, m * dx + c, color = 'orange', linewidth = 2, label = "Fitted Line")
ax2.legend(loc = 'upper left')
# 儲存、顯示圖片
fig.savefig('double_slit_y_d.svg')
fig.savefig('double_slit_y_d.png')
fig.show()
```
<br></br>
基本的繪圖技巧及設定請參考〈[Matplotlib 繪圖技巧：在同一張圖上繪製兩個函數、兩張圖並列](https://hackmd.io/s/Bkd0EMyCm)〉，以下只說明這個程式當中首次用到的功能。
1. 使用 **d, y = np.loadtxt("interference_data.csv", delimiter = ',', usecols = (0, 1), unpack = True)** 從 interference_data.csv 讀取資料，資料以 ',' 分隔，將第0、1欄的資料分別存入 d、y。**np.loadtxt** 的語法如下
```python
np.loadtxt('檔名', delimiter = '分隔符號', usecols = (欄位編號), unpack = [True 或 False])
```
其中欄位編號是從0開始，unpack = True 會將讀取的資料分開分別存入不同的變數中。
2. 假設最接近直線的方程式為 y = mx + c，若改用陣列型式可以寫成 y = Ap，其中 A = [x 1]、p = [m c]。所以我們需要先產生陣列 A
```python
A = np.vstack([dx, np.ones(len(dx))]).T
```
A 的內容為
```
[[ 10. 1. ]
[ 5. 1. ]
[ 3.33333333 1. ]
[ 2.5 1. ]
[ 2. 1. ]]
```
再用 **m, c = np.linalg.lstsq(A, y, rcond = -1)[0]** 計算斜率 m、截距 c。
3. 使用 **plt.plot** 指令繪圖，如果只要畫數據點而不連線，點的大小為10，參數設定為 **marker = 'o', markersize = 10, linestyle = ''**；如果要畫最接近直線，縱軸的資料設定為 **m * dx + c**。
以下圖的數據為例，m = 0.618223336853、c = 0.030780095037，m 的理論值為 $\lambda L = 633 \times 10^{-6} \times 1000 = 0.633$，誤差約為 2.37%，這是因為我們在產生數據時有加上一點誤差值。
<img height="100%" width="100%" src="https://i.imgur.com/cCh0T9Q.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">左：雙狹縫干涉 Δy - d 關係圖；右：雙狹縫干涉 Δy - 1/d 關係圖</div>
<br></br>
## 結語
這是一些使用 numpy 及 Matplotlib 的筆記，我覺得這樣的功能對於處理高中物理實驗數據已經夠用了，如果想要更深入研究 numpy 及 Matplotlib 可以參考官方網站的說明文件。
<br></br>
## 參考資料
1. **numpy.arange**: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.arange.html
2. **numpy.random.random**: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.random.html
3. **numpy.vstack**: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.vstack.html
4. **numpy.ndarray.T**: https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.ndarray.T.html
5. **numpy.savetxt**: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.savetxt.html
6. **numpy.loadtxt**: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.loadtxt.html
7. **使用numpy跟sympy實作Linear regression**: http://terrence.logdown.com/posts/314392-simple-linear-regressionnumpy
8. **numpy.linalg.lstsq**: https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.linalg.lstsq.html
9. **matplotlib.pyplot.plot**: https://matplotlib.org/api/_as_gen/matplotlib.pyplot.plot.html
---
###### tags:`Physics`、`Python`
Sign in with Wallet
Connect another wallet