數值法求解靜態電磁場邊值問題

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

1 引言

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

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

2 基于有限差分法的靜態電磁場邊值問題

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

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

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

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

圖1 靜態電磁場邊值問題網格劃分

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

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

3 三角分解法

3.1 計算原理

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

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

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

3.2 編程思路和計算結果

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

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

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

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

表1 直接法運行時間

圖2 采用直接法計算得到的網格電位分布圖，

4 迭代法

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

4.1 高斯-賽德爾迭代法

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

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

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

進行迭代，就可計算出待求的網格電位。

4.2 超松弛爾迭代法

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

4.3 迭代法計算結果

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

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

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

5 比較與結論

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

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

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

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

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

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

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

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

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

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

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

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

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