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

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

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

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

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

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

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

7、展開(kāi)線性化近似于解出x記為，則 （n=0,1,2）（2） 幾何意義過(guò)切線與求交點(diǎn)，解出，則2.2.2解題過(guò)程令，有，那么根據(jù)Newton迭代法建立迭代公式NY開(kāi)始x0=0.5e=0.0001結(jié)束x-x0e輸出x0=x+2*e2.3流程圖2.4Matlab程序設(shè)計(jì)仿真源程序見(jiàn)附錄二2.5Matlab仿真結(jié)果x=0.5671第三章 靜電場(chǎng)問(wèn)題的計(jì)算3.1問(wèn)題描述長(zhǎng)直接地金屬槽，如圖3-2所示，其側(cè)壁和底面電位為零，頂蓋電位為，求槽內(nèi)電位，并繪出電位分布圖。 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)行泰勒級(jí)數(shù)展開(kāi)：把以上四式相加，在相加的過(guò)程中，h的所有奇次方項(xiàng)都抵消了。得到的結(jié)果的精度為h的二次項(xiàng)。 （2）由于場(chǎng)中任意點(diǎn)都滿足泊松方程：式中為場(chǎng)源，則式（2）可變?yōu)椋?（3）對(duì)于無(wú)源場(chǎng)，則二維拉普拉斯方程的有限差分形式為： （4）上式表示任一點(diǎn)的電位等于圍繞它的四個(gè)等間距點(diǎn)的電位的平均值，距離h越小則結(jié)果越精確，用式（4）可以近似的求解二維拉普拉斯方程。邊界條件： （2）解題過(guò)程：在直角坐標(biāo)系中，金屬槽中的電位函數(shù)滿足拉普拉斯方程：其邊界條

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

10、(1)其初始條件為： （2）現(xiàn)在設(shè)時(shí)間步長(zhǎng)為，空間步長(zhǎng)為,如圖9.3所示，將平面均分為的網(wǎng)格，并使 則有： 對(duì)節(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í)際上也是與恒溫源相連的。也就是說(shuō)，對(duì)于絕熱壁應(yīng)滿足： （ ） ( )上述邊界條件的差分近似式為：即： （ ） ( ) （5）對(duì)于與恒溫源相連的邊界

11、，在熱傳導(dǎo)過(guò)程中始終有恒定的熱流，?？扇w一化值，例如高溫?zé)嵩纯扇　?”，而低溫?zé)嵩纯扇　?”。按圖9.3的情況，邊界條件還有： 綜合上述初值、邊值混合問(wèn)題，并設(shè)初始時(shí)刻各點(diǎn)溫度均為零，則上述差分格式可歸納為： （6）可以證明，對(duì)于二維熱傳導(dǎo)方程，若滿足則差分格式式（4）或式（6）就是穩(wěn)定的差分格式，一般的講，對(duì)于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.4

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

13、120s時(shí)的場(chǎng)分布圖4.5 150s時(shí)的場(chǎng)分布圖4.6 200s時(shí)的場(chǎng)分布圖4.7 1000s時(shí)的場(chǎng)分布第五章 矩量法在靜電場(chǎng)邊值問(wèn)題計(jì)算中的應(yīng)用5.1問(wèn)題描述利用矩量法求無(wú)界空間中邊長(zhǎng)為2a的正方形導(dǎo)電薄板的電容。5.2原理分析一塊正方形導(dǎo)體板，如上圖所示。邊長(zhǎng)為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ù)此電荷密度由逼近，平

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

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

16、場(chǎng)的積分方程來(lái)計(jì)算給定媒質(zhì)中的場(chǎng)的分布。這次的實(shí)踐，使我對(duì)Matlab的使用變得熟練了，并且在報(bào)告的寫(xiě)作過(guò)程中也熟練掌握了數(shù)學(xué)公式的錄入，文章的排版等技能。總的來(lái)說(shuō)，這次實(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何紅雨.電磁場(chǎng)數(shù)值計(jì)算法與MATLAB實(shí)現(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)生格點(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=-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);%畫(huà)出電場(chǎng)線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);%畫(huà)出電場(chǎng)線contour(x, y, V, linspace(min(V(:), max(V(:), 180);%畫(huà)出等位線 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(電偶極子的場(chǎng)和等位線);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è)恰好為一個(gè)周期 v1(1,j)=0;end v1(:,1)=0;v2=v1;maxt=1;t=0;k=0

21、;%初始化while(maxt0.00001) %迭代精度 k=k+1;%計(jì)算迭代總次數(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);%計(jì)算梯度quiver(xx,yy,Gx,Gy,0.5,r) %根據(jù)梯度畫(huà)箭頭axis(-3.5,hx+6.5,-2,15)plot(1,1,hx,hx,1,1,hy,hy,1,1,k)%畫(huà)導(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(靜電場(chǎng)點(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);%定義場(chǎng)矩陣u(M1+2:M2,1)