南 南 京 郵 電 大 學 實

　驗 驗 報 告 課程名稱:

　 計算物理實踐

　專 專

　業:

　應用物理學

　 學 學

　號:

　 姓 姓

　名:

　 完成日期:

　2014 年 7 月

　南郵計算物理實踐報告

　 目

　錄 第一章

　簡單物理實驗的模擬及實驗數據處理 „„„„„„„„„1

　1、1 問題描述„„„„„„„„„„„„„„„„„„„„„„„„„1 1、2 原理分析„„„„„„„„„„„„„„„„„„„„„„„„„1

　 1、2、1 特殊情況„„„„„„„„„„„„„„„„„„„„„„1

　 1、2、2 一般情況„„„„„„„„„„„„„„„„„„„„„„3 1、3Matlab 程序仿真„„„„„„„„„„„„„„„„„„„„„„4 1、4Matlab 仿真結果„„„„„„„„„„„„„„„„„„„„„„4 第二章

　方程組的解 „„„„„„„„„„„„„„„„„„„„„„5

　2、1 問題描述„„„„„„„„„„„„„„„„„„„„„„„„„5 2、2 原理分析„„„„„„„„„„„„„„„„„„„„„„„„„5

　 2、2、1 迭代公式的建立及其幾何意義„„„„„„„„„„„„„5

　 2、2、2 解題過程„„„„„„„„„„„„„„„„„„„„„„5 2、3 流程圖„„„„„„„„„„„„„„„„„„„„„„„„„„6 2、4Matlab 程序仿真„„„„„„„„„„„„„„„„„„„„„„6 2、5Matlab 仿真結果„„„„„„„„„„„„„„„„„„„„„„6 第三章

　靜電場問題的計算 „„„„„„„„„„„„„„„„„„„7

　3、1 問題描述„„„„„„„„„„„„„„„„„„„„„„„„„7 3、2 原理分析„„„„„„„„„„„„„„„„„„„„„„„„„7 3、3Matlab 程序仿真„„„„„„„„„„„„„„„„„„„„„„9 3、4Matlab 仿真結果„„„„„„„„„„„„„„„„„„„„„„9 第四章

　熱傳導方程與波動方程的差分解法 „„„„„„„„„„10

　4、1 問題描述„„„„„„„„„„„„„„„„„„„„„„„„„10 4、2 原理分析„„„„„„„„„„„„„„„„„„„„„„„„„10 4、3 解題步驟„„„„„„„„„„„„„„„„„„„„„„„„„13 4、4Matlab 程序仿真„„„„„„„„„„„„„„„„„„„„„„13 4、5Matlab 仿真結果„„„„„„„„„„„„„„„„„„„„„„13 第五章

　矩量法在靜電場邊值問題計算中的應用 „„„„„„„„16

　5、1 問題描述„„„„„„„„„„„„„„„„„„„„„„„„„16 5、2 原理分析„„„„„„„„„„„„„„„„„„„„„„„„„16 5、3Matlab 程序仿真„„„„„„„„„„„„„„„„„„„„„„18 5、4Matlab 仿真結果„„„„„„„„„„„„„„„„„„„„„„18

　結束語 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„19

　參考文獻 „„„„„„„„„„„„„„„„„„„„„„„„„„„„20

　附錄一 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„21

　附錄二 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„22

　附錄三 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„23

　附錄四 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„25

　附錄五 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„26

　南郵計算物理實踐報告 第一章

　簡單物理實驗的模擬及實驗數據處理 1、1 問題描述 模擬電偶極子的場與等位線。

　設在 ) , ( b a 處有電荷 q ? ,在 ) , ( b a ? ? 處有電荷 q ? 。那么在電荷所在平面上任何 一 點 的 電 勢 與 場 強 分 別 為 )1 1(4) , (0 ? ?? ?r rqy x V??, V E ?? ??。

　其 中2 2 2 2) ( ) ( , ) ( ) ( b y a x r b y a x r ? ? ? ? ? ? ? ?? ?,9019 104 ??? ? 。

　又 設 電 荷62 10 q?? ? , 5 . 1 ? a , 5 . 1 ? ? b 。

　1、2 原理分析 電偶極子就是指一對等值異號的點電荷相距一微小距離所構成的電荷系統,它就是一種常見的場源存在形式。

　1、2、1特殊情況 圖(1)表示中心位于坐標系原點上的一個電偶極子,它的軸線與Z軸重合,兩個點電荷q 與-q 間的距離為L。此電偶極子在場點 P 處產生的電位等于兩個點電荷在該點的電位之與,即

　)1 1(4) (0 ? ?? ?r rqr???

　(1) 其中 r ? 與 r ? 分別就是q 與-q 到 P 點的距離。

　圖(1)

　電偶極子

　一般情況下,我們關心的就是電偶極子產生的遠區場,即負偶極子到場點的距離r 遠遠大于偶極子長度L的情形,此時可以的到電偶極子的遠區表達式

　 204c o s) (rqlr???? ?

　(2) 可見電偶極子的遠區電位與 ql 成正比,與r的平方成反比,并且與場點位置矢量r與z軸的夾角β有關。

　為了便于描述電偶極子,引入一個矢量 p ,模為qL ,方向由-q 指向q ,稱之為此電偶極子的電矩矢量,簡稱為偶極矩,記作

　 p ql ?

　(3) 此時(2)式又可以寫成

　 20204 4cos) (rperqlrr?? ???? ? ?

　(4) 電偶極子的遠區電場強度可由(4)式求梯度得到。因電位

　只就是坐標r 與β的函數,于就是有

　 2 30 0cos sin2 4rp pE e er r?? ??? ??? ??? ? ?

　(5) 從(4)式與(5)式可以瞧到,電偶極子的遠區電位與電場分別與r的平方與r的三次方成反比。因此,其電位與場強隨距離 的下降比單個點電荷更為迅速,這就是由于兩個點電荷q與-q的作用在遠區相互抵消的緣故。

　根據(4)式,電偶極子的等電位面方程可由 204cos) (rqlr???? ?

　 為定值得到。

　將電力線微分方程寫成球坐標形式,并注意此時電場只有r與 ? 兩個分量,則有:

　 ??c crErdEdr?

　 (6) 把電場表達式(5)帶入上式,得:

　 ???? ?22sin) (sinsincos2d drdr? ?

　(7) 解上式得:

　南郵計算物理實踐報告 Cr? ?2sin1

　(8) 式(8)即就是電偶極子遠區場的電力線方程。

　圖(2)繪出了電偶極子 ? 為常數的平面內(8)式取不同的常數所對應的等電位線與電場線。

　圖(2)

　電偶極子的場與等電位線 說明:圖中準確的只就是電力線的形狀,電力線的疏密并不嚴格與場強成正比,只就是疏的地方場強小些,密的地方場強大些而已。

　1、2、2一般情況 前面討論了電偶極子的中點位于坐標系原點且偶極矩方向為Z的情況。對于中點不在原點與偶極矩非Z的方向的一般情況,通過與前面類似的推導,可以得到遠區的電位: 204cos) (rqlr???? ?

　 (9) 其中,r就是電偶極子中心指向場點P的相對單位位置矢量,偶極矩P=qL,L的方向依然規定為從-q到q 。

　經推導還可得到遠區場的電場強度表達式:

　03 00304sin4cos 2???????rprrpV E ? ? ?? ?

　(10) 由上式可以瞧出,電偶極子的電場線均分布于由r、θ構成的平面上,并且任意一個平面上的電場線分布都相同。

　從以上幾種不同情況下電偶極子在空間激發的電場結果來瞧,電場強度與p 成正比,與源點到場點的距離3r 成反比,電偶極子在遠處的性質就是由其電偶極矩

　來表征的,電偶極矩就是電偶極子的重要特征。

　設電荷所在平面上任意一點的電勢為

　 )1 1(4) , (0 ? ?? ?r rqy x V??

　(11) 其中

　2 2 2 2) ( ) ( , ) ( ) ( b y a x r b y a x r ? ? ? ? ? ? ? ?? ?

　 (12) 因此,只要給定空間任意一點的位置坐標P(x,y),就可以算出這一點的電位。

　1、3Matlab 程序設計仿真 源程序見附錄一 1、4Matlab 仿真結果

　第二章

　方程組的解法 2、1問題描述 用牛頓法解方程 1 0xxe ? ? ,精度自設。

　2、2原理分析

　2、2、1 迭代公式的建立及其幾何意義 (1)建立公式 將 (x) f 在nx 點 Taylor 展開 """ 2(x )(x) (x ) (x )(x x ) (x x ) ...2!nn n n nff f f ? ? ? ? ? ?

　"(x) (x ) (x )(x x )n n nf f f ? ? ? ——Taylor 展開線性化 (x) 0 f ? 近似于"(x ) (x )(x x ) 0n n nf f ? ? ?

　解出 x 記為1 nx?,則1"(x )(x )nn nnfx xf?? ?

　(n=0,1,2....) （2）幾何意義

　過 ( , ( ))n nx f x 切線 ( ) "( )( )n n ny f x f x x x ? ? ? 與 0 y ? 求交點,解出1 nx x?? ,則1"(x )(x )nn nnfx xf?? ?

　2、2、2 解題過程 令 1 ) ( ? ?xxe x f ,有x xxe e x f ? ? ) ( " ,那么根據 Newton 迭代法建立迭代公式 1"(x ) 1(x )xnn n nx xnf xex x xf e xe??? ? ? ?? 2、3流程圖

　N Y 開始 x0=0、5 e=0、0001 00 00001xx xx ex xe x e?? ?? x-x0>e

　2、4Matlab程序設計仿真 源程序見附錄二 2、5Matlab仿真結果 x=0、5671

　第三章

　靜電場問題的計算 3、1問題描述 長直接地金屬槽,如圖 3-2 所示,其側壁與底面電位為零,頂蓋電位為x ? ? sin 100 ? ,求槽內電位,并繪出電位分布圖。

　 3、2原理分析 (1)原理分析: 二維拉普拉斯方程

　 ) , ( ) , (2y x f y xyy xx? ? ? ? ? ? ?

　 (1) 有限差分法的網格劃分,通常采用完全有規律的分布方式,這樣可使每個離散點上得到相同形式的差分方程,有效的提高解題速度,經常采用的就是正方形網格劃分。

　設 網 格 節點 (i,j) 的電 位 為j i,?, 其 上 下 左右 四 個 節點 的 電 位分 別 為。

　， ， ，j i j i j i j i , 1 , 1 1 , 1 , ? ? ? ?? ? ? ?在 h 充分小的情況下,可以j i,?為基點進行泰勒級數展開: ? ?????????? ??333222, 1 ,6121hyhyhyj i j i? ? ?? ?

　? ?????????? ??333222, 1 ,6121hyhyhyj i j i? ? ?? ? ? ?????????? ??333222, , 16121hxhxhxj i j i? ? ?? ? ? ?????????? ??333222, , 16121hxhxhxj i j i? ? ?? ? 把以上四式相加,在相加的過程中,h 的所有奇次方項都抵消了。得到的結果的精度為 h 的二次項。

　2 22, 1 , 1 1, 1, ,2 24 ( )i j i j i j i j i jhx y? ?? ? ? ? ?? ? ? ?? ?? ? ? ? ? ?? ?

　 (2) 由于場中任意點 ( , ) i j 都滿足泊松方程: 2 222 2= ( , ) F x yx y? ??? ?? ? ?? ? 式中 ( , ) F x y 為場源,則式(2)可變為:

　2, , 1 , 1 1, 1,1( ) ( , )4 4i j i j i j i j i jhF x y ? ? ? ? ?? ? ? ?? ? ? ? ?

　(3) 對于無源場, ( , ) 0 F x y ? ,則二維拉普拉斯方程的有限差分形式為:

　 ) (41, 1 , 1 1 . 1 , , j i j i j i j i j i ? ? ? ?? ? ? ? ? ? ? ? ?

　(4) 上式表示任一點的電位等于圍繞它的四個等間距點的電位的平均值,距離 h越小則結果越精確,用式(4)可以近似的求解二維拉普拉斯方程。

　邊界條件: ( , ) 0(0, ) (a, ) ( ,0) 0( ,b) 100sinxx yyx yy y x Vx xV? ? ?? ? ?? ?? ? ? ? ??? ? ?????

　 (2)解題過程: 在直角坐標系中,金屬槽中的電位函數 ? 滿足拉普拉斯方程: 2 22 20x y? ? ? ?? ?? ? 其邊界條件滿足混合型邊值問題的邊界條件: ( , ) 0(0, ) (a, ) ( ,0) 0( ,b) 100sinxx yyx yy y x Vx xV? ? ?? ? ?? ?? ? ? ? ??? ? ????? 取步長 1 h ? , , x y 方向上的網格數為 16, 10 m n ? ? ,共有160個網孔與17 11 187 ? ? 個節點,其中槽內的節點(電位待求點)有 15 9 135 ? ? 個,邊界節點52個,設迭代精度為610 ? ,利用MATLAB編程求解。

　3、3Matlab程序設計仿真 源程序見附錄三 3、4Matlab仿真結果

　 第四章

　熱傳導方程與波動方程的差分解法 4、1問題描述 求有限空間內的熱傳導問題:2 22 2( , ,0)u u ut x yu x y xy??? ?? ??? ? ????? 的數值解,邊界條件如教材中圖9、2所示,其她參數可以自取,將計算結果圖形化。

　 4、2原理分析 二維熱傳導方程的初、邊值混合問題與一維的類似,在確定差分格式并給出定解條件后,按時間序號分層計算,只就是每一層就是由二維點陣組成,通常稱為網格。

　內部無熱源均勻介質中二維熱傳導方程為: 2 22 2( )u u ut x y?? ? ?? ?? ? ?

　( 0 x l ? ?

　0 y s ? ?

　0 t T ? ? )

　 (1) 其初始條件為: (x,y,0) (x,y) u ? ?

　 (2) 現在設時間步長為 ? ,空間步長為 h ,如圖9、3所示,將 xOy 平面均分為M N ? 的網格,并使 Nh l ?

　 Mh s ? 則有: t k ? ?

　0,1,2.... k ?

　x ih ?

　0,1,2,...,N i ?

　y jh ?

　0,1,2,..., j M ?

　對節點 ( , ) i j ,在 k 時刻(即 k ? 時刻)有: , , , , 1 , ,2, , 1, , , , 1, ,k2 22, , , 1, , , , 1,2 222i j k i j k i j ki j k i j k i j k i ji j k i j k i j k i j ku u utu u u ux hu u u uy h??? ?? ?? ? ?????? ? ? ???????? ?????

　(3) 將差分格式(3)代入偏微分方程(1)中,可得: , , 1 ,j,k 1, , 1, , , 1, , 1,(1 4 ) ( )i j k i i j k i j k i j k i j ku u u u u u ? ?? ? ? ? ?? ? ? ? ? ?

　(4) 式中2h??? ?

　式(4)為二維熱傳導方程的顯式差分格式,運用式(4)與邊界條件就可以由初始條件逐次計算出任意時刻溫度的分布。

　下面討論邊界條件: 如圖9、3所示陰影部分,即在 0 x ? 邊界的10 y M h ? ? 與2M h y Mh ? ? 區域以及整個 x Nh ? , 0 y Mh ? ? 邊界均為絕熱壁;而在 0 x ? 邊界的1 2M h y M h ? ? 區域為與恒溫熱源相連的口。

　0 y ? 與 y Mh ? 兩邊界溫度始終為0,實際上也就是與恒溫源相連的。也就就是說,對于絕熱壁應滿足: 0, ,0j kux??? (1 21,2,..., 1, 1,..., 1 j M M M ? ? ? ?

　 1,2,3,... k ? ) , ,0N j kux???

　( 1,2,..., 1 j M ? ?

　1,2,3,... k ? ) 上述邊界條件的差分近似式為: 1, , 0, ,0j k j ku uh??

　, , 1,j,k0N j k Nu uh???

　即:0, , 1, , j k j ku u ?

　(1 21,2,..., 1, 1,..., 1 j M M M ? ? ? ?

　 1,2,3,... k ? )

　, , 1, , N j k N j ku u??

　 ( 1,2,..., 1 j M ? ?

　1,2,3,... k ? )

　 (5) 對于與恒溫源相連的邊界,在熱傳導過程中始終有恒定的熱流,?？扇w一化值,例如高溫熱源可取“1”,而低溫熱源可取“0”。按圖9、3的情況,邊界條件還有: 0, ,1j ku ?

　 1 2,...,M j M ?

　,0,k , ,0i i M ku u ? ?

　 0,1,2,..., i N ?

　綜合上述初值、邊值混合問題,并設初始時刻各點溫度均為零,則上述差分格式可歸納為: , , 1 , , 1, , 1, , , 1, , 1,, ,00, , 1, , 1 2, , 1, ,,0, , ,(1 4 ) ( )0 0,1,2,..., ; 0,1,2,...,1,2,..., 1, 1,..., 1; 1,2,3,...1,2,..., 1; 1,2,3,...i j k i j k i j k i j k i j k i j ki jj k j kN j k N j ki k i M ku u u u u uu i N j Mu u j M M M ku u j M ku u? ?? ? ? ? ??? ? ? ? ? ?? ? ?? ? ? ? ? ?? ? ? ??0, , 1 20 0,1,2,...,1 ,...,j ki Nu j M M????????? ??? ? ??

　(6) 可以證明,對于二維熱傳導方程,若滿足 214 h??? ? ?

　則差分格式式(4)或式(6)就就是穩定的差分格式,一般的講,對于n維拋物線型微分方程差分格式穩定的充分條件就是: 212 h n??? ? ?

　4、3解題步驟 1.給定 ? 、 h 、 ? 與 T 以及 XN 與 YM ,題目中已知 0.5 h ? ,14? ? , T 的值分別取0s,10s,100s,120s,150s,200s與1000s, XN 與 YM 取18與16; 2.計算XNNh? 為36;YMMh? 為32;2h ???? 為0、05; k 的上界T?; 3.計算初值與邊值:, ,0( , )i ju ih jh ? ? ;0, , 1 (, )j ku g k jh ? ? ;, , 2 (, )N j ku g k jh ? ? ; ,0, 3 (, )i ku g k ih ? ? ;, , 4 (, )i N ku g k ih ? ? ;

　4.用差分格式計算, , 1 i j ku?; 4、4Matlab程序設計仿真 源程序見附錄四 4、5Matlab仿真結果 通過Matlab畫出0s 到1000s 之間的一些溫度場的分布圖,如下圖4、1—圖4、7分別為0s,10s,100s,120s,150s,200s,1000s的溫度場分布圖。

　結論:很明顯可以瞧出,溫度呈整體下降的趨勢。由于低溫熱源的范圍比高溫熱源的更大,所以熱量的流入大于流出?？梢詳喽?只要時間足夠長,整個溫度場除高溫熱源外,其她地方的溫度都要與低溫熱源相同(設為0)。1000s 時,如圖4、7所示的場分布與無限長時間之后的場分布就已經很接近了。

　圖4、1

　 0s時的場分布

　圖4、2

　10s時的場分布

　 圖4、3

　100s時的場分布

　圖4、4

　120s時的場分布

　圖4、5

　150s時的場分布

　 圖4、6

　200s時的場分布

　圖4、7

　 1000s時的場分布

　第五章

　矩量法在靜電場邊值問題計算中的應用 5、1問題描述 利用矩量法求無界空間中邊長為2a的正方形導電薄板的電容。

　5、2原理分析

　一塊正方形導體板,如上圖所示。邊長為 2a 米,位于 z=0 平面,中心坐標在原點,設 ( , ) x y ? 表示導電板上面電荷密度,板的厚度為零,則空間任意一點的靜電位就是

　( x , y , z )0( , )4a aa ad dR? ? ?? ? ???? ?? ??

　(1) 式中2 2 2 1/2[(x ) (y ) z ] R ? ? ? ? ? ? ? , ( , ) ? ? ? 為待求的面電荷密度。

　邊界條件:0(x,y,0) ? ? ?

　( , x a y a ? ? ) 導體板電容:1( , )a aa aqC d dV V? ?? ? ?? ?? ?? ? 算子方程:00( , )( )4a aa ad d LR? ? ?? ? ? ???? ?? ?? ? 算子:014a aa aL d dR? ???? ?? ?? （1）將導體板分為 N 個均勻小塊nS ? ,并選基函數為分域脈沖函數。

　1Nn nnp ? ??? ?

　其中1 S0nnnpS? ????在 上在其它 上

　 (2) 將式(2)代入式(1)得1Nmn nnV l ??? ?

　m=1,2,3,…N

　 (3) 式中2 214 ( ) ( )a amna al d dx y? ??? ? ?? ??? ? ?? ? 據此電荷密度由逼近,平行板電容相應地近似為: 111Nn n mn nn mnC S l SV???? ? ? ?? ?

　(4)

　若令 2 2 / b a N ? 表示的邊長,由nS ? 本身面上的單位電荷密度在其中心處產生的電位就是: 2 21 2(0.8814)4b bmnb bbl dx dyx y????? ?? ??? ? （2）用點匹配法選權函數為 (x x ) (y y )m m mw ? ? ? ? ?, (x ,y )m m為mS ?

　的中心點,求內積: , (p ) (x x ) (y )L(p )dxdymn m n m m nx ay al w L y ? ????? ?? ? ???

　2 21(p )|4 (x ) (y )mn nmn n r rm ml L d d? ?? ??? ? ??? ?? ?? ? ?? ?

　(5) mnl 就是nS ? 處單位均勻電荷密度( 1np ? )在mS ? 處中心的電位。

　0 0, (x x )(y y )m m m mx ay ag w g dxdy ? ? ????? ?? ? ? ???

　 ??011..1mg ?? ?? ?? ?? ? ?? ?? ?? ?? ? 式(5)適用于 m n ? 時mnl 的求解,當 m=n 時 2 2001 2ln(1 2)4b bmnb bbl d d ? ????? ? ?? ?? ? ??? ? (6) 其中22abN?

　(3)矩陣求逆解得: ? ? ? ? ? ?1n mn ml g ??? 1Nn nnp ? ??? ?

　5、3Matlab程序設計仿真 源程序見附錄五 5、4Matlab仿真結果 當邊長 2a=10 時,電容 C=7、9556e-010 由公式推導可知:C 的變化與 a 成正比; 有實驗驗證可知:C的變化也與a成正比。

　 結束語 經過這次計算物理學實驗周的學習,我認識到自己對于以前學習過的一些課程掌握得還不夠透徹,Matlab 編程語言的運用也不夠熟練。通過這次實驗也很好的鞏固了以前學習的一些知識點,并且使我了解了如何利用計算機來模擬與計算一些物理問題。這次實驗讓我認識到數理方程的實用性,掌握了利用差分代替微分來求解波動方程、熱傳導方程、拉普拉斯方程等的基本原理與方法。

　本次實踐涉及到的二維拉普拉斯方程以及二維熱傳導方程的解題方法,都就是先將連續的方程以及邊界條件離散化,再用計算機進行計算,因為計算機智能對離散的數值進行計算。

　對于非線性方程的求解往往就是采用迭代的方法求解,本次實踐主要涉及了Newton 迭代法的重要思想,也就是將連續的方程離散化后再進行計算。

　矩量法主要分為三個步驟:(1)離散化;(2)取樣檢測;(2)矩陣求逆;適用于場源分布不確定的情況,用未知場的積分方程來計算給定媒質中的場的分布。

　這次的實踐,使我對 Matlab 的使用變得熟練了,并且在報告的寫作過程中也熟練掌握了數學公式的錄入,文章的排版等技能。

　總的來說,這次實踐帶給了我很多的收獲。

　參考文獻 [1]陳鍾賢、計算物理學、哈爾濱工業大學出版社、2001、3 [2]楊振華,酈志新、數學實驗、科學出版社、2010、2 [3]林亮,吳群英、數值分析方法與實驗:基于MATLAB實現、高等教育出版社、2012、9

　[4]李慶楊,王能超,易大義、數值分析、華中科技大學出版社、2006、7 [5]鐘季康,鮑鴻吉、大學物理習題計算機解法—MATLAB編程應用、機械工業出版社、2008、1 [6]何紅雨、電磁場數值計算法與MATLAB實現、華中科技大學出版社 附錄一: close all; clear; clc; k = 9e+9;

　e_p = 2e-6;

　e_n = -e_p; d = -10:0、1:10; [x, y] = meshgrid(d);%產生格點矩陣 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_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,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, "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;end fname=inline("x*exp(x)-1"); dfname=inline("exp(x)+x*exp(x)"); x0=0、5; x=x0;x0=x+2*e; tic while abs(x0-x)>e

　 x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); end toc 附錄三: hx=17;hy=11;%設置網格 v1=ones(hy,hx);%設置二維數組 for j=1:hx%設置邊界條件

　 v1(hy,j)=100*sin(pi*(2*(j-1)/(hx-1)));%假設恰好為一個周期

　 v1(1,j)=0; end

　 v1(:,1)=0; v2=v1;maxt=1;t=0;k=0;%初始化 while(maxt>0、00001) %迭代精度

　 k=k+1;%計算迭代總次數

　 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(t>maxt) maxt=t;

　end

　 end

　end

　 v2(2:hy-1,hx)=v2(2:hy-1,hx-1);%右邊界邊界條件

　 v1=v2; end subplot(1,2,1),mesh(v2) %3D 網格圖 axis([0,17,0,14,-20,100]) subplot(1,2,2),contour(v2,16)

　hold on x=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") %根據梯度畫箭頭 axis([-3、5,hx+6、5,-2,15]) plot([1,1,hx,hx,1],[1,hy,hy,1,1],"k")%畫導體框 text(hx/2-2,hy+0、6,"\phi=100sin(\pix)","fontsize",11);%上標注 text(hx/2-1,0、5,"\phi=0","fontsize",11);%下標注 text(-1、8,hy/2,"\phi=0","fontsize",11);%左標注 text(hx+0、2,hy/2,"\partial\phi/\partialn=0","fontsize",11);% 右標注 title("靜電場點位分布圖 "); hold off 附錄四: N=36;M=32;M1=12;M2=20;D=1;H=0、5;T=0、05;time=10;%初始參數定義 u=zeros(M+1,N+1);%定義場矩陣 u(M1+2:M2,1)=ones(M2-M1-1,1);%邊界條件 for i=2:M for j=2:N u(i,j)=(i-1)*H*(j-1)*H;%初始條件 end end u2=u;%差分方程運算開始 for k=1:time/T%k為時間步數

　for i=2:M

　for j=2:N

　u2(i,j)=(1-4*D*T/H/H)*u(i,j)+D*T/H/H*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1));

　 end

　end

　for i=1:M+1

　for j=1:N+1

　u(i,j)=u2(i,j);

　end

　u(i,N+1)=u2(i,N);

　end

　for i=1:M1+1

　u(i,1)=u2(i,2);

　end

　for i=M2+1:M+1

　u(i,1)=u2(i,2);

　end end%差分方程運算結束 mesh(u)%畫圖 xlabel("X Axis"),ylabel("Y Axis"),zlabel("Temperature"),title("Thermal Field Distribution") 附錄五: a=10; N=100; n1=sqrt(N); ltt=ones(N,N); b=a/n1; e1=1e-9; E=1/36/pi*e1; %介電常數 for i=1:n1

　 %獲取各小塊中心坐標

　 for j=1:n1

　 k=n1*(i-1)+j;

　 x(k)=(2*i-1)*b;

　 y(k)=(2*j-1)*b;

　 end end for m=1:N

　 for n=1:N

　if m==n

　 ltt(m,n)=2*b/pi/E*0、8814;

　 else

　 ltt(m,n)=b^2/pi/E/sqrt((x(m)-x(n))^2+(y(m)-y(n))^2);

　 end

　 end end L1=ltt; L2=inv(L1); Lsum=sum(sum(L2)); C=4*b^2*Lsum

推薦訪問： 物理 實踐 計算