MATLAB在量子力學(xué)學(xué)習(xí)中的應(yīng)用

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

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

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

1 常用極限的差錯

(1)

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

sinx～x，(x→0)

(2)

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

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

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

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

>>help limit

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

>>figure,ezplot(y),grid

>>hold on,stem(0,L)

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

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

2 常微分方法的求法

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

(11)

可化為標(biāo)準(zhǔn)形式：

(12)

根據(jù)公式可解得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

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

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

3 歐拉公式的驗證

MATLAB可對復(fù)數(shù)進(jìn)行計算，很容易驗證歐拉公式.例如：

eix=cosx+ isinx

(13)

利用MATLAB的指令：

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

>>simplify(y)

化簡的結(jié)果為0，這就驗證了歐拉公式.i是MATLAB默認(rèn)的虛數(shù)單位.對于著名公式：

eiπ+1=0

(14)

輸入MATLAB的函數(shù)：

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

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

>>exp(i*pi)+1

只能得到很小的復(fù)數(shù)值.

4 常用不定積分的計算

>>int(sin(x))

回車立即可得

-cos(x)

>>int(cot(x))

回車立即可得

log(sin(x))

5 常用定積分的計算

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

(15)

在命令窗口輸入指令：

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

回車立即可得

pi/32

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

(16)

用int函數(shù)還能驗證三角函數(shù)的正交性.

6 常用的級數(shù)求和公式的計算

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

(17)

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

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

結(jié)果為

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)

其中，參數(shù)positive表示變量n和p是正數(shù).結(jié)果是

piecewise(1

這是條件函數(shù)或分段函數(shù)，意思是當(dāng)1

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

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

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

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

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

結(jié)果為

n + n*(n -1)

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

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

>>gradient(r)

結(jié)果為

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函數(shù)還能對數(shù)值矩陣求數(shù)值梯度.拉普拉斯算符為del2，只能對數(shù)值矩陣求值.

8 狄拉克函數(shù)的計算

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

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

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

(20)

在命令窗口輸入指令

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

回車可得1.指令中，eps是MATLAB的默認(rèn)變量，表示數(shù)值之間的最小間隔，其值為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)

說明狄拉克函數(shù)是偶函數(shù):

δ(-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)

也不能根據(jù)極限驗證下式：

(26)

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

9 伽瑪函數(shù)和貝塔函數(shù)的計算

1)伽瑪函數(shù)定義為

(27)

在命令窗口輸入指令：

>>syms x t positive

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

回車可得

gamma(x)

輸入指令：

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

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

圖3 Γ函數(shù)的曲線和解析延拓曲線

輸入指令

>>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))

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

2) B函數(shù)定義為

(30)

B函數(shù)與Γ函數(shù)的關(guān)系為

(31)

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

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

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

圖4 B函數(shù)的曲面

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

(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函數(shù)表示的結(jié)果，n不分奇數(shù)和偶數(shù)，因而更簡單.用如下指令可以畫出曲線：

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

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

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

10 特殊函數(shù)的微分方程的求解

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

1)厄密微分方程為

y″-2xy′+2ny=0

(33)

用下列函數(shù)可以求方程的解：

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

結(jié)果很長.用函數(shù)hermiteH可以計算厄密多項式的值.

2) 貝塞爾微分方程為

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

(34)

用下列函數(shù)可以求方程的解：

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

結(jié)果為

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

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

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

(35)

用下列函數(shù)可以求方程的解：

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

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

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

4) 拉蓋爾微分方程為

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

(36)

用下列函數(shù)可以求方程的解：

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

結(jié)果也很長.用函數(shù)laguerreL可以計算拉蓋爾多項式的值.

11 MATLAB在量子力學(xué)中的應(yīng)用舉例

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

(37)

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

(38)

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

(39)

再設(shè)

(40)

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

(41)

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

(42)

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

(43)

其中，Nn是歸一化常數(shù)：

(44)

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

(45)

利用MATLAB編程如下：

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

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

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

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

figure,ezplot(Pn)%創(chuàng)建圖形窗口,畫波函數(shù)曲線

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

grid on,axis tight%加網(wǎng)格,貼框

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

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

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

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

12 結(jié)論

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