1、實(shí)驗(yàn)名稱：用有限差分法解靜電場(chǎng)邊值 時(shí)間：2011/5/7地點(diǎn)：電氣工程學(xué)院 姓名 學(xué)號(hào):20094470 用有限差分法解靜電場(chǎng)邊值問題一、目的1掌握有限差分法的原理與計(jì)算步驟；2理解并掌握求解差分方程組的超松弛迭代法，分析加速收斂因子a的作用；3學(xué)會(huì)用有限差分法解簡(jiǎn)單的二維靜電場(chǎng)邊值問題，并編制計(jì)算程序。二、方法原理有限差分法是數(shù)值計(jì)算中應(yīng)用得最早而又相當(dāng)簡(jiǎn)單、直觀的一種方法。應(yīng)用有限差分法通常所采取的步驟是： 采用一定的網(wǎng)格分割方式離散化場(chǎng)域。 進(jìn)行差分離散化處理。用離散的、只含有限個(gè)未知數(shù)的差分方程組，來近似代替場(chǎng)域內(nèi)具有連續(xù)變量的偏微分方程以及邊界上的邊界條件（也包括場(chǎng)域內(nèi)不同媒質(zhì)分界

2、面上的銜接條件）。 結(jié)合選定的代數(shù)方程組的解法，編制計(jì)算機(jī)程序，求解由上面所得對(duì)應(yīng)于待求邊值問題的差分方程組，所得解答即為該邊值問題的數(shù)值解?，F(xiàn)在，以靜電場(chǎng)邊值問題為例，說明有限差分法的應(yīng)用。f(s)為邊界點(diǎn)s的點(diǎn)函數(shù)，二位場(chǎng)域D和邊界L示于圖5.1-1中。圖5.1-1 有限差分的網(wǎng)格分割1 離散化場(chǎng)域應(yīng)用有限差分法時(shí)，首先需從網(wǎng)格劃分著手決定離散點(diǎn)的分布方式。通常采用完全有規(guī)律的方式，這樣在每個(gè)離散點(diǎn)上可得出相同形式的差分方程，有效地提高解題速度。如圖5.1-1所示，現(xiàn)采用分別與x，y軸平行的等距（步距為h）網(wǎng)格線把場(chǎng)域D分割成足夠多的正方形網(wǎng)格。各個(gè)正方形的頂點(diǎn)（也即網(wǎng)格線的交點(diǎn)）稱為網(wǎng)格

3、的結(jié)點(diǎn)。這樣，對(duì)于場(chǎng)域內(nèi)典型的內(nèi)結(jié)點(diǎn)0，它與周圍相鄰的結(jié)點(diǎn)1、2、3和4構(gòu)成一個(gè)所謂對(duì)稱的星形。2差分格式造好網(wǎng)格后，需把上述靜電場(chǎng)邊值問題中的拉普拉斯方程（1）式離散化。設(shè)結(jié)點(diǎn)0上的電位值為j0。結(jié)點(diǎn)1、2、3和4上的電位值相應(yīng)為j1、j2、j3和j4，則基于差分原理的應(yīng)用，拉普拉斯方程（1）式在結(jié)點(diǎn)0處可近似表達(dá)為j1+j2+j3+j4-4j1=0 （3）這就是規(guī)則正方形網(wǎng)格內(nèi)某點(diǎn)的電位所滿足的拉普拉斯方程的差分格式，或差分方程。對(duì)于場(chǎng)域內(nèi)的每一個(gè)結(jié)點(diǎn)，關(guān)系式（3）式都成立，都可以列出一個(gè)相同形式的差分方程。但是，對(duì)于近鄰邊界的結(jié)點(diǎn)，其邊界不一定正好落在正方形網(wǎng)格的結(jié)點(diǎn)上，而可能如圖5.1

4、-2所示。其中1、2為邊界線上的結(jié)點(diǎn)，p、q為小于1的正數(shù)。仿上所述，可推得對(duì)這些近鄰邊界結(jié)點(diǎn)的拉普拉斯方程的差分格式為 （4）式中：1和2分別是給定邊界條件函數(shù)f (s)在對(duì)應(yīng)邊界點(diǎn)處的值，是已知的。圖5.1-2 近鄰邊界的結(jié)點(diǎn)3邊界條件的近似處理 為了求解給定的邊值問題，還必須對(duì)邊界條件，以及具體問題中可能存在的分界面上的銜接條件，進(jìn)行差分離散化處理，以構(gòu)成相應(yīng)的差分邊值問題。這里，我們只考慮正方形網(wǎng)格分割下的邊界條件的近似處理。第一類邊界條件如果網(wǎng)格結(jié)點(diǎn)正好落在邊界L上，因此對(duì)應(yīng)于邊界條件（2）式的離散化處理，就是把點(diǎn)函數(shù)f (s)的值直接賦予對(duì)應(yīng)的邊界結(jié)點(diǎn)。如果邊界L不通過網(wǎng)格分割時(shí)所

5、引進(jìn)的結(jié)點(diǎn)（例如圖5.1-2中的1、2結(jié)點(diǎn)是邊界線L與網(wǎng)格線的交點(diǎn)，并不是網(wǎng)格分割時(shí)所引進(jìn)的網(wǎng)格結(jié)點(diǎn)），那末在緊鄰邊界的結(jié)點(diǎn)的差分格式應(yīng)選用（4）式，這時(shí)，把點(diǎn)函數(shù)f (s)的值直接賦予邊界線L與網(wǎng)格線的交點(diǎn)1和2。第二類邊界條件應(yīng)當(dāng)指出，從實(shí)際電場(chǎng)問題的分析出發(fā)，如圖5.1-3所示，以電力線為邊界的第二類齊次邊界條件是常見的一種情況。 （5）這時(shí)，可沿著場(chǎng)域邊界外側(cè)安置一排虛設(shè)的網(wǎng)格結(jié)點(diǎn)，顯然，對(duì)于邊界結(jié)點(diǎn)0，由于該處，故必有j1=j3，因此相應(yīng)于邊界條件（5）式的差分計(jì)算格式為2j1+j2+j4-4j0=0 （6） 圖5.1-3 第二類齊次邊界的一種情況 圖5.1-4對(duì)稱線上結(jié)點(diǎn)的差分格式

6、同樣，在許多工程問題中，常常能夠判定待求電場(chǎng)具有某些對(duì)稱性質(zhì)，這樣只需要計(jì)算某一對(duì)稱部分的場(chǎng)就能完全決定整個(gè)場(chǎng)的分布。為此，還必須導(dǎo)出位于場(chǎng)的對(duì)稱線上的結(jié)點(diǎn)所滿足的差分計(jì)算格式。以對(duì)稱線與網(wǎng)格結(jié)點(diǎn)相重合為例（見圖5.1-4），設(shè)線為一對(duì)稱線，對(duì)于位于對(duì)稱線上的任一結(jié)點(diǎn)0，由拉普拉斯方程（因?qū)ΨQ性，必有j1=j3）可得相應(yīng)的差分計(jì)算格式是 2j1+j2+j4-4j0=0 （7）媒質(zhì)分界面上的銜接條件在此選取兩種情況進(jìn)行差分離散化的處理。分界面與網(wǎng)格線相重合的情況；設(shè)分界面L與網(wǎng)格線相重合，如圖5.1-5所示，在兩種媒質(zhì)和中電位都滿足拉普拉斯方程。容易導(dǎo)得，兩種媒質(zhì)分界面上銜接條件在結(jié)點(diǎn)0的差分格

7、式為 （8）其中。 圖5.1-5 分界面與網(wǎng)格線相重合 圖5.1-6 分界面L對(duì)網(wǎng)格呈對(duì)角線形態(tài)分界面對(duì)于網(wǎng)格呈對(duì)角線形態(tài)的情況：如圖5.1-6所示，分界面L對(duì)于網(wǎng)格呈對(duì)角線形態(tài)，在兩種媒質(zhì)和中電位j都滿足拉普拉斯方程。容易導(dǎo)得，兩種媒質(zhì)分界面上銜接條件在結(jié)點(diǎn)0的差分格式為 （9）其中 ?？傊?，類似以上的分析處理方法，可以逐個(gè)導(dǎo)得各種類型的邊界條件和銜接條件差分離散化的計(jì)算格式。限于篇幅，在此不再展開。4差分方程組的求解在對(duì)場(chǎng)域D內(nèi)各個(gè)結(jié)點(diǎn)（包括所有場(chǎng)域內(nèi)點(diǎn)和有關(guān)的邊界結(jié)點(diǎn)）逐一列出對(duì)應(yīng)的差分方程，組成差分方程組后，就可選擇一定的代數(shù)解法，以算出各離散結(jié)點(diǎn)上待求的電位值。注意到差分方程組的系數(shù)

8、一般是有規(guī)律的，且各個(gè)方程都很簡(jiǎn)單，包含的項(xiàng)數(shù)不多（最多不超過5項(xiàng)），因此，對(duì)于有限差分法，通常都采用逐次近似的迭代方法求解。在迭代法的應(yīng)用中，為加速迭代解收斂速度，一般采用的是超松弛迭代法。由于編寫計(jì)算機(jī)程序的需要，每一網(wǎng)格結(jié)點(diǎn)的位置由雙下標(biāo)（i，j）予以識(shí)別，如圖5.1-7所示。對(duì)于差分方程（3）式，采用超松弛迭代法（規(guī)定迭代的運(yùn)算順序是：從左下角開始做起，即i小的先做；對(duì)固定的i，j小的先做。），則關(guān)于結(jié)點(diǎn)0迭代到第（n+1）次時(shí)的近似值，應(yīng)由如下迭代公式算得） （10）圖5.1-7 結(jié)點(diǎn)的雙下標(biāo)（i，j）標(biāo)號(hào)式中：a 稱為加速收斂因子，其取值范圍是1a1e-5) k=k+1 maxt=

9、0 for i=2:hy-1 for j=2:hx-1 v2(i,j)=v1(i,j)+alpha/4*(v1(i+1,j)+v1(i,j+1)+v2(i-1,j)+v2(i,j-1)-4*v1(i,j); t=abs(v2(i,j)-v1(i,j); if(tmaxt)maxt=t;end end end v1=v2; end clf subplot(1,2,1),mesh(v2) axis(0,41,0,21,0,100) subplot(1,2,2),contour(v2,15) hold on x=1:1:hx;y=1:1:hy xx,yy=meshgrid(x,y); Gx,Gy=g

10、radient(v2,0.6,0.6); quiver(xx,yy,Gx,Gy,-0.8,r) axis(-1.5,hx+2.5,-2,13) plot(1,1,hx,hx,1,1,hy,hy,1,1,k) text(hx/2,0.3,0v,fontsize,11); text(hx/2-0.5,hy+0.5,100v,fontsize,11); text(hx+0.3,hy/2,0v,fontsize,11); hold off輸出數(shù)據(jù)：010010010100100100100100100100100100100100100100100100100049.45168.6777.42581.

11、9884.58986.17487.1687.7688.08388.18688.08387.7687.1686.17484.58981.9877.42568.6749.4510029.13447.80559.05165.90770.20172.94774.70775.79576.38876.57776.38875.79574.70772.94770.20165.90759.05147.80529.1340019.28134.36545.06652.39557.36260.70862.92464.32565.09865.34565.09864.32562.92460.70857.36252.395

12、45.06634.36519.2810013.62325.3134.45141.24546.14449.59751.95853.48254.33454.60854.33453.48251.95849.59746.14441.24534.45125.3113.623009.901818.826.18631.98836.37339.57841.82943.3144.14844.41944.14843.3141.82939.57836.37331.98826.18618.89.9018007.183613.80419.50224.1527.78330.51332.47233.7834.52834.77234.52833.7832.47230.51327.78324.1519.50213.8047.1836005.02899.7313.86817.32820.09522.21823.76524.81225.41425.61125.41424.81223.76522.21820.09517.32813.8689.735.0289003.20196.21948.911211.19713.05414.49715.56116.28616.7061