南郵計(jì)算物理實(shí)踐報(bào)告

1、 南 京 郵 電 大 學(xué)實(shí)驗(yàn) 報(bào) 告課程名稱： 計(jì)算物理實(shí)踐 專 業(yè)： 應(yīng)用物理學(xué) 學(xué) 號：姓 名：完成日期： 2014 年7月 23 / 26目 錄第1章 簡單物理實(shí)驗(yàn)的模擬及實(shí)驗(yàn)數(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章 靜電場問題的計(jì)算73.1問題描述73.2原理分析73.3Matlab程序仿真93.4

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

17、 = 9e+9; e_p = 2e-6; e_n = -e_p;d = -10:0.1:10;x, y = meshgrid(d);%產(chǎn)生格點(diǎ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=-Inf) = V1_min; V1(V1<V1_min) = V1_

18、min;V2(V2=Inf) = V2_max; V2(V2>V2_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

19、,max(d)/3*2;sy=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, '

20、ro', x_p, y_p, 'r+', 'MarkerSize', 16); axis(min(d), max(d), min(d), max(d);title('電偶極子的場和等位線');hold off;附錄二：function x=newton(fname,dfname,x0,e)if nargin<4,e=1e-4;endfname=inline('x*exp(x)-1');dfname=inline('exp(x)+x*exp(x)');x0=0.5;x=x0;x0=x+2*e;ticwh

21、ile abs(x0-x)>e 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è)恰好為一個(gè)周期 v1(1,j)=0;end v1(:,1)=0;v2=v1;maxt=1;t=0;k=0;%初始化while(maxt>0.00001) %迭代精度 k=k+1;%計(jì)算迭代總次數(shù) maxt=0; for i=2:hy-1 for j=2:hx

22、-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(t>maxt) 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)subplot(1,2,2),contour(v2,16) hold onx=1:1:hx;y=1:1:hy;某,yy=meshgrid(x,y);Gx

23、,Gy=gradient(v2,0.6,0.6);%計(jì)算梯度quiver(某,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('靜電場點(diǎn)位分布圖 ');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)=