數(shù)值法求解靜態(tài)電磁場邊值問題

電子科技大學物理電子學院 薛 冰

1 引言

在科學和工程中的許多問題，常常歸結(jié)為求解方程或方程組。由于計算量大，大多數(shù)這類問題都不能采用手工計算，而需要利用計算機，采用數(shù)值方法進行計算[1-2]。而其中線性方程組求解不但在科學和工程中常常涉及到，而且數(shù)值計算方法其它分支的一些研究也往往歸結(jié)為此類問題[3]。求解線性方程組的數(shù)值方法主要有兩類：直接法和迭代法。在沒有舍入誤差的時候，采用直接法可以精確求得方程組的解。這類方法主要有高斯消去法、三角分解法、平方根法和追趕法等[3]。迭代法采用極限過程來逐步逼近準確解，包括雅克比迭代法、高斯-賽德爾迭代法及超松弛迭代法等。它對計算機內(nèi)存要求低，更適合于求解高階方程組。

本文采用三角分解法、高斯-賽德爾迭代法及超松弛迭代法，針對用有限差分法求解靜態(tài)電磁場邊值問題得到的差分方程組進行求解。采用不同網(wǎng)格步長離散化待求區(qū)域，得到規(guī)模不同的差分方程組。采用不同數(shù)值方法求解方程組，可以發(fā)現(xiàn)不同方法有不同的運算效率及適用范圍。

2 基于有限差分法的靜態(tài)電磁場邊值問題

有限差分法是一種求解微分方程的數(shù)值方法，由于其邏輯清晰，方法簡單，自上世紀五十年代以來得到了廣泛應用[4]。為求解由微分方程定解問題所構(gòu)造的數(shù)學模型，有限差分法將定解區(qū)域（場區(qū)）離散化，并以各離散點上函數(shù)的差商來近似該點的導數(shù)，使待求的微分方程轉(zhuǎn)化為差分方程組。求解差分方程組，即可得到各離散點處的待求函數(shù)值。

下面我們首先采用一個簡單的靜態(tài)電磁場邊值問題，來說明有限差分法的求解過程[5-6]。然后用數(shù)值法求解得到的方程組。

圖1所示的一個長直接地金屬矩形槽，其橫截面為正方形D： ，其側(cè)壁和地板均接地，電位為0，頂蓋與側(cè)壁絕緣，電位為V0，求該區(qū)域中的電位分布。

由于該槽沿長度方向為均勻分布，故可以將其簡化成二維電磁場問題。區(qū)域中無電荷源，故該區(qū)域中的電位符合二維拉普拉斯分布。設(shè)電位為 ，則有：

圖1 靜態(tài)電磁場邊值問題網(wǎng)格劃分

要采用有限差分法求解，首先將該區(qū)域離散化，沿x和y方向均勻劃分網(wǎng)格。因為求解區(qū)域為正方形，故兩個方向的網(wǎng)格步長都設(shè)為h，每邊的網(wǎng)格數(shù)均為N+1個。區(qū)域內(nèi)的電位用網(wǎng)格點上的電位 來表示，其中分別表示網(wǎng)格在x和y方向的序號，如圖1所示。

(1)式可轉(zhuǎn)化成如下方程組[7]：

3 三角分解法

3.1 計算原理

三角分解法，又稱為LU分解法。對線性方程組 ，其中系數(shù)矩陣A可以分解為兩個矩陣的乘積 ，其中L是單位下三角矩陣，U是上三角矩陣。

根據(jù)矩陣乘法，可以依次求得矩陣U的第k行和L的第k列。原方程組就轉(zhuǎn)化三角方程組：

最后通過回代過程可以方便地求出方程組的解。

3.2 編程思路和計算結(jié)果

采用MATLAB編程求解此問題[8-9]。設(shè)步長為h，則方程組的系數(shù)矩陣為N×N(N=L/h-1)階的方陣。當h較小時，網(wǎng)格數(shù)很多，矩陣規(guī)模就會很龐大。若采用普通的方法來存貯系數(shù)矩陣，則需要N×N階的數(shù)組空間。本問題中，(2)式左邊的系數(shù)矩陣是一個稀疏陣，非零元素很少且為帶狀排列。為節(jié)省內(nèi)存資源，我們采用對角存儲法中的等帶寬存貯法來存放矩陣元素[10-11]。因為對帶狀稀疏矩陣進行LU分解，不會在帶狀結(jié)構(gòu)之外引入非零填入元，因此求解算法不會有任何改變，而且非零元集中存儲，還可以提高計算效率[10]。

對角線存儲法將各非零對角線作為數(shù)組單元進行存儲。對于具有d條非零對角線的N×N稀疏矩陣A，可用兩個數(shù)組表示：一個是N×d的值數(shù)組E，它的每一列存儲了矩陣A的一條對角線元素，列數(shù)d為矩陣A的非零對角線的數(shù)量；另一個是1×d的偏移數(shù)組t=[l1,l2,…,ld]，它依次存儲值數(shù)組E中每一列所對應的對角線相對于主對角線的偏移量，其中在主對角線下方的為負值，上方的為正值。

假設(shè)此問題中，正方形槽的邊長L=1m，采用筆記本電腦（聯(lián)想電腦，CPU為酷睿i 5，主頻2.6GHz，內(nèi)存2GB）進行計算。選取3種步長，分別是h=0.1m，0.05m和0.01m，將計算結(jié)果列于表1。由于每次計算的運行時間有少量偏差，故所列數(shù)值為5次計算的均值。而當網(wǎng)格步長為h=0.01m時，此時步長最小，網(wǎng)格數(shù)量為99×99個，方程組階數(shù)較高，此電腦無法運算，故沒有得到結(jié)果。

圖2給出了網(wǎng)格步長為0.1m和0.05m時的電位分布圖。很顯然，此結(jié)果符合槽中的電位分布情況。

表1 直接法運行時間

圖2 采用直接法計算得到的網(wǎng)格電位分布圖，

4 迭代法

對于線性方程組Ax=b，設(shè)其同解方程組為x=Bx+f，若初始解為x(0)，其迭代格式可以表示為：

4.1 高斯-賽德爾迭代法

對于線性方程組Ax=b，令A=L+D+U，其中，L是矩陣A對角線以下元素構(gòu)成的嚴格下三角矩陣，U是矩陣A對角線以上元素構(gòu)成的嚴格上三角矩陣，D是A的對角元構(gòu)成的矩陣。則可構(gòu)造迭代格式為：

其中k為迭代次數(shù)，此即高斯-賽德爾迭代格式。

對于式(2)，可構(gòu)造如下的高斯-賽德爾迭代格式：

進行迭代，就可計算出待求的網(wǎng)格電位。

4.2 超松弛爾迭代法

把高斯－賽德爾法的迭代值作為中間值與 加權(quán)求平均，就得到超松弛迭代格式：

4.3 迭代法計算結(jié)果

應用MATLAB編程，用與直接法相同的電腦求解。選擇同樣的3種步長，并設(shè)迭代收斂的精度要求為。采用超松弛迭代時，根據(jù)計算得到的松弛因子為。用高斯-賽德爾迭代法和超松弛迭代法計算的迭代次數(shù)和時間（5次運行時間的均值）如表2所示。當網(wǎng)格步長為0.01m時，采用迭代法能夠比較快速地完成計算。采用細網(wǎng)格，可以更準確地模擬空間的電位分布情況，只是運算時間會更長。

為了研究松弛因子對迭代收斂速度的影響，在步長為h=0.1m時，還比較了時的計算時間迭代收斂時間分別為112.3315(s)，11.07(s)和68.0965(s)，可見松弛因子的取值對收斂速度影響很大。同時也可發(fā)現(xiàn)，在 的最優(yōu)值附近取值時，超松弛迭代法都比高斯-賽德爾迭代法的收斂速度更快。

圖4 步長為0.01m時的電位分布圖

5 比較與結(jié)論

本文采用直接法和迭代法對一個靜態(tài)電磁場邊值問題進行求解。從求解情況可知，當矩陣規(guī)模較小時，采用直接法比迭代法的運算效率更高。但當矩陣規(guī)模超過102時，直接法的計算速度明顯降低，且如果矩陣進一步增大，在普通電腦上甚至可能無法計算。而兩種迭代法的求解結(jié)果表明，采用超松弛迭代法改進了高斯-賽德爾迭代法的收斂速度，選取合適的松弛因子更有利于提高收斂速度。

表2 迭代法計算結(jié)果（超松弛迭代法的松弛因子為2）

[1]張飛飛,馬群,黃家慶,佟曉君.求解非線性方程組的二分法[A].科技創(chuàng)新導報,2009.

[2]柳輝.解非線性方程的牛頓迭代法及其應用[A].重慶工學院學報(自然科學版),2007,08.

[3]孫志忠,吳宏偉,袁慰平,聞?wù)鸪?計算方法與實習[M].第5版,東南大學出版社,2011.07.

[4]許秀娟.兩類拋物型方程的有限差分法[D].哈爾濱工業(yè)大學,2009.

[5]宋燎原,王平,張海峰等.靜態(tài)電磁場邊值問題計算方法[J].大學物理,2007,08.

[6]祝昆.靜態(tài)電磁場邊值問題的求解方法[J].六盤水師范高等專科學校學報,2008,06.

[7]王秉中.計算電磁學[M].第1版,科學出版社, 2007.07.

[8]雷亞平,肖洪祥,匡晚成等.基于MATLAB的電磁場數(shù)值分析[J].電子測試,2007,10.

[9]祝昆,邱學云.基于MATLAB的靜態(tài)場邊值問題求解方法[J].文山師范高等專科學校學報,2009,02.

[10]張永杰,孫秦.稀疏矩陣存儲技術(shù)[A].長春理工大學學報,2006,09.

[11]姚遠,劉鵬,王輝,等.基于稀疏矩陣存儲的狀態(tài)表壓縮算法[J].計算機應用,2010,08.