有限差分法模擬一維諧振子

TOC\o"1-5"\h\z\o"CurrentDocument"第1章概述 4\o"CurrentDocument"第2章有限差分方法 5\o"CurrentDocument"2.1有限差分法基本思想 5\o"CurrentDocument"2.2差分方程組的求解 3\o"CurrentDocument"2.2.1高斯-賽德?tīng)柕?32.2.2逐次超松弛法 6\o"CurrentDocument"第3章求解諧振子的微分方程 7\o"CurrentDocument"3.1一維諧振子 7\o"CurrentDocument"3.2二維各向同性諧振子 9\o"CurrentDocument"第4章總結(jié) 10\o"CurrentDocument"參考文獻(xiàn) 11\o"CurrentDocument"附錄 12\o"CurrentDocument"附1一維線性諧振子的程序設(shè)計(jì) 12\o"CurrentDocument"附1.1基態(tài)一維線性諧振子 12\o"CurrentDocument"附1.2第一激發(fā)態(tài)一維線性諧振子 13\o"CurrentDocument"附1.3第二激發(fā)態(tài)一維線性諧振子 13\o"CurrentDocument"附2二維線性諧振子的程序設(shè)計(jì) 14第1章概述微分方程和積分微分方程數(shù)值解的方法?；舅枷胧前堰B續(xù)的定解區(qū)域用有限個(gè)離散點(diǎn)構(gòu)成的網(wǎng)格來(lái)代替，這些離散點(diǎn)稱作網(wǎng)格的節(jié)點(diǎn)；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來(lái)近似； 把原方程和定解條件中的微商用差商來(lái)近似，積分用積分和來(lái)近似，于是原微分方程和定解條件就近似地代之以代數(shù)方程組，即有限差分方程組，解此方程組就可以得到原問(wèn)題在離散點(diǎn)上的近似解。然后再利用插值方法便可以從離散解得到定解問(wèn)題在整個(gè)區(qū)域上的近似解。有限差分法可廣泛用來(lái)求解偏微分方程的近似解，在電磁場(chǎng)中求解點(diǎn)位函數(shù)的拉普拉斯方程時(shí)，可采用有限差分法的基本思想是：用網(wǎng)格將場(chǎng)域進(jìn)行分割，再把拉普拉斯方程用以各網(wǎng)格點(diǎn)處的點(diǎn)位作為未知數(shù)的差分方程式來(lái)進(jìn)行代換，將求解拉普拉斯方程解得問(wèn)題變?yōu)榍舐?lián)立差分方程組的解得問(wèn)題［1］，在差分網(wǎng)格非常多和情況下，利用并行計(jì)算方法對(duì)其進(jìn)行區(qū)域分解，每個(gè)進(jìn)程負(fù)責(zé)運(yùn)算一部分區(qū)域，區(qū)域邊界之間進(jìn)行必要地通信可有效提高計(jì)算速度，解決更大規(guī)模的問(wèn)題。往往只討論它在靜態(tài)場(chǎng)中的應(yīng)用，即泊松方程或拉普拉斯方程的有限差分形式，很少涉及到它在時(shí)諧場(chǎng)(即亥姆霍茲方程)中的應(yīng)用。本文重點(diǎn)討論亥姆霍茲方程的有限差分形式以及它在時(shí)諧場(chǎng)中的應(yīng)用。同時(shí)，有限差分法(finitedifferencemethod)是基于差分原理的一種數(shù)值計(jì)算方法，在求解微分方程定解問(wèn)題中廣泛應(yīng)用。有限差分法是以差分原理為基礎(chǔ)的一種數(shù)值計(jì)算法。它用離散的函數(shù)值構(gòu)成的差商來(lái)近似逼近相應(yīng)的偏導(dǎo)數(shù)，而所謂的差商則是基于差分的應(yīng)用的數(shù)值微分表達(dá)式。用離散的只含有有限個(gè)未知量的差分方程組去近似代替連續(xù)變量的微分方程和定解條件，并把差分方程組的姐作為威風(fēng)方程定解問(wèn)題的近似解.有限差分法可以處理幾乎所有形式的勢(shì)函數(shù)，且主程序不依賴于勢(shì)函數(shù)的具體形式，對(duì)于多數(shù)兩字體都可以進(jìn)行相對(duì)準(zhǔn)確的計(jì)算。因此，將有限差分法應(yīng)用于量子力學(xué)本征值問(wèn)題的計(jì)算，有助于相對(duì)準(zhǔn)確地進(jìn)行量子體系和形象直觀地教學(xué)研究［23］。量子力學(xué)教程中隊(duì)一維無(wú)限深勢(shì)阱、線性諧振子、氫原子等量子體系的薛定諤方程進(jìn)行了嚴(yán)格的求解，得到了描述體系狀態(tài)的波函數(shù)和能量的精確解。多數(shù)量子體系的哈密頓算數(shù)比較復(fù)雜，薛定諤方程不能嚴(yán)格求解，因此，研究和發(fā)展薛定諤方程的數(shù)值計(jì)算方法具有重要意義⑷。

第2章有限差分方法2.1有限差分法基本思想有限差分法是解偏微分方程的主要數(shù)值方法之一，其基本思想是把連續(xù)的問(wèn)題離散化，即首先對(duì)求解區(qū)域作網(wǎng)格剖分，用有限個(gè)網(wǎng)格節(jié)點(diǎn)代替連續(xù)區(qū)域；其次將微分算子離散化，從而把微分方程的定解問(wèn)題化為代數(shù)方程組的求解問(wèn)題，解方程組就可以得到原問(wèn)題在離散點(diǎn)上的近似解。然后再利用插值方法便可以從離散解得到定解問(wèn)題在整個(gè)區(qū)域上的近似解。參照文獻(xiàn)［5］，給出有限差分法數(shù)值計(jì)算的基本思想：（1） 區(qū)域的離散或子區(qū)域的劃分。（2） 插值函數(shù)的選擇。（3） 方程組的建立。（4） 方程組的求解。2.2差分方程組的求解利用有限差分要面臨求解的問(wèn)題。在實(shí)際應(yīng)用中，可以采用逐次迭代的迭代方法求解。這里介紹兩種常用的迭代方法：高斯賽-德?tīng)柕ê椭鸫纬沙诘ā?.2.1高斯-賽德?tīng)柕ㄖ校╧+1）=1中（k+1）+中（k中（k+1）=1中（k+1）+中（k+1）+中（k）+中（k）+h2二i,j4 -(2-1)%I其中i,j=1,2,(2-1)Xi-1,j 'i,jT 'i+1,j 'i,j+1 8J反復(fù)迭代（k=0,1......），一直進(jìn)行到對(duì)所有節(jié)點(diǎn)滿足下列條件為止。(2-2)"以+1）一平以）＜W(2-2)式中，W是預(yù)定的最大允許誤差。在高斯-賽德?tīng)柕?，網(wǎng)格節(jié)點(diǎn)一般按“自然順序”排列，即先“從左到右”，再“從上到下”排列，如圖1-1所示；

圖1-1圖1-1網(wǎng)格節(jié)點(diǎn)排列2.2.2逐次超松弛法則可得將式(2-1)在同一點(diǎn)上相鄰兩次迭代的差值計(jì)為r(n)(i,j)則可得、1R(n)(l、1R(n)(l,j)=中(k+1)—中(k)=?i，j l,j 4—(P(k+1)+中(k+1)+中(k)+中(k)+4"l—1,j i,j-1 l+1,j l,j+1P(k+1)=P(k)+ap(n)a4(2-4)P(k+1)=P(k)+ap(n)a4(2-4)= P(k+1)+P(k+1)+P(k)+P(k)+h2—7 4p(k)4"i—1，ji，j—1i+1，ji，j+1 8i，j|其中a稱為“加速收斂因子”，是一個(gè)供選擇的參數(shù)，其值在1<a<2之間，該方法的快慢與a有著明顯的關(guān)系。實(shí)踐表明，如果a選得好，可以較快的加速迭代的速度。經(jīng)驗(yàn)表明正方形場(chǎng)域由正方形網(wǎng)格劃分，每邊節(jié)點(diǎn)數(shù)為(P+1)時(shí)，最佳的收斂因子為a0(2-5)

第3章求解諧振子的微分方程3.1一維諧振子對(duì)一維線性諧振子，其能量本證方程為h2d2甲1 廣(3-1)+口①2x2Q=Eq(3-1)2曰dx22邊界條件(3-2)去長(zhǎng)度單位為力s'2(3-2)去長(zhǎng)度單位為力s'2引入無(wú)量綱參量&=oxa=以①.:方，則式則式(3-1)化為無(wú)量綱形式令V=qG)，ii式令V=qG)，ii式(3-2)可寫為二qG)=&2dg2q^=Xq(E)式(2-3)可寫為db+g2h-，Ei+i=w(3-3)(3-4)V0=0,V0=0,v=0(3-5)考慮式(3-5)，式(3-4)可寫為SVi以Vi，其中Vm=2VJ、m-1/,2,2、+g2+101wx2（A&*也2211w2 +g（Ag）2 m-21（A&）2S為三對(duì)角矩陣用Matlab編程計(jì)算可同時(shí)得到矩陣S

步長(zhǎng)值時(shí)一維諧振子能量S為三對(duì)角矩陣用Matlab編程計(jì)算可同時(shí)得到矩陣S

步長(zhǎng)值時(shí)一維諧振子能量E的數(shù)值計(jì)m安N=0N=1N=2500.20.99752.98744.96731000.10.99942.99694.99192000.050.99982.99924.9980表1一維諧振子能量（以"2單）計(jì)算結(jié)果比較精確，且計(jì)算精度隨求解區(qū)域內(nèi)格點(diǎn)數(shù)目增加而增高，主量子1 3 5數(shù)為0、1、2時(shí)。能量精確解分別為-方①、-方①、-方①，簡(jiǎn)并度f(wàn)為1.表1」 」 」中N為量子數(shù)，取m=200，△戶0.05時(shí)波函數(shù)圖形如圖3-1所示.表1中N表1中N為量子數(shù)，取m=200，^e=0.05時(shí)波函數(shù)圖形如圖3-1所示，Matlab程序見(jiàn)附1.1，附1.2，附1.3.圖3-1一維諧振子波函數(shù)3.2二維各向同性諧振子對(duì)于二維各向同性諧振子，薛定諤方程為邊界條件一宗言+等)+2岫2邊界條件一宗言+等)+2岫2+y"網(wǎng)XT3,中G,y)=0；yT3,中(x,y)=0(3-6)(3-7)取長(zhǎng)度單位呻，能量單位為力W；2，引入無(wú)量綱參量&=ax,匚=ay,a=寸呻：力，入=E；2力①,則式(3-6)可寫為(3-8)+豈邱,?)+￡*2從,?)=抑&,°)(3-8)令^=中(，°)&=i-—△&,°=j-—A°,1<i<m-1;1<j<n-1,式(3-8)可寫1,j ijiI2JjI2J為：一r八巾-r：\? +12La：\+h八1+(&2+°2)0-? -f人八巾 =叫(A°力 i,j-1 (A&力 i-1,j |_ 1也&力(A°>fe J J i,j 也&力 i+1,j 也°力 i,j+1 i,j(3-9)式(3-7)可寫為80,=0,8.0=0,0,=0,8.=0 (3-10)則在是(2-10)下，式(2-11)可寫為D?.1+C?.+D?.1=人0. (3-11)102\r?102\r?)0,1,j?2,j:，?=:1j?(A°)2/?—-Lj"m-1,jJC=2^^-+二］+(2C=2^^-+二］+(2+C@)1j1一@2)0)_山:晶通P01'2J+0+C2/1j輸 再+@1一也1)(尸+^2+匚2

1j/則式（3-11）可進(jìn)步寫成矩陣式=人e，其中jjfCD0 0)DC00S二：:::00.CDk00.DC)ks為三對(duì)角塊矩陣，是大型稀疏矩陣，用Matlab編程計(jì)算可同時(shí)得到稀疏矩陣s的本征值和本征矢.表2給出了取不同格點(diǎn)、步長(zhǎng)值時(shí)二維諧振子能量E的數(shù)值計(jì)算結(jié)果與精確解的比較，從表2可以得到與表1相似的結(jié)論，主量子數(shù)為0、1、2時(shí)，能量精確解分別為力①、2力①、3力①，簡(jiǎn)并度分別為1、2、3.mn氐N=0N=1N=220200.50.51.96953.90465.777020200.50.5—3.91765.805420200.50.5——5.881040400.250.251.99593.98005.952540400.250.25—4.01455.987040400.250.25——6.116880800.1250.1252.00404.00005.993380800.1250.125—4.05076.043780800.1250.125——6.2166表2二維諧振子的能量（以"2單位）取m=n=20，Ae=A^0.5二維諧振子的波函數(shù)圖形如圖3-3所示，其中橫向的2個(gè)坐標(biāo)分別為無(wú)量綱參量旨?？v向坐標(biāo)為波函數(shù)平。,。），Matlab程序見(jiàn)附2。-0；25圖3-3二維諧振子波函數(shù)第4章總結(jié)計(jì)算物理學(xué)中，有限差分法（簡(jiǎn)稱差分法），它以概念清晰，方法簡(jiǎn)單，直觀的特點(diǎn)應(yīng)用于數(shù)值分析領(lǐng)域，因而應(yīng)用廣泛。無(wú)論是常微分方程還是偏微分方程，各種類型的二階線性方程，以全高階或非線性方程，均可利用差分法轉(zhuǎn)換為代數(shù)方程組，然后利用計(jì)算機(jī)求其數(shù)值解。有限差分法是以差分原理為基礎(chǔ)的一種數(shù)值解法，通過(guò)學(xué)習(xí)計(jì)算物理學(xué)的有限差分法，我們對(duì)很多問(wèn)題都可以得到足夠高的計(jì)算精度。有限差分法是我們較容易掌握的數(shù)值解法，我們?cè)谖锢韺W(xué)中，可以利用有限差分法解任何偏微分方程。利用有限差分法解邊值問(wèn)題時(shí)，首先將求解區(qū)域分為很多個(gè)網(wǎng)格和節(jié)點(diǎn)，并用差商代替微商，然后，使區(qū)域中的偏微分方程轉(zhuǎn)化為以節(jié)點(diǎn)的數(shù)值為未知量的差分方程組，最后，解該方程組便可得到各離散點(diǎn)待求的數(shù)值解。該數(shù)值解是近似解，但逼近區(qū)域的真實(shí)解。在我們的實(shí)驗(yàn)設(shè)計(jì)中，如果離散化的點(diǎn)選擇的足夠密的話，我們得到的數(shù)值解與真實(shí)解的誤差就能減小到可接受的程度。參考文獻(xiàn)[1]吳連坳，井孝功，丁慧明等.徑向薛定諤方程的有限差分解法[J],吉林大學(xué)自然科學(xué)報(bào)，1994，(03):67-70.⑵張志涌.精通MATLAB6.5北京：北京航空航天大學(xué)出版社，2003.劉建軍，翟利學(xué).有限差分法解能量本征方程，北京工業(yè)大學(xué)學(xué)報(bào)，2008,(34)：325-328.曾謹(jǐn)言.量子力學(xué)[M].北京：科學(xué)出版社,，2000.附錄附1一維線性諧振子的程序設(shè)計(jì)附1.1基態(tài)一維線性諧振子clcclearM=200;H=20;h=H/(2*M);R=1./h人2;fori=1:2*M-1forj=1:2*M-1A(i,j)=0;endendforp=1:2*M-1A(p,p)=2*R+((p-M)*h)人2;endfori=1:2*M-2A(i,i+1)=-R;endfori=2:2*M-1A(i,i-1)=-R;end[v,d]=eig(A);d=eig(A);y=v(:,1);z=linspace(-H/2,H/2,2*M-1);y11=y.A2;M2=trapz(z,y11);y12=1/M2*y11;%plot(z,y12);y=(1/M2)A0.5*y;plot(z,y);holdon附1.2第一激發(fā)態(tài)一維線性諧振子clcclearM=200;H=20;h=H/(2*M);R=1./hA2;fori=1:2*M-1forj=1:2*M-1A(i,j)=0;endendforp=1:2*M-1A(p,p)=2*R+((p-M)*h)A2;endfori=1:2*M-2A(i,i+1)=-R;endfori=2:2*M-1A(i,i-1)=-R;end[v,d]=eig(A);d=eig(A);y=v(:,2);z=linspace(-H/2,H/2,2*M-1);y11=y.A2;M2=trapz(z,y11);y12=1/M2*y11;%plot(z,y12);y=(1/M2)A0.5*y;plot(z,y);holdon附1.3第二激發(fā)態(tài)一維線性諧振子clcclearM=200;H=20;h=H/(2*M);R=1./hA2;fori=1:2*M-1forj=1:2*M-1A(i,j)=0;endendforp=1:2*M-1A(p,p)=2*R+((p-M)*h)A2;endfori=1:2*M-2A(i,i+1)=-R;endfori=2:2*M-1A(i,i-1)=-R;end[v,d]=eig(A);d=eig(A);y=v(:,3);z=linspace(-H/2,H/2,2*M-1);y11=y.A2;M2=trapz(z,y11);y12=1