數(shù)值法求解靜態(tài)電磁場(chǎng)邊值問(wèn)題

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

1 引言

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

本文采用三角分解法、高斯-賽德?tīng)柕俺沙诘ǎ槍?duì)用有限差分法求解靜態(tài)電磁場(chǎng)邊值問(wèn)題得到的差分方程組進(jìn)行求解。采用不同網(wǎng)格步長(zhǎng)離散化待求區(qū)域，得到規(guī)模不同的差分方程組。采用不同數(shù)值方法求解方程組，可以發(fā)現(xiàn)不同方法有不同的運(yùn)算效率及適用范圍。

2 基于有限差分法的靜態(tài)電磁場(chǎng)邊值問(wèn)題

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

下面我們首先采用一個(gè)簡(jiǎn)單的靜態(tài)電磁場(chǎng)邊值問(wèn)題，來(lái)說(shuō)明有限差分法的求解過(guò)程[5-6]。然后用數(shù)值法求解得到的方程組。

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

由于該槽沿長(zhǎng)度方向?yàn)榫鶆蚍植?，故可以將其?jiǎn)化成二維電磁場(chǎng)問(wèn)題。區(qū)域中無(wú)電荷源，故該區(qū)域中的電位符合二維拉普拉斯分布。設(shè)電位為 ，則有：

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

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

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

3 三角分解法

3.1 計(jì)算原理

三角分解法，又稱(chēng)為L(zhǎng)U分解法。對(duì)線(xiàn)性方程組 ，其中系數(shù)矩陣A可以分解為兩個(gè)矩陣的乘積 ，其中L是單位下三角矩陣，U是上三角矩陣。

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

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

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

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

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

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

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

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

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

4 迭代法

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

4.1 高斯-賽德?tīng)柕?/h3>

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

其中k為迭代次數(shù)，此即高斯-賽德?tīng)柕袷健?/p>

對(duì)于式(2)，可構(gòu)造如下的高斯-賽德?tīng)柕袷剑?/p>

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

4.2 超松弛爾迭代法

把高斯－賽德?tīng)柗ǖ牡底鳛橹虚g值與 加權(quán)求平均，就得到超松弛迭代格式：

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

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

為了研究松弛因子對(duì)迭代收斂速度的影響，在步長(zhǎng)為h=0.1m時(shí)，還比較了時(shí)的計(jì)算時(shí)間迭代收斂時(shí)間分別為112.3315(s)，11.07(s)和68.0965(s)，可見(jiàn)松弛因子的取值對(duì)收斂速度影響很大。同時(shí)也可發(fā)現(xiàn)，在 的最優(yōu)值附近取值時(shí)，超松弛迭代法都比高斯-賽德?tīng)柕ǖ氖諗克俣雀臁?/p>

圖4 步長(zhǎng)為0.01m時(shí)的電位分布圖

5 比較與結(jié)論

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

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

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

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

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

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

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

[6]祝昆.靜態(tài)電磁場(chǎng)邊值問(wèn)題的求解方法[J].六盤(pán)水師范高等專(zhuān)科學(xué)校學(xué)報(bào),2008,06.

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

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

[9]祝昆,邱學(xué)云.基于MATLAB的靜態(tài)場(chǎng)邊值問(wèn)題求解方法[J].文山師范高等專(zhuān)科學(xué)校學(xué)報(bào),2009,02.

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

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