電子科技大學(xué)物理電子學(xué)院 薛 冰

1 引言

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

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

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

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

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

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

由于該槽沿長度方向?yàn)榫鶆蚍植迹士梢詫⑵浜喕啥S電磁場問題。區(qū)域中無電荷源，故該區(qū)域中的電位符合二維拉普拉斯分布。設(shè)電位為 ，則有：

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

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

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

3 三角分解法

3.1 計(jì)算原理

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

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

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

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

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

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

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

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

表1 直接法運(yùn)行時間

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

4 迭代法

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

4.1 高斯-賽德爾迭代法

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

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

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

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

4.2 超松弛爾迭代法

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

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

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

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

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

5 比較與結(jié)論

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

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

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

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

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

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

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

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

[7]王秉中.計(jì)算電磁學(xué)[M].第1版,科學(xué)出版社, 2007.07.

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

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

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

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