Matplotlib 繪圖技巧：讀取數據及線性擬合

作者：王一哲 日期：2018/11/23

前言

對於物理實驗而言，我們在處理數據時通常會繪製 xy 散佈圖，將應變變因的數值畫在 y 軸、操縱變因的數值畫在 x 軸，如果數據點分布在一條斜直線上，就可以利用最小平方法找出最接近直線的斜率和截距，再利用物理定律解釋斜率和截距的成因。但如果數據點不是分布在一條斜直線上，就要先想辦法處理數據，圖形化為直線。

以雙狹縫干涉的實驗為例，假設雷射光波長為

$\lambda$、狹縫間距為

$d$、狹縫與屏幕間的距離為

$L$，則干涉的寬度

$$\mathrm{\Delta }y=\frac{\lambda L}{d}$$
如果我們固定

$\lambda$和

$d$，改變

$L$並測量

$\mathrm{\Delta }y$，再畫出

$\mathrm{\Delta }y-L$關係圖，這些數據點應該會分布在一條斜直線上，而且斜率 =

$\frac{\lambda }{d}$、截距 = 0。

如果我們固定

$\lambda$和

$L$，改變

$d$並測量

$\mathrm{\Delta }y$，再畫出

$\mathrm{\Delta }y-d$關係圖，這些數據點會分布在一條曲線上。如果我們猜測

$\mathrm{\Delta }y$和

$d$成反比，應該要計算

$1/d$，再畫出

$\mathrm{\Delta }y-1/d$關係圖，如果數據點分布在一條斜直線上，才能證明

$\mathrm{\Delta }y$和

$d$成反比。由理論式可知，直線的斜率 =

$\lambda L$、截距 = 0。

左：雙狹縫干涉 Δy - d 關係圖；右：雙狹縫干涉 Δy - 1/d 關係圖

產生數據資料

我們需要有數據資料才能繪圖，我們利用以下的程式經由理論計算再加上一點隨機誤差，產生

$\mathrm{\Delta }y$和

$d$的數據並儲存成csv檔。如果已經有實驗數據的同學可以直接跳到繪圖的部分。如果想要產生數據並直接繪圖的同學，可以將以下的兩個程式合併。

假設我們使用氦氖雷射進行實驗，則

$\lambda =633\mathrm{nm}$，狹縫與屏幕間的距離設定為

$L=1\mathrm{m}$，狹縫間距為

$d=0.1,0.2,...,0.5\mathrm{mm}$。為了使數據大一點，以下的長度單位皆為 mm。

"""
 雙狹縫干涉: 產生數據用的程式
 日期: 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')

1. 由於 lambda 是 Python 的保留字，改用 length 作為波長的變數名稱。
2. 使用 np.arange(0.1, 0.6, 0.1) 產生陣列 [0.1, 0.2, … , 0.5] 並指定給 d。語法為
   ​​​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 教學：分析光的干涉、繞射照片〉。產生的資料檔如下
   ​​​0.10,6.24
   ​​​0.20,3.03
   ​​​0.30,2.15
   ​​​0.40,1.54
   ​​​0.50,1.31
   

numpy.vstack

首先引入函式庫 import numpy as np，再用 np.arange 產生陣列 A、B，我們可以看到 np.vstack([A, B]) 的效果是將 A、B 以垂直方向疊起來，C.T 則是將陣列 C 當中的元素行、列對調。

A = np.arange(1, 4, 1)
B = np.arange(4, 7, 1)
C = np.vstack([A, B])
D = C.T

這些變數的內容如下

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)) 產生的陣列會是

[[d1, d2, d3, d4, d5],
 [y1, y2, y3, y4, y5]]

但是我們需要的格式是

[[d1, y1],
 [d2, y2],
 [d3, y3],
 [d4, y4],
 [d5, y5]]

所以需要使用 np.vstack((d, y)).T 將陣列轉置。

讀取數據檔並繪圖

我們使用以下的程式讀取將剛才產生的數據檔，將數據處理過後作圖、計算最接近直線的斜率和截距。

"""
 雙狹縫干涉: 讀取數據檔並繪圖的程式
 日期: 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()

基本的繪圖技巧及設定請參考〈Matplotlib 繪圖技巧：在同一張圖上繪製兩個函數、兩張圖並列〉，以下只說明這個程式當中首次用到的功能。

1. 使用 d, y = np.loadtxt("interference_data.csv", delimiter = ',', usecols = (0, 1), unpack = True) 從 interference_data.csv 讀取資料，資料以 ',' 分隔，將第0、1欄的資料分別存入 d、y。np.loadtxt 的語法如下
   ​​​np.loadtxt('檔名', delimiter = '分隔符號', usecols = (欄位編號), unpack = [True 或 False])
   

   其中欄位編號是從0開始，unpack = True 會將讀取的資料分開分別存入不同的變數中。
2. 假設最接近直線的方程式為 y = mx + c，若改用陣列型式可以寫成 y = Ap，其中 A = [x 1]、p = [m c]。所以我們需要先產生陣列 A
   ​​​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%，這是因為我們在產生數據時有加上一點誤差值。

左：雙狹縫干涉 Δy - d 關係圖；右：雙狹縫干涉 Δy - 1/d 關係圖

結語

這是一些使用 numpy 及 Matplotlib 的筆記，我覺得這樣的功能對於處理高中物理實驗數據已經夠用了，如果想要更深入研究 numpy 及 Matplotlib 可以參考官方網站的說明文件。

參考資料

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