MATLAB在量子力學學習中的應用

周群益，莫云飛，周麗麗，侯兆陽

(1.廣州理工學院，廣東 廣州 510540；2.長沙學院 電子信息與電氣工程學院，湖南 長沙 410022；3.贛南醫學院 信息工程學院，江西 贛州 341000；4.長安大學 理學院應用物理系，陜西 西安 710064)

量子力學中應用了大量的數學公式，朱涵睿等作者總結了一些常用數學工具.這是一項有意義的工作.可惜，公式中出現了一些明顯的差錯[1]，有必要更正.MATLAB是一個強大的數學工具，可以簡單地解決許多計算問題.

1 常用極限的差錯

(1)

因此，sinx與x是等價無窮小量，表示為

sinx～x，(x→0)

(2)

其他8個函數的極限和等價無窮小量分別為：

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

檢驗極限正確與否的最簡單的方法就是用MATLAB計算.例如，對于比正弦函數，在命令窗口用下列指令即可計算極限：

>>syms x,y=sin(x)/x,L=limit(y,x,0)

其中，“>>”是MATLAB的提示符.指令syms x定義符號變量；y表示符號表達式；limit函數求極限，第1個參數y是表達式，第2個參數x是自變量，第3個參數0表示自變量的極限；L表示函數的極限.這些指令和函數可通過help指令了解其功能.例如，在命令窗口輸入：

>>help limit

回車就能說明函數limit的功能和格式.用如下指令還能畫出曲線：

>>figure,ezplot(y),grid

>>hold on,stem(0,L)

其中，指令figure開創圖形窗口；指令ezplot是簡易繪圖指令，第1個參數是表達式，自變量的默認范圍是[-2π,2π]；指令grid加網格；指令hold on保留已畫的曲線；指令stem畫桿線，第1是參數是橫坐標，第2個參數是縱坐標.如圖1所示，表頭顯示比正弦函數，此函數是偶對稱曲線，當x→0時的極限是1.

圖1 比正弦函數的曲線和x趨于0時的極限

2 常微分方法的求法

文獻[1]列舉了各類常微分方程的求法公式，這都是高等數學教材的內容.其實，對于有顯函數解的常微分方程，可以用dsolve函數直接求解.例如，一階線性微分方程：

(11)

可化為標準形式：

(12)

根據公式可解得y=(ex+C)(x+ 1)2.

利用MATLAB的指令：

>>y=dsolve(′(x+1)*Dy=2*y+exp(x)*(x+1)^3′,′x′)

回車立即可得

y =

exp(x)*(x + 1)^2 + C1*(x + 1)^2

函數中的第1個參數是微分方程，第2個參數是自變量，兩項都放在單引號中.系統默認變量為t，所以第二項′x′是必不可少的.

不論是伯努利方程還是常系數線性微分方程，不論是齊次還是非齊次微分方程，有顯函數解的常微分方程都能用dsolve函數求解.

3 歐拉公式的驗證

MATLAB可對復數進行計算，很容易驗證歐拉公式.例如：

eix=cosx+ isinx

(13)

利用MATLAB的指令：

>>y=cos(x)+i*sin(x)-exp(i*x)

>>simplify(y)

化簡的結果為0，這就驗證了歐拉公式.i是MATLAB默認的虛數單位.對于著名公式：

eiπ+1=0

(14)

輸入MATLAB的函數：

>>exp(i*sym(pi))+1

結果就是0.注意：pi是MATLAB的是默認的變量，表示圓周率，本來是一個近似的數值，sym(pi)將數值轉化為符號，可以得到精確的符號結果.直接用如下函數：

>>exp(i*pi)+1

只能得到很小的復數值.

4 常用不定積分的計算

>>int(sin(x))

回車立即可得

-cos(x)

>>int(cot(x))

回車立即可得

log(sin(x))

5 常用定積分的計算

MATLAB的函數int很容易計算常用的定積分.例如，欲檢驗無窮積分公式：

(15)

在命令窗口輸入指令：

>>int(x^2/(1+x^2)^4,0,Inf)

回車立即可得

pi/32

就驗證了式(15).用同樣的方法可以驗證無窮積分公式：

(16)

用int函數還能驗證三角函數的正交性.

6 常用的級數求和公式的計算

MATLAB的級數求和函數為symsum，如果要驗證公式：

(17)

可以在命令窗口用如下指令：

>>syms n,symsum(1/n^2,n,1,Inf)

結果為

pi^2/6

用同樣的方法可以驗證公式：

(18)

如果要驗證公式：

(19)

可以在命令窗口用如下指令：

>>symsum(1/(2*n)^2,1,Inf)

>>symsum(1/(2*n-1)^2,1,Inf)

>>syms n p positive,symsum(1/n^p,n,1,Inf)

其中，參數positive表示變量n和p是正數.結果是

piecewise(1

這是條件函數或分段函數，意思是當1

>>figure,ezplot(zeta(p),[1,10]),grid

如圖2所示，表頭顯示函數ζ(p)，p不限于偶數，也不限于正整數，ζ(p)隨著p的增加而減小，當p→+∞時，ζ(p)→1.

圖2 ζ(p)函數的曲線和極限

如果在命令窗口用如下指令：

>>syms i,symsum(2*i+1,0,n-1)

結果為

n + n*(n -1)

7 梯度算符和拉普拉斯算符

>>syms x y z,r=sqrt(x^2+y^2+z^2)

>>gradient(r)

結果為

x/(x^2 + y^2 + z^2)^(1/2)

y/(x^2 + y^2 + z^2)^(1/2)

z/(x^2 + y^2 + z^2)^(1/2)

gradient函數還能對數值矩陣求數值梯度.拉普拉斯算符為del2，只能對數值矩陣求值.

8 狄拉克函數的計算

MATLAB的狄拉克函數為dirac，在命令窗口輸入指令：

>>syms x a,int(dirac(x-a),-Inf,Inf)

回車可得1，這就驗證了公式：

(20)

在命令窗口輸入指令

>>int(dirac(x-a),a-eps,a+eps)

回車可得1.指令中，eps是MATLAB的默認變量，表示數值之間的最小間隔，其值為2.2204e-16，即2.204×10-16.這就驗證了如下公式：

(21)

在命令窗口輸入指令：

>>syms f(x),int(f(x)*dirac(x-a),x,-Inf,Inf)

回車可得

f(a)

這就驗證了如下公式：

(22)

在命令窗口輸入指令：

>>dirac(-x)

回車可得

dirac(x)

說明狄拉克函數是偶函數:

δ(-x) =δ(x)

(23)

在命令窗口輸入指令：

>>d=f(x)*dirac(x-a)-f(a)*dirac(x-a)

>>simplify(d)

回車可得0，這就驗證了如下公式：

f(x)δ(x-a)=f(a)δ(x-a)

(24)

不過，用MATLAB不能驗證下式：

(25)

也不能根據極限驗證下式：

(26)

這是因為δ(x)函數是一種廣義函數.

9 伽瑪函數和貝塔函數的計算

1)伽瑪函數定義為

(27)

在命令窗口輸入指令：

>>syms x t positive

>>int(exp(-t)*t^(x-1),t,0,Inf)

回車可得

gamma(x)

輸入指令：

>>figure,ezplot(gamma(x)),grid

可畫出函數曲線.如圖3所示，Γ函數在區間(0,+∞)是先降后升的曲線，通過解析延拓，Γ函數可以擴展到小于零的區域.

圖3 Γ函數的曲線和解析延拓曲線

輸入指令

>>gamma(sym(1/2))

回車可得

pi^(1/2)

>>syms p,g=gamma(p)*gamma(1-p)

>>simplify(g)

回車可得

pi/sin(pi*p)

這就驗證了公式：

(28)

為了驗證定積分：

(29)

在命令窗口輸入指令：

>>syms x n a positive

>>int(x^(2*n)*exp(-a*x^2),0,Inf)

回車可得

gamma(n + 1/2)/(2*a^(n + 1/2))

這個結果比式(29)的右邊的公式更簡單.

2) B函數定義為

(30)

B函數與Γ函數的關系為

(31)

B函數是一個二元函數，用如下指令可以畫出B函數的曲面：

>>syms x y,figure,ezsurf(beta(x,y),[0,2]),alpha(0.5)

其中，ezsurf指令畫二元函數的曲面，第1個參數是二元函數，第2個參數是兩個自變量的范圍；alpha(0.5)指令使曲面半透明.如圖4所示，B函數隨著x或y的增加而減小，當x或y→+∞時，B→0.

圖4 B函數的曲面

為了檢驗正弦冪函數和余弦冪函數的定積分公式：

(32)

在命令窗口輸入指令：

>>syms x n positive,In=int(sin(x)^n,x,0,pi/2)

或者

>>In=int(cos(x)^n,x,0,pi/2)

回車可得

beta(1/2,n/2 + 1/2)/2

這是用B函數表示的結果，n不分奇數和偶數，因而更簡單.用如下指令可以畫出曲線：

>>figure,ezplot(In,[0,10]),grid

如圖5所示，n不限于整數，In隨著n的增加而減小，當n→+∞時，In→0.

圖5 正弦和余弦冪函數的定積分曲線In

10 特殊函數的微分方程的求解

在量子力學的計算中要解一些微分方程并用到一些特殊函數[4-6].

1)厄密微分方程為

y″-2xy′+2ny=0

(33)

用下列函數可以求方程的解：

>>y=dsolve(′D2y-2*x*Dy+2*n*y′,′x′)

結果很長.用函數hermiteH可以計算厄密多項式的值.

2) 貝塞爾微分方程為

x2y″ +xy′ +(x2-n2)y=0

(34)

用下列函數可以求方程的解：

>>y=dsolve(′x^2*D2y+x*Dy+(x^2-n^2)*y′,′x′)

結果為

C1*besselj(n,x) + C2*bessely(n,x)

其中，besselj(n,x)表示第n階第一類貝塞爾函數，bessely(n,x)表示第n階第二類貝塞爾函數.貝塞爾多項式的值可用函數besselj和bessely計算.

3) 連帶勒讓德微分方程為

(35)

用下列函數可以求方程的解：

>>y=dsolve(′(1-x^2)*D2y-2*x*Dy+

(l*(1+l)-m^2/(1-x^2))*y′,′x′)

結果也很長.用函數legendre可以計算連帶勒讓德多項式的值.

4) 拉蓋爾微分方程為

xy″+(1-x)y′+ny=0

(36)

用下列函數可以求方程的解：

>>y=dsolve(′x*D2y+(1-x)*Dy+n*y′,′x′)

結果也很長.用函數laguerreL可以計算拉蓋爾多項式的值.

11 MATLAB在量子力學中的應用舉例

以線性諧振子為例[4]，其薛定諤方程為

(37)

其中，x是位移，ψ是定態波函數，?=h/2π是約化普朗克常數，m是振子的質量，E是能量.U是其勢能：

(38)

其中，ω是振子的圓頻率.設一個由能量決定的無量綱參數：

(39)

再設

(40)

α的倒數具有長度的量綱.還設一個無量綱的位移ξ=αx，則薛定諤方程可簡化為

(41)

要使方程有解，能量必須取分立值：

(42)

其中，n是自然數0,1,2,….線性諧振子的波函數為

(43)

其中，Nn是歸一化常數：

(44)

而Hn是n階厄密多項式.概率密度為

(45)

利用MATLAB編程如下：

n=8;syms x%量子數,定義符號變量

Hn=hermiteH(n,x)%厄密多項式

Nn=sqrt(1/sqrt(sym(pi))/2^n/factorial(n));%歸一化常數

Pn=Nn*exp(-x^2/2)*Hn%波函數

figure,ezplot(Pn)%創建圖形窗口,畫波函數曲線

hold on,ezplot(Pn^2)%保持圖像,畫概率密度曲線

grid on,axis tight%加網格,貼框

當量子數n取8時，如圖6所示，波函數是波動的，有8個零點，其值大約在-0.5到0.6之間；概率密度也是波動的，其值不小于零.當n取不同的正整數時，可得類似的曲線(圖略).增加一些指令，還能將量子化的概率密度與經典粒子的概率密度進行比較.

圖6 一維線性諧振子的波函數和概率密度

勒讓德多項式和拉蓋爾多項式在氫原子角向概率密度和徑向概率密度的計算中具有重要應用.

量子力學的計算比較復雜，我們精選了10個范例[7]，從黑體輻射到氫原子光譜的規律，從德布羅意波長到隧道效應，全部做到了可視化.特別是氫原子的量子力學處理，對于不同的量子數，用MATLAB繪制的概率密度的等值面圖，形象描述其分布規律.本書有200多個MATLAB程序，其中有10多個關于量子力學的程序，可供讀者參考.

12 結論

在量子力學的學習中要用大量數學工具，手工計算是必不可少的，更重要的是充分利用計算機語言這種工具.MATLAB是一種計算和繪圖功能都十分強大的工具，不但能夠檢錯，還可以解決許多復雜的計算問題，被譽為“工科神器”.我們建議在大一就開設MATLAB課程，并且將各門學科與MATLAB編程方法結合起來，提高學生提出問題和解決問題的能力.