一個拿電腦畫函數圖像的事情，以及一點瞎扯

來自專欄不懂數學的人的物理之路9 人贊了文章

仔細想想自己大學已經過去一半了，然而至今沒做出過什麼像樣的東西，反倒是能不由自主地感受到自己每天都在變菜的事實，之前院里的老師告訴我應當在某一個自己已經選好的道路上堅定地走下去有所建樹，想想也不無道理。開始接觸編程模擬模擬以及MD的各種東西也快一年了，雖然lammps上手沒有多久，還處在基本上不會的階段，但是和一年前比感覺上變菜對時間的一階導似乎變小了不少hhhhh.

說起來，去年這個時候和院里丁老師說自己想進實驗室做項目的最開始，其實是想去量子方向的課題的，包括也十分想去某旦的近代物理實驗室，看了很久的各種版本的量子力學，在這麼多學得都很爛的科目里初等量子力學算是唯一還算可以的科目了，但是不覺間已經了解了這麼久的MD和計算模擬，也稍微有點沒有順從初心的感覺。嘛說到初心這種東西，三年前的現在我應該剛從睿達回來，一邊看著科大的教授在我的筆記本上寫的「麥克斯韋方程組是宇宙的一組史詩」一邊做著考進省隊保送清北最後研究理論物理做黑洞和引力波的春秋大夢呢，誰知道最後一路連崩晃著晃著晃到211去了呢，一晃竟然都過去這麼久了。

前兩天突發奇想感覺自己應該寫點東西，算是對這段時間的工作總結，於是開了這個專欄，有什麼寫什麼吧。

OK廢話完了以後也該開始寫點正經東西了，很多人小時候都有做科學家的夢想，然而大部分都會被現實勸退，當然也有些不知死活的像我這樣的留下來了。以前有個學長說的還不錯，大學數學學得好可以做理論研究，中學數學學得好可以做實驗，像我們這種只有小學數學還可以的（存疑）就只能去玩玩電腦做做計算了。（雖然直到前兩天我發現自己連電腦都玩不好x

就我個人看來計算機的好處其實有兩個，一個在於直接地反映出計算結果，比如你要是想知道本徵解那個一堆加和的是個什麼樣子，那簡單啊，畫個圖看看就行了。還有一個就是算一些不解析的事情。其實說到底解析這件事情，還是因為人類本身的上限問題，雖然很多解析解的得到花了很多技巧費了很大事情，有時候打了幾頁草稿紙才算出來一個表達式，但是其實考慮一些本質上的問題，如果我能把每個點的值都通過迭代計算的方法得出來的話，我又為什麼要知道解析解長什麼樣子呢？而這件人類做不到的事情，（算的夠快的）計算機便可以幫助實現一下。

這裡的問題來自於去年院里老師布置的一道作業，因為是第一稿，就拿它來糊弄一下，大半夜的也不是很想再去構思什麼新的內容，問題很簡單，code也不難，本質上屬於我說的第一類問題，我給個表達式然後看看它是什麼結果。我會在結尾附上源碼，總體來說沒什麼挑戰性，雖然這道題本身我做了很久，但是就這裡的這個方法而言，我其實只花了40分鐘-1個小時就從計算到寫code全部搞定了。

首先是問題：

電子入射到能量為5eV，寬度為1nm的勢壘。繪製隧穿幾率與入射電子能量的曲線。

然後姑且我也算是個學物理的，那就簡單地從薛定諤方程入手解一下這個問題[1]。

該一維空間中的方勢壘可以表達為：

這樣的話我們就可以將各個部分的薛定諤方程寫出來，因為該問題不含時，所以我們直接用定態薛定諤方程套一下：

至於一維定態薛定諤方程嘛，簡單地分離變數化簡一下就可以得到了，這裡就不多bb了。

於是我們有：

先考慮E>U0的情況，這時各個區域的特解可以寫為：

我們從一個稍微物理一些的角度分析一下這個解，畢竟這是不懂數學的人的物理教程。

首先我們基於這樣一個前提：電子的波函數在勢壘的位置是可以發生透射的，這樣才會有3區解存在，然後波矢前的正負號反映了行波的方向，這樣的話顯然，在勢壘的左邊和勢壘的內部都會存在向兩個方向傳播的波，即一個前行波和一個反射波，於是1和2均有兩項作為解的表達式，而3區因為只有透射波，所以只有一項作為解的表達式。

Ok接下來是大家都非常喜歡的小學數學解題過程之解方程組得出係數

由勢壘邊界處的波函數及其微商連續條件，可以得出如下方程：

因為我們需要得到的是隧穿幾率與電子能量的曲線，基於概率幅的定義，其實這就是C方和A方的比值，所以雖然上方程組是不定方程組，我們依舊可以通過它得到我們需要的精確表達式。

簡單的化簡之後我們有：

其中k1和k2是波數，我在這裡還是給出它們的表達式吧：

同樣的我們可以得出E<U0時T的表達式：

其中

剩下來要做的就是把題給數值和一些常數代入進去，然後畫一個還算好看的圖就行了。

拿python大體畫了一下，大概長這麼個樣子：

從中還是可以反映出一些信息的，比如在粒子的能量約為2.5eV的時候會出現隧穿概率突然增大的情況(這件事情我在兩個老師那裡聽到的解釋並不太一樣，雖然我個人認為是共振透射的影響，但是我覺得這其實是一個可以深究的問題)，而在粒子能量略大於5eV時也並非會100%地發生隧穿。這可能就是量子力學所具有的迷人之處之一吧，通過理論計算出的結果並不是完全地符合我們的「想當然」，甚至於可能會顛覆我們的想像。

參考文獻：

[1] 蘇汝鏗. 量子力學[M]. 北京:高等教育出版社, 2002. 52-55

最後附上python的源碼，寫的很醜莫怪，畢竟不專業。

import matplotlib.pyplot as pltfrom matplotlib.pyplot import plot,savefigimport numpy as npimport mathU=5 L=1m=9.1h=6.63/(2*math.pi)i=0step=20000 step=int(step)E=np.zeros((1,step)) k1=np.zeros((1,step)) k2=np.zeros((1,step)) k3=np.zeros((1,step)) A=np.zeros((1,step)) B=np.zeros((1,step))C=np.zeros((1,step))D=np.zeros((1,step))T=np.zeros((1,step))for i in range(0,step-1): E[0,0]=0 k1[0,i]=math.sqrt(2*m*E[0,i])/h if E[0,i]<U: k3[0,i]=math.sqrt(2*m*(U-E[0,i]))/h C[0,i]=(k1[0,i]*k1[0,i]-k3[0,i]*k3[0,i])**2*math.sinh(k3[0,i]*L)**2 D[0,i]=4*k1[0,i]*k1[0,i]*k3[0,i]*k3[0,i] T[0,i]=D[0,i]/(C[0,i]+D[0,i]) else: k2[0,i]=math.sqrt(2*m*(E[0,i]-U))/h A[0,i]=4*k1[0,i]*k1[0,i]*k2[0,i]*k2[0,i] B[0,i]=(k1[0,i]*k1[0,i]-k2[0,i]*k2[0,i])**2*math.sin(k2[0,i]*L)**2 T[0,i]=A[0,i]/(A[0,i]+B[0,i]) E[0,i+1]=E[0,i]+0.001plt.xlim(-1,20)plt.ylim(-0.1,1.1)plt.xlabel("Energy/J")plt.ylabel("Tpossblity")plt.plot(E,T,--*b) plt.show()

推薦閱讀：