魏健達(dá)，張江敏

(福建師范大學(xué)物理與能源學(xué)院，福建 福州 350117)

現(xiàn)代技術(shù)的進(jìn)步使工程師們需要承接更復(fù)雜且昂貴的項目，這些項目非?？粗叵到y(tǒng)的可靠性和安全性．比如：建筑工程、航天飛行、核能開發(fā)等，需要技術(shù)人員對這些復(fù)雜的物理系統(tǒng)建立相應(yīng)的數(shù)學(xué)模型用于計算．而自然科學(xué)(如熱力學(xué)、流體力學(xué)、固體物理等)總是需要借助偏微分方程[1-3]來對物理系統(tǒng)進(jìn)行描述，但是往往理論上解析求解一些復(fù)雜的偏微分方程是非常困難的，需要借助工具軟件做數(shù)值仿真．仿真求解最常用的一種方法就是有限元方法[4-5]．有限元一詞最早由R.W.Clough提出[6]，現(xiàn)在科學(xué)家們普遍認(rèn)為有限元[7-8]的根源可以追溯到3個獨立的研究小組，即應(yīng)用數(shù)學(xué)家Courant；物理學(xué)家Synge以及工程師Argyris、Kelsey．Taylor的著作The Finite Element Method和麻省理工學(xué)院Klaus-Jürgen Bathe教授的著作Finite Element Procedures中均運用有限元方法解決了許多固體力學(xué)和流體力學(xué)的問題．自 20 世紀(jì) 60 年代以來，有限元方法被廣泛應(yīng)用,比如：1973年工程師謝干權(quán)和李建華攻克的三維有限元程序應(yīng)用于對全國大型水壩的計算；在生物力學(xué)和手術(shù)治療方面李杰，Deng等[9-10]采用有限元方法研究頂點旋轉(zhuǎn)扳法對頸椎引力的影響，達(dá)到減小頸椎間盤應(yīng)力集中的效果．

在量子力學(xué)的一維勢阱問題中，人們不會考慮邊界對于勢阱的影響．但是，當(dāng)研究二維勢阱問題時，就會出現(xiàn)許多不同的邊界(比如：圓形邊界、橢圓形邊界以及本文要用有限元方法仿真的多角形邊界)，這些邊界對求解的問題的影響是巨大的．本文從二維無限深方勢阱出發(fā)，用有限元方法去研究不同邊界勢阱中粒子的本征問題和概率問題，將有限元這種數(shù)學(xué)中解偏微分方程的方法與量子力學(xué)中的問題結(jié)合起來，提供了研究量子力學(xué)問題的新思路，為進(jìn)一步探索量子世界開辟新的道路．

有限元方法是把連續(xù)的結(jié)構(gòu)離散成有限個單元和有限個節(jié)點，再對其組合體進(jìn)行分析[11-12]，其本質(zhì)是將大型物理系統(tǒng)細(xì)分為更小、更簡單的部分．它通過空間維度上特定的空間離散化來實現(xiàn)，將物理系統(tǒng)離散化即可精確表示復(fù)雜的幾何形狀，又可以描述材料特性，從而輕松表示整體解決方案以及精確描述局部的現(xiàn)象．本文利用有限元方法，首先將多角形勢阱網(wǎng)格化成矩形單元；然后借助雅可比矩陣轉(zhuǎn)化坐標(biāo)系，耦合每個矩形單元；最終實現(xiàn)薛定諤方程和概率密度函數(shù)的求解．

1 有限元方法

1.1 二維薛定諤方程的弱形式[13]

在二維無限深勢阱中薛定諤方程和邊界條件表示如下：

(1)

(2)

1.2 物理坐標(biāo)系與自然坐標(biāo)系的轉(zhuǎn)換

使用物理坐標(biāo)系可以直接計算波函數(shù)圖像，但對于數(shù)值計算并不方便．因此，將波函數(shù)分割成物理坐標(biāo)系下規(guī)律不明顯的四邊形小塊，再利用形函數(shù)[17-19]將物理坐標(biāo)系中規(guī)律不明顯的四邊形小塊(如圖1(a))轉(zhuǎn)化成自然坐標(biāo)系中的矩形單元(圖1(b))．

圖1 坐標(biāo)的變換Fig.1 Transformation of coordinates

其物理坐標(biāo)系的x和y與自然坐標(biāo)系下的ξ和η之間的轉(zhuǎn)化公式是：

(3)

而ψ對ξ與η的偏導(dǎo)數(shù)，可以通過J矩陣和ψ對x與y求偏導(dǎo)數(shù)聯(lián)系起來.

(4)

其中J矩陣可以由形函數(shù)推導(dǎo)出來，在本文中稱它是聯(lián)系物理坐標(biāo)系與自然坐標(biāo)系的Jacobian矩陣[20-21].

1.3 數(shù)值積分多項式

(5)

2 方勢阱中的解析解方法與數(shù)值方法

2.1 粒子在方勢阱中的解析解方法

人們對于求解方勢阱問題比較熟悉．其中一個粒子束縛在二維方勢阱中的問題，可以通過二維不含時的薛定諤方程[22-24](6)得出．

(6)

其中V在(-a≤x≤a，-a≤y≤a)范圍內(nèi)取值是0，在(-a≤x≤a，-a≤y≤a)范圍之外取值是+∞，得到方程的本征波函數(shù)和本征值[25-26]分別是：

(7)

圖2 基態(tài)粒子在方勢阱中的概率 密度分布圖(解析解)Fig.2 Probability density distribution of ground state particles in square potential well (analytical solution)

2.2 差分方法計算方勢阱中粒子的概率密度

差分方法是解微分方程的數(shù)值方法，具體而言，是把微分用有限差分代替，把導(dǎo)數(shù)用有限差商代替的方法．它將基本方程和邊界條件近似的用差分方程表示，將求解微分方程的問題改成求解代數(shù)方程的問題．下文用差分方法來求解無限深方勢阱中粒子的概率密度．在笛卡爾坐標(biāo)系中，方勢阱的取值范圍是(-1≤x≤1，-1≤y≤1)，將圖像中x方向和y方向平均分成42份，總共的節(jié)點為1 764個．其中取每份的dx=dy=0.048 7．通過解二維方勢阱中的薛定諤方程，得到粒子在不同能量下的概率密度分布圖，如圖3．

圖3 不同能量的粒子在方勢阱中的概率密度分布圖(差分法)Fig.3 Probability density distribution of particles with different energies in square potential well (difference method)

在方勢阱中，本文對從差分方法和解析解方法中計算出的能量進(jìn)行對比．表1可以看出，差分方法有較高的精度，即使到了第三激發(fā)態(tài)，誤差也僅僅為0.4%．

表1 基于方勢阱中差分方法與解析解之能量誤差Tab.1 Energy errors between the difference method and analytical solution based on square potential well

2.3 網(wǎng)格化對有限元方法精度的影響

劃分勢阱是有限元方法的核心步驟，選擇合理的網(wǎng)格化可以使得有限元方法需要更少的計算卻有更高的精確度．將無限深方勢阱網(wǎng)格化，以圖4為例，取邊心距a=1，a被平均截成的Na份(此處取Na=19)，其中每一份的長度da=0.052 6．圓心角2π被平均分成Nθ份(此處取Nθ=40)，其中每一份的度數(shù)dθ=0.078 5，見圖4．

圖4 方勢阱的網(wǎng)格化圖像Fig.4 Grid image of a square potential well

下文通過對Na和Nθ的改變，來觀察何時粒子的能量趨近于穩(wěn)定，以粒子在方勢阱中的基態(tài)能量為例．在計算中，令Nθ不變，取Nθ為321，Na從10開始，間隔為5，直至增加到85．網(wǎng)格化圖像中的交點也從3 210個開始增加，之后交點的個數(shù)分別為4 815，6 420，8 025，9 630，11 235，12 840，14 445，16 050，17 655，19 260，20 865，22 470，25 680，27 285．如圖5，當(dāng)Na開始增加時候，有限元方法解出的基態(tài)能量迅速下降，直到Na增加至85，基態(tài)能量仍然沒有保持穩(wěn)定．接下去的計算中，令Na不變，并且取Na為40，Nθ從50開始，以20為間隔，直至增加到300．網(wǎng)格化圖像的交點從3 240個開始增加，之后交點的個數(shù)分別為4 840，6 440，8 040，9 640，11 240，12 840，14 440，16 040，17 640，19 240，20 840，22 440，24 040．如圖6，圖6中Nθ數(shù)量從50開始增加，當(dāng)數(shù)量達(dá)到100和150時，計算得到的基態(tài)能量不隨著Nθ增加而改變精度，當(dāng)Nθ達(dá)到220之后，基態(tài)能量保持穩(wěn)定，基本不隨Nθ的變化而變化．

圖5 基態(tài)能量隨log10Na變化關(guān)系曲線 Fig.5 Curve of ground state energy

圖6 基態(tài)能量 隨log10Nθ變化關(guān)系曲線 Fig.6 Curve of ground state energy

對比圖5與圖6，有限元方法計算出的基態(tài)能量與Na成正相關(guān)程度大于其計算出的基態(tài)能量與Nθ的正相關(guān)程度．Nθ對有限元方法計算基態(tài)能量影響較小，Nθ增加了250，數(shù)值計算僅僅只改變0.000 607．相比之下Na增加了75，有限元方法的計算基態(tài)能量的結(jié)果從2.540 365變化為2.048 552，變化量為0.49 1813．

2.4 有限元方法求解粒子在方勢阱中的概率密度

本文用有限元方法得出方勢阱中基態(tài)粒子概率密度分布圖(圖7)與利用解析解得到的粒子在方勢阱中基態(tài)粒子概率密度分布圖(圖2)進(jìn)行比較．在圖7和圖2中，可見，當(dāng)坐標(biāo)和同時等于0時，z有最大的取值，即方勢阱中粒子被測量到的概率最大；隨著坐標(biāo)|x|和|y|的增大，在z方向的數(shù)值減小，即方勢阱中粒子被測量到的概率減小．

圖7 基態(tài)粒子在方勢阱中的概率 密度分布圖(有限元方法)Fig.7 Probability density distribution of ground state particles in square potential well (finite element method)

表2給出了有限元方法和解析解[27-29]的誤差分析．由表2可知，從基態(tài)到第三激發(fā)態(tài)，有限元方法與解析解誤差保持在4.04%到4.14%之間．有限元方法相對于差分法比較，取大致相同的格點數(shù)量，但是它的精度卻遠(yuǎn)遠(yuǎn)不如差分法來的準(zhǔn)確(差分方法與解析解誤差保持在0.05%到0.40%)．在對比中可以看見差分方法與有限元方法求得的第一激發(fā)態(tài)(簡并)圖像不同，那是因為簡并態(tài)的波函數(shù)是兩個“基矢”波函數(shù)的疊加，計算機取的系數(shù)不同，數(shù)值模擬的圖像也不相同，這種情況也只會出現(xiàn)在簡并態(tài)中．但是有限元方法的優(yōu)點在于它可以解出不同邊界的勢阱中粒子的概率密度．在圖8、圖9中，本文利用有限元方法數(shù)值計算出了正多角形勢阱中不同能量下粒子概率密度．甚至正五角形(非凸)勢阱中粒子概率密度分布圖(圖10)也可以用有限元方法求解．

表2 基于方勢阱中有限元方法與解析解之能量誤差Tab.2 Energy errors between finite element method and analytical solution based on square potential well

圖9 不同能量的粒子在五角形勢阱(凸)中的概率密度分布圖(有限元方法)Fig.9 Probability density distribution of particles with different energies in pentagonal well (convex) (finite element method)

圖10 不同能量的粒子在五角形勢阱(非凸)中的概率密度分布圖(有限元方法)Fig.10 Probability density distribution of particles with different energies in pentagonal situation well (nonconvex) (finite element method)

3 粒子在多角形勢阱中的概率密度

3.1 粒子在正三角形勢阱中的解析解

(8)

為了避免方程中粒子的質(zhì)量與量子數(shù)出現(xiàn)混淆，文中將粒子的質(zhì)量用μ代替，且令a=1．在(8)中m和n只能取大于0的整數(shù)，并且在取值上有m≥2n的限制．另外當(dāng)m=2n時，粒子處于非簡并態(tài)，例如當(dāng)m=2，n=1時，粒子在正三角形勢阱中處于非簡并的基態(tài)．粒子在正三角形勢阱中波函數(shù)為：

(9)

通過用matlab計算，經(jīng)過平移后，繪制出正三角形勢阱中基態(tài)粒子解析解圖像，如圖11所示.

從圖11可以看出，當(dāng)x和y的取值趨近于0的時候，概率密度在z方向上的取值增大．當(dāng)x=y=0時，概率密度在z方向上有最大值，即粒子在正三角形勢阱中出現(xiàn)的概率最大；隨著坐標(biāo)|x|和|y|的增大，在z方向上的數(shù)值減小，粒子出現(xiàn)的概率逐漸減?。划?dāng)?shù)搅藙葳宓倪吔缟蠒r，粒子出現(xiàn)的概率幾乎為0．

圖11 基態(tài)粒子在正三角形勢阱中的 概率密度分布圖(解析解)Fig.11 Probability density distribution of ground state particles in a positive triangular potential well (analytic solution)

3.2 有限元方法計算粒子在正三角形勢阱中的概率密度

有限元方法解得正三角形勢阱中基態(tài)粒子的概率密度分布圖，如圖12所示.

圖12 基態(tài)粒子在正三角形勢阱中 的概率密度分布圖(有限元方法)Fig.12 Probability density distribution of ground state particles in positive triangular potential well (finite element method)

將有限元方法得出的正三角形勢阱中基態(tài)粒子的概率密度分布圖(圖12)與利用解析解得到的粒子在正三角形勢阱中基態(tài)粒子的概率密度分布圖(圖11)進(jìn)行比較．當(dāng)坐標(biāo)x和y同時等于0時，兩圖像在方向都取值最大，即正三角形勢阱中粒子被測量到的概率最大；隨著坐標(biāo)|x|和|y|的增大，在z方向的數(shù)值減小；當(dāng)?shù)搅藙葳宓倪吔缟蠒r，粒子出現(xiàn)的概率幾乎為0．粒子在正三角形勢阱中其他能量下的概率密度分布圖見圖13．

圖13 不同能量的粒子在三角形勢阱中的概率密度分布圖(有限元方法)Fig.13 Probability density distribution of particles with different energies in triangular potential well (finite element method)

3.3 五角形的網(wǎng)格化方法與有限元方法計算粒子在五角形勢阱的概率密度

邊界條件對于二維薛定諤方程的解的具體形式有本質(zhì)的影響，在凸五角形勢阱與非凸的五角形勢阱中，粒子的概率密度是難以用解析解或者用差分法求解的．但是，有限元計算為這類問題提供了解決辦法．以下是本文網(wǎng)格化五角形勢阱(凸、非凸)的方式，如圖14、圖15．

圖14 五角形勢阱(凸)的網(wǎng)格化圖像Fig.14 Grid image of pentagonal well(convex)

圖15 五角形勢阱(非凸)的網(wǎng)格化圖像Fig.15 Grid image of pentagonal (nonconvex) well

在圖14中，取邊心距a=1，a被平均截成Na份(此處Na=19)，其中每一份的長度da=0.052 6．圓心角2π被平均分成Nθ份(此處Nθ=50)，其中每一份dθ=0.062 8．在圖15中，取邊心距a=1，a被平均截成Na份(此處Na=19)，其中每一份的長度da=0.052 6．圓心角2π被平均分成Nθ份(此處Nθ=50)，其中每一份dθ=0.062 8．有限元方法計算粒子在五角形勢阱(凸、非凸)中基態(tài)和激發(fā)態(tài)的概率密度分布圖如圖9、圖10所示．

4 結(jié)束語

本文對不同邊界勢阱中的粒子概率密度分布圖實現(xiàn)了可視化，計算出不同邊界勢阱中粒子的能量，通過將有限元方法與解析解的對比，證實了有限元方法的可行性，為有限元方法和量子力學(xué)的結(jié)合提供了重要的實例．與差分方法相比，雖然有限元方法的精度略低，但是它能計算除了方勢阱之外的其他邊界形貌的勢阱，這是非常有價值的．在今后的工作中，將利用有限元方法去處理更多粒子在不同形狀勢阱中的概率密度問題，以及將這種方法應(yīng)用到其他量子力學(xué)問題．為豐富量子力學(xué)研究方法，探索更多有趣且有用的量子現(xiàn)象鋪平了道路．