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

(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àng)有意義的工作.可惜，公式中出現(xiàn)了一些明顯的差錯(cuò)[1]，有必要更正.MATLAB是一個(gè)強(qiáng)大的數(shù)學(xué)工具，可以簡單地解決許多計(jì)算問題.

1 常用極限的差錯(cuò)

(1)

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

sinx～x，(x→0)

(2)

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

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

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

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

>>help limit

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

>>figure,ezplot(y),grid

>>hold on,stem(0,L)

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

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

2 常微分方法的求法

文獻(xiàn)[1]列舉了各類常微分方程的求法公式，這都是高等數(shù)學(xué)教材的內(nèi)容.其實(shí)，對于有顯函數(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個(gè)參數(shù)是微分方程，第2個(gè)參數(shù)是自變量，兩項(xiàng)都放在單引號中.系統(tǒng)默認(rèn)變量為t，所以第二項(xiàng)′x′是必不可少的.

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

3 歐拉公式的驗(yàn)證

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

eix=cosx+ isinx

(13)

利用MATLAB的指令：

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

>>simplify(y)

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

eiπ+1=0

(14)

輸入MATLAB的函數(shù)：

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

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

>>exp(i*pi)+1

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

4 常用不定積分的計(jì)算

>>int(sin(x))

回車立即可得

-cos(x)

>>int(cot(x))

回車立即可得

log(sin(x))

5 常用定積分的計(jì)算

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

(15)

在命令窗口輸入指令：

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

回車立即可得

pi/32

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

(16)

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

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

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

(17)

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

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

結(jié)果為

pi^2/6

用同樣的方法可以驗(yàn)證公式：

(18)

如果要驗(yàn)證公式：

(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→+∞時(shí)，ζ(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ù)的計(jì)算

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

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

回車可得1，這就驗(yàn)證了公式：

(20)

在命令窗口輸入指令

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

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

(21)

在命令窗口輸入指令：

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

回車可得

f(a)

這就驗(yàn)證了如下公式：

(22)

在命令窗口輸入指令：

>>dirac(-x)

回車可得

dirac(x)

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

δ(-x) =δ(x)

(23)

在命令窗口輸入指令：

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

>>simplify(d)

回車可得0，這就驗(yàn)證了如下公式：

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

(24)

不過，用MATLAB不能驗(yàn)證下式：

(25)

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

(26)

這是因?yàn)棣?x)函數(shù)是一種廣義函數(shù).

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

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)

這就驗(yàn)證了公式：

(28)

為了驗(yàn)證定積分：

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

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

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

(30)

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

(31)

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

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

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

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

為了檢驗(yàn)正弦冪函數(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→+∞時(shí)，In→0.

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

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

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

1)厄密微分方程為

y″-2xy′+2ny=0

(33)

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

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

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

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ù).貝塞爾多項(xiàng)式的值可用函數(shù)besselj和bessely計(jì)算.

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可以計(jì)算連帶勒讓德多項(xiàng)式的值.

4) 拉蓋爾微分方程為

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

(36)

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

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

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

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

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

(37)

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

(38)

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

(39)

再設(shè)

(40)

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

(41)

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

(42)

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

(43)

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

(44)

而Hn是n階厄密多項(xiàng)式.概率密度為

(45)

利用MATLAB編程如下：

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

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

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

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

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

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

12 結(jié)論

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