南郵計算物理實踐報告

1、南 京 郵 電 大 學(xué)實驗 報 告課程名稱： 計算物理實踐 專 業(yè)： 應(yīng)用物理學(xué) 學(xué) 號： 姓 名： 完成日期： 2014 年7月 目 錄第1章 簡單物理實驗的模擬及實驗數(shù)據(jù)處理11.1問題描述11.2原理分析1 1.2.1特殊情況1 1.2.2一般情況31.3Matlab程序仿真41.4Matlab仿真結(jié)果4第2章 方程組的解52.1問題描述52.2原理分析5 2.2.1迭代公式的建立及其幾何意義5 2.2.2解題過程52.3流程圖62.4Matlab程序仿真62.5Matlab仿真結(jié)果6第3章 靜電場問題的計算73.1問題描述73.2原理分析73.3Matlab程序仿真93.4Matlab

2、仿真結(jié)果9第4章 熱傳導(dǎo)方程和波動方程的差分解法104.1問題描述104.2原理分析104.3解題步驟134.4Matlab程序仿真134.5Matlab仿真結(jié)果13第5章 矩量法在靜電場邊值問題計算中的應(yīng)用165.1問題描述165.2原理分析165.3Matlab程序仿真185.4Matlab仿真結(jié)果18結(jié)束語19參考文獻(xiàn)20附錄一21附錄二22附錄三23附錄四25附錄五26第一章 簡單物理實驗的模擬及實驗數(shù)據(jù)處理1.1問題描述模擬電偶極子的場和等位線。設(shè)在處有電荷，在處有電荷。那么在電荷所在平面上任何一點的電勢和場強分別為，。其中，。又設(shè)電荷，。1.2原理分析電偶極子是指一對等值異號的點電

3、荷相距一微小距離所構(gòu)成的電荷系統(tǒng)，它是一種常見的場源存在形式。1.2.1特殊情況圖（1）表示中心位于坐標(biāo)系原點上的一個電偶極子，它的軸線與Z軸重合，兩個點電荷q 和-q 間的距離為L。此電偶極子在場點 P 處產(chǎn)生的電位等于兩個點電荷在該點的電位之和，即 （1）其中與分別是q 和-q 到 P 點的距離。圖（1） 電偶極子一般情況下，我們關(guān)心的是電偶極子產(chǎn)生的遠(yuǎn)區(qū)場，即負(fù)偶極子到場點的距離r 遠(yuǎn)遠(yuǎn)大于偶極子長度L的情形，此時可以的到電偶極子的遠(yuǎn)區(qū)表達(dá)式 （2）可見電偶極子的遠(yuǎn)區(qū)電位與成正比，與r的平方成反比，并且和場點位置矢量r與z軸的夾角有關(guān)。為了便于描述電偶極子，引入一個矢量，模為qL ，方向

4、由-q 指向q ，稱之為此電偶極子的電矩矢量，簡稱為偶極矩，記作 （3）此時（2）式又可以寫成 （4）電偶極子的遠(yuǎn)區(qū)電場強度可由（4）式求梯度得到。因電位 只是坐標(biāo)r 和的函數(shù)，于是有 （5）從（4）式和（5）式可以看到，電偶極子的遠(yuǎn)區(qū)電位和電場分別與r的平方和r的三次方成反比。因此，其電位和場強隨距離的下降比單個點電荷更為迅速，這是由于兩個點電荷q和-q的作用在遠(yuǎn)區(qū)相互抵消的緣故。根據(jù)（4）式，電偶極子的等電位面方程可由 為定值得到。將電力線微分方程寫成球坐標(biāo)形式，并注意此時電場只有r和兩個分量，則有： （6）把電場表達(dá)式（5）帶入上式，得： （7）解上式得： （8）式（8）即是電偶極子遠(yuǎn)區(qū)

5、場的電力線方程。圖（2）繪出了電偶極子為常數(shù)的平面內(nèi)（8）式取不同的常數(shù)所對應(yīng)的等電位線和電場線。 圖（2） 電偶極子的場與等電位線說明：圖中準(zhǔn)確的只是電力線的形狀，電力線的疏密并不嚴(yán)格與場強成正比，只是疏的地方場強小些，密的地方場強大些而已。1.2.2一般情況前面討論了電偶極子的中點位于坐標(biāo)系原點且偶極矩方向為Z的情況。對于中點不在原點和偶極矩非Z的方向的一般情況，通過與前面類似的推導(dǎo)，可以得到遠(yuǎn)區(qū)的電位： （9）其中，r是電偶極子中心指向場點P的相對單位位置矢量，偶極矩P=qL，L的方向依然規(guī)定為從-q到q 。經(jīng)推導(dǎo)還可得到遠(yuǎn)區(qū)場的電場強度表達(dá)式： （10）由上式可以看出，電偶極子的電場線

6、均分布于由r、構(gòu)成的平面上，并且任意一個平面上的電場線分布都相同。從以上幾種不同情況下電偶極子在空間激發(fā)的電場結(jié)果來看，電場強度與p 成正比，與源點到場點的距離成反比，電偶極子在遠(yuǎn)處的性質(zhì)是由其電偶極矩來表征的，電偶極矩是電偶極子的重要特征。設(shè)電荷所在平面上任意一點的電勢為 （11）其中 （12）因此，只要給定空間任意一點的位置坐標(biāo)P（x，y），就可以算出這一點的電位。1.3Matlab程序設(shè)計仿真源程序見附錄一1.4Matlab仿真結(jié)果第二章 方程組的解法2.1問題描述用牛頓法解方程，精度自設(shè)。2.2原理分析2.2.1迭代公式的建立及其幾何意義（1）建立公式將在點Taylor展開Taylor

7、展開線性化近似于解出x記為，則 （n=0,1,2）（2） 幾何意義過切線與求交點，解出，則2.2.2解題過程令，有，那么根據(jù)Newton迭代法建立迭代公式NY開始x0=0.5e=0.0001結(jié)束x-x0e輸出x0=x+2*e2.3流程圖2.4Matlab程序設(shè)計仿真源程序見附錄二2.5Matlab仿真結(jié)果x=0.5671第三章 靜電場問題的計算3.1問題描述長直接地金屬槽，如圖3-2所示，其側(cè)壁和底面電位為零，頂蓋電位為，求槽內(nèi)電位，并繪出電位分布圖。 3.2原理分析（1）原理分析：二維拉普拉斯方程 （1）有限差分法的網(wǎng)格劃分，通常采用完全有規(guī)律的分布方式，這樣可使每個離散點上得到相同形式的差

8、分方程，有效的提高解題速度，經(jīng)常采用的是正方形網(wǎng)格劃分。 設(shè)網(wǎng)格節(jié)點（i,j）的電位為，其上下左右四個節(jié)點的電位分別為在h充分小的情況下，可以為基點進(jìn)行泰勒級數(shù)展開：把以上四式相加，在相加的過程中，h的所有奇次方項都抵消了。得到的結(jié)果的精度為h的二次項。 （2）由于場中任意點都滿足泊松方程：式中為場源，則式（2）可變?yōu)椋?（3）對于無源場，則二維拉普拉斯方程的有限差分形式為： （4）上式表示任一點的電位等于圍繞它的四個等間距點的電位的平均值，距離h越小則結(jié)果越精確，用式（4）可以近似的求解二維拉普拉斯方程。邊界條件： （2）解題過程：在直角坐標(biāo)系中，金屬槽中的電位函數(shù)滿足拉普拉斯方程：其邊界條

9、件滿足混合型邊值問題的邊界條件：取步長，方向上的網(wǎng)格數(shù)為，共有160個網(wǎng)孔和個節(jié)點，其中槽內(nèi)的節(jié)點（電位待求點）有個，邊界節(jié)點52個，設(shè)迭代精度為，利用MATLAB編程求解。3.3Matlab程序設(shè)計仿真源程序見附錄三3.4Matlab仿真結(jié)果第四章 熱傳導(dǎo)方程和波動方程的差分解法4.1問題描述求有限空間內(nèi)的熱傳導(dǎo)問題：的數(shù)值解，邊界條件如教材中圖9.2所示，其他參數(shù)可以自取，將計算結(jié)果圖形化。4.2原理分析二維熱傳導(dǎo)方程的初、邊值混合問題與一維的類似，在確定差分格式并給出定解條件后，按時間序號分層計算，只是每一層是由二維點陣組成，通常稱為網(wǎng)格。內(nèi)部無熱源均勻介質(zhì)中二維熱傳導(dǎo)方程為： ( )

10、(1)其初始條件為： （2）現(xiàn)在設(shè)時間步長為，空間步長為,如圖9.3所示，將平面均分為的網(wǎng)格，并使 則有： 對節(jié)點，在時刻（即時刻）有： （3）將差分格式（3）代入偏微分方程（1）中，可得： （4）式中式（4）為二維熱傳導(dǎo)方程的顯式差分格式，運用式（4）和邊界條件就可以由初始條件逐次計算出任意時刻溫度的分布。下面討論邊界條件：如圖9.3所示陰影部分，即在邊界的和區(qū)域以及整個，邊界均為絕熱壁；而在邊界的區(qū)域為與恒溫?zé)嵩聪噙B的口。和兩邊界溫度始終為0，實際上也是與恒溫源相連的。也就是說，對于絕熱壁應(yīng)滿足： （ ） ( )上述邊界條件的差分近似式為：即： （ ） ( ) （5）對于與恒溫源相連的邊界

11、，在熱傳導(dǎo)過程中始終有恒定的熱流，?？扇w一化值，例如高溫?zé)嵩纯扇　?”，而低溫?zé)嵩纯扇　?”。按圖9.3的情況，邊界條件還有： 綜合上述初值、邊值混合問題，并設(shè)初始時刻各點溫度均為零，則上述差分格式可歸納為： （6）可以證明，對于二維熱傳導(dǎo)方程，若滿足則差分格式式（4）或式（6）就是穩(wěn)定的差分格式，一般的講，對于n維拋物線型微分方程差分格式穩(wěn)定的充分條件是：4.3解題步驟1. 給定、和以及和，題目中已知，的值分別取0s,10s,100s，120s，150s，200s和1000s，和取18和16；2. 計算為36；為32；為0.05；的上界；3. 計算初值和邊值：；4. 用差分格式計算；4.4

12、Matlab程序設(shè)計仿真源程序見附錄四4.5Matlab仿真結(jié)果通過Matlab畫出0s 到1000s 之間的一些溫度場的分布圖，如下圖4.1圖4.7分別為0s,10s,100s,120s,150s,200s,1000s的溫度場分布圖。結(jié)論：很明顯可以看出，溫度呈整體下降的趨勢。由于低溫?zé)嵩吹姆秶雀邷責(zé)嵩吹母?，所以熱量的流入大于流出?？梢詳喽?，只要時間足夠長，整個溫度場除高溫?zé)嵩赐?，其他地方的溫度都要與低溫?zé)嵩聪嗤ㄔO(shè)為0）。1000s 時，如圖4.7所示的場分布與無限長時間之后的場分布就已經(jīng)很接近了。圖4.1 0s時的場分布圖4.2 10s時的場分布圖4.3 100s時的場分布圖4.4

13、120s時的場分布圖4.5 150s時的場分布圖4.6 200s時的場分布圖4.7 1000s時的場分布第五章 矩量法在靜電場邊值問題計算中的應(yīng)用5.1問題描述利用矩量法求無界空間中邊長為2a的正方形導(dǎo)電薄板的電容。5.2原理分析一塊正方形導(dǎo)體板，如上圖所示。邊長為2a米，位于z=0平面，中心坐標(biāo)在原點，設(shè)表示導(dǎo)電板上面電荷密度，板的厚度為零，則空間任意一點的靜電位是 （1）式中，為待求的面電荷密度。邊界條件： （）導(dǎo)體板電容：算子方程：算子：（1） 將導(dǎo)體板分為N個均勻小塊，并選基函數(shù)為分域脈沖函數(shù)。 其中 （2）將式（2）代入式（1）得 m=1,2,3,N （3）式中據(jù)此電荷密度由逼近，平

14、行板電容相應(yīng)地近似為： （4） 若令 表示的邊長，由本身面上的單位電荷密度在其中心處產(chǎn)生的電位是：（2） 用點匹配法選權(quán)函數(shù)為，為的中心點，求內(nèi)積： （5）是處單位均勻電荷密度（）在處中心的電位。 式（5）適用于時的求解，當(dāng)m=n時 （6）其中（3）矩陣求逆解得： 5.3Matlab程序設(shè)計仿真源程序見附錄五5.4Matlab仿真結(jié)果當(dāng)邊長2a=10時，電容C=7.9556e-010由公式推導(dǎo)可知：C的變化和a成正比；有實驗驗證可知：C的變化也和a成正比。結(jié)束語經(jīng)過這次計算物理學(xué)實驗周的學(xué)習(xí)，我認(rèn)識到自己對于以前學(xué)習(xí)過的一些課程掌握得還不夠透徹，Matlab編程語言的運用也不夠熟練。通過這次實

15、驗也很好的鞏固了以前學(xué)習(xí)的一些知識點，并且使我了解了如何利用計算機來模擬和計算一些物理問題。這次實驗讓我認(rèn)識到數(shù)理方程的實用性，掌握了利用差分代替微分來求解波動方程、熱傳導(dǎo)方程、拉普拉斯方程等的基本原理和方法。本次實踐涉及到的二維拉普拉斯方程以及二維熱傳導(dǎo)方程的解題方法，都是先將連續(xù)的方程以及邊界條件離散化，再用計算機進(jìn)行計算，因為計算機智能對離散的數(shù)值進(jìn)行計算。對于非線性方程的求解往往是采用迭代的方法求解，本次實踐主要涉及了Newton迭代法的重要思想，也是將連續(xù)的方程離散化后再進(jìn)行計算。矩量法主要分為三個步驟：（1）離散化；（2）取樣檢測；（2）矩陣求逆；適用于場源分布不確定的情況，用未知

16、場的積分方程來計算給定媒質(zhì)中的場的分布。這次的實踐，使我對Matlab的使用變得熟練了，并且在報告的寫作過程中也熟練掌握了數(shù)學(xué)公式的錄入，文章的排版等技能?？偟膩碚f，這次實踐帶給了我很多的收獲。參考文獻(xiàn)1陳鍾賢.計算物理學(xué).哈爾濱工業(yè)大學(xué)出版社.2001.32楊振華，酈志新.數(shù)學(xué)實驗.科學(xué)出版社.2010.23林亮，吳群英.數(shù)值分析方法與實驗：基于MATLAB實現(xiàn).高等教育出版社.2012.94李慶楊，王能超，易大義.數(shù)值分析.華中科技大學(xué)出版社.2006.75鐘季康，鮑鴻吉.大學(xué)物理習(xí)題計算機解法MATLAB編程應(yīng)用.機械工業(yè)出版社.2008.16何紅雨.電磁場數(shù)值計算法與MATLAB實現(xiàn).

17、華中科技大學(xué)出版社附錄一：close all; clear; clc;k = 9e+9; e_p = 2e-6; e_n = -e_p;d = -10:0.1:10;x, y = meshgrid(d);%產(chǎn)生格點矩陣a=1.5,b=-1.5;x_n = -a; y_n = -b;x_p = a; y_p = b; V1 = k * e_n ./ sqrt(x-x_n).2 + (y-y_n).2); V2 = k * e_p ./ sqrt(x-x_p).2 + (y-y_p).2); V1_min = k * e_n /0.1; V2_max = k * e_p /0.1;V1(V1=-I

18、nf) = V1_min; V1(V1V2_max) = V2_max;V = V1 + V2;E_x, E_y = gradient(-V);hold on; grid on;t=linspace(-pi, pi, 25);px = 0.1 * cos(t) + x_p;py = 0.1 * sin(t) + y_p;streamline(x, y, E_x, E_y, px, py);%畫出電場線sx=min(d)/3*2,min(d),min(d),min(d),min(d)/3*2,min(d),min(d),min(d)/3*1,0,max(d)/3*1,max(d)/3*2;sy=

19、min(d),min(d)/3*1, 0,max(d)/3*1, max(d),max(d)/3*2,max(d),max(d),max(d),max(d),max(d);streamline(x, y, E_x, E_y, sx, sy);%畫出電場線contour(x, y, V, linspace(min(V(:), max(V(:), 180);%畫出等位線 plot(x_n, y_n, ro, x_n, y_n, r-, MarkerSize, 16); plot(x_p, y_p, ro, x_p, y_p, r+, MarkerSize, 16); axis(min(d), ma

20、x(d), min(d), max(d);title(電偶極子的場和等位線);hold off;附錄二：function x=newton(fname,dfname,x0,e)if nargine x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);endtoc附錄三：hx=17;hy=11;%設(shè)置網(wǎng)格v1=ones(hy,hx);%設(shè)置二維數(shù)組for j=1:hx%設(shè)置邊界條件 v1(hy,j)=100*sin(pi*(2*(j-1)/(hx-1);%假設(shè)恰好為一個周期 v1(1,j)=0;end v1(:,1)=0;v2=v1;maxt=1;t=0;k=0

21、;%初始化while(maxt0.00001) %迭代精度 k=k+1;%計算迭代總次數(shù) maxt=0; for i=2:hy-1 for j=2:hx-1 v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1)/4;%拉普拉斯方程差分形式 t=abs(v2(i,j)-v1(i,j); if(tmaxt) maxt=t; end end end v2(2:hy-1,hx)=v2(2:hy-1,hx-1);%右邊界邊界條件 v1=v2;endsubplot(1,2,1),mesh(v2) %3D網(wǎng)格圖axis(0,17,0,14,-20,100)subp

22、lot(1,2,2),contour(v2,16) hold onx=1:1:hx;y=1:1:hy;xx,yy=meshgrid(x,y);Gx,Gy=gradient(v2,0.6,0.6);%計算梯度quiver(xx,yy,Gx,Gy,0.5,r) %根據(jù)梯度畫箭頭axis(-3.5,hx+6.5,-2,15)plot(1,1,hx,hx,1,1,hy,hy,1,1,k)%畫導(dǎo)體框text(hx/2-2,hy+0.6,phi=100sin(pix),fontsize,11);%上標(biāo)注text(hx/2-1,0.5,phi=0,fontsize,11);%下標(biāo)注text(-1.8,hy/2,phi=0,fontsize,11);%左標(biāo)注text(hx+0.2,hy/2,partialphi/partialn=0,fontsize,11);% 右標(biāo)注title(靜電場點位分布圖 );hold off附錄四：N=36;M=32;M1=12;M2=20;D=1;H=0.5;T=0.05;time=10;%初始參數(shù)定義u=zeros(M+1,N+1);%定義場矩陣u(M1+2:M2,1)