前言
對於物理實驗而言，我們在處理數據時通常會繪製 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。
左：雙狹縫干涉 Δy - d 關係圖；右：雙狹縫干涉 Δy - 1/d 關係圖
產生數據資料
我們需要有數據資料才能繪圖，我們利用以下的程式經由理論計算再加上一點隨機誤差，產生 $\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')
- 由於 lambda 是 Python 的保留字，改用 length 作為波長的變數名稱。
- 使用 np.arange(0.1, 0.6, 0.1) 產生陣列 [0.1, 0.2, ... , 0.5] 並指定給 d。語法為
但是產生的數據不包含結束值。np.arange[起始值, 結束值, 間距] - 使用 np.random.random(size = len(d)) - 0.5 產生長度與 d 相等、數值介於 -0.5 到 +0.5 之間的浮點數。
- 使用 np.vstack((d, y)).T 將陣列 d 和 y 合併後再轉置，由於這個部分比較不容易理解，詳細的說明在下一段。
- 使用 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]]
讀取數據檔並繪圖
我們使用以下的程式讀取將剛才產生的數據檔，將數據處理過後作圖、計算最接近直線的斜率和截距。
"""
雙狹縫干涉: 讀取數據檔並繪圖的程式
日期: 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 繪圖技巧：在同一張圖上繪製兩個函數、兩張圖並列〉，以下只說明這個程式當中首次用到的功能。
- 使用 d, y = np.loadtxt("interference_data.csv", delimiter = ',', usecols = (0, 1), unpack = True) 從 interference_data.csv 讀取資料，資料以 ',' 分隔，將第0、1欄的資料分別存入 d、y。np.loadtxt 的語法如下
其中欄位編號是從0開始，unpack = True 會將讀取的資料分開分別存入不同的變數中。np.loadtxt('檔名', delimiter = '分隔符號', usecols = (欄位編號), unpack = [True 或 False]) - 假設最接近直線的方程式為 y = mx + c，若改用陣列型式可以寫成 y = Ap，其中 A = [x 1]、p = [m c]。所以我們需要先產生陣列 A
A 的內容為A = np.vstack([dx, np.ones(len(dx))]).T
再用 m, c = np.linalg.lstsq(A, y, rcond = -1)[0] 計算斜率 m、截距 c。[[ 10. 1. ] [ 5. 1. ] [ 3.33333333 1. ] [ 2.5 1. ] [ 2. 1. ]] - 使用 plt.plot 指令繪圖，如果只要畫數據點而不連線，點的大小為10，參數設定為 marker = 'o', markersize = 10, linestyle = ''；如果要畫最接近直線，縱軸的資料設定為 m * dx + c。
左：雙狹縫干涉 Δy - d 關係圖；右：雙狹縫干涉 Δy - 1/d 關係圖
結語
這是一些使用 numpy 及 Matplotlib 的筆記，我覺得這樣的功能對於處理高中物理實驗數據已經夠用了，如果想要更深入研究 numpy 及 Matplotlib 可以參考官方網站的說明文件。
參考資料
- numpy.arange: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.arange.html
- numpy.random.random: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.random.random.html
- numpy.vstack: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.vstack.html
- numpy.ndarray.T: https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.ndarray.T.html
- numpy.savetxt: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.savetxt.html
- numpy.loadtxt: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.loadtxt.html
- 使用numpy跟sympy實作Linear regression: http://terrence.logdown.com/posts/314392-simple-linear-regressionnumpy
- numpy.linalg.lstsq: https://docs.scipy.org/doc/numpy-1.15.1/reference/generated/numpy.linalg.lstsq.html
- matplotlib.pyplot.plot: https://matplotlib.org/api/_as_gen/matplotlib.pyplot.plot.html
HackMD 版本連結：https://hackmd.io/s/ByN63TEC7
沒有留言:
張貼留言