湖南大學(xué)電氣與信息工程學(xué)院 林昊宇

基于MATLAB的超松弛迭代法的電位研究

湖南大學(xué)電氣與信息工程學(xué)院 林昊宇

有限差分法作為靜態(tài)場(chǎng)電位求解問題的有效數(shù)值解法，它將場(chǎng)域劃分為有限個(gè)網(wǎng)格節(jié)點(diǎn)，用差分方程組求得近似解。本文由此出發(fā)，詳細(xì)介紹了用有限差分法求解場(chǎng)域電位的基本原理，并詳細(xì)闡述了幾種常用的迭代方法，并用實(shí)例重點(diǎn)分析了運(yùn)用超松弛迭代法求解場(chǎng)域電位的方法。本文運(yùn)用MATLAB語言進(jìn)行計(jì)算，簡(jiǎn)單、易讀，擺脫了傳統(tǒng)C語言程序的不足，最后通過繪制電位分布圖像仿真驗(yàn)證了算法和程序的正確性與有效性。

有限差分法；超松弛迭代法；靜態(tài)場(chǎng)電位問題；MATLAB

一、引言

隨著計(jì)算機(jī)技術(shù)、實(shí)驗(yàn)研究技術(shù)的迅速發(fā)展，電磁場(chǎng)學(xué)科在教研、工程上的應(yīng)用也在逐步加深。而經(jīng)典電磁學(xué)理論作為電磁技術(shù)發(fā)展的理論基礎(chǔ)，其核心麥克斯韋方程組的重要性不言而喻。因此在工程上求解電磁場(chǎng)的基本任務(wù)，便是根據(jù)電磁場(chǎng)域的特性建立數(shù)學(xué)模型，利用麥克斯韋方程列出方程，利用邊界條件，求解出磁場(chǎng)的分布。

電位的計(jì)算則是靜電場(chǎng)計(jì)算的基本問題。只要計(jì)算出電位，包括電場(chǎng)強(qiáng)度在內(nèi)的電磁場(chǎng)物理量都可以由其求得。在以往的區(qū)域電位求解中，主要有解析法和數(shù)值法兩大類。解析法所求得的結(jié)果較為精準(zhǔn)，但是在實(shí)際的工程問題中往往因?yàn)闂l件過于復(fù)雜而無法求得。故在此情況下，用數(shù)值法求解電磁場(chǎng)電位便是良方。近年來，隨著計(jì)算機(jī)技術(shù)的發(fā)展和各種算法技術(shù)的進(jìn)步，數(shù)值方法有了很大的飛躍，主要內(nèi)容包括：有限差分法、有限元法、蒙特卡洛法等。本文將使超松弛迭代法，通過MATLAB研究區(qū)域磁場(chǎng)的電位問題。

二、關(guān)于超松弛迭代法

超松弛迭代法是在有限差分法的基礎(chǔ)上為提高收斂速度，采用加權(quán)平均而得到的新算法。因此我們需要先了解運(yùn)用有限差分法求解區(qū)域磁場(chǎng)的基本原理。所謂有限差分法（finite difference method），是一種求偏微分（或常微分）方程和方程組定解問題的數(shù)值解的方法，它以差分原理為基礎(chǔ)，具有簡(jiǎn)單、直觀的特點(diǎn)，是最早廣泛應(yīng)用于電磁場(chǎng)數(shù)值分析領(lǐng)域的方法。用此方法時(shí)，采用離散化思想，計(jì)算的步驟通常如下：

采用一定的網(wǎng)格劃分格式，把實(shí)際連續(xù)的的磁場(chǎng)離散為有限的多個(gè)點(diǎn)。

使用差分原理，對(duì)磁場(chǎng)域內(nèi)偏微分方程以及磁場(chǎng)域邊界條件進(jìn)行差分離散化處理，即用差商替代偏導(dǎo)數(shù)，寫出相應(yīng)的差分格式。

結(jié)合代數(shù)方程組的求解，可通過計(jì)算機(jī)程序，求出由以上的步驟得到的代求邊界問題的差分方程的解，即可求出網(wǎng)格節(jié)點(diǎn)的位函數(shù)值，進(jìn)一步求得電位分布。

下面以二維泊松方程的第一類邊值問題為例介紹用有限元差分法求解區(qū)域電位的基本原理。

二維靜電場(chǎng)的邊值條件滿足：

將如圖1的場(chǎng)域分成足夠小的正方形網(wǎng)絡(luò)，網(wǎng)絡(luò)之間的距離為步長(zhǎng)h，節(jié)點(diǎn)A,B,C,D,E上的電位分別用φ0,φ1,φ2,φ3和φ4表示。

圖1 區(qū)域磁場(chǎng)網(wǎng)格劃分

設(shè)函數(shù)φ在x0處可微，則沿x方向在x0處的泰勒公式展開為：

由(3)+(4)得：

同理，在y方向上有：

將式(5)、(6)代入式(1)的邊界條件中，得到泊松方程的五點(diǎn)差分格式：

式(8)表明，任意點(diǎn)的電位等于它周圍四個(gè)點(diǎn)電位的平均值。當(dāng)求解區(qū)域很大時(shí)，劃分的網(wǎng)格點(diǎn)很多，那么求解的方程組中，未知數(shù)也將很多。此時(shí)采用迭代法較為簡(jiǎn)便。下面將就常用的幾種迭代法作簡(jiǎn)要分析。

1.簡(jiǎn)單迭代法

在圖1中，將包含界別在內(nèi)的節(jié)點(diǎn)均以下標(biāo)(i,j)表示，用上標(biāo)n表示某點(diǎn)電位的第n次迭代值。由式(8)得出點(diǎn)(i,j)的第n+1次電位的計(jì)算公式為：

在以式(9)進(jìn)行迭代計(jì)算時(shí)，迭代順序可按先行后列，或先列后行進(jìn)行，逐級(jí)求出近似值。迭代過程如遇到邊界節(jié)點(diǎn)時(shí)，代入邊界值或邊界差分格式，直到連續(xù)兩次迭代求得的電位差值在允許誤差范圍內(nèi)，結(jié)束迭代。

2.高斯-賽德迭代法

為了節(jié)約計(jì)算時(shí)間，對(duì)簡(jiǎn)單迭代法進(jìn)行適當(dāng)改進(jìn)，即每計(jì)算出一個(gè)節(jié)點(diǎn)的高一次近似值，就立即用它參與其他節(jié)點(diǎn)的差分方程的計(jì)算，它的表達(dá)式可以寫為：

由于高一級(jí)迭代值的提前使用，高斯-賽德迭代法比簡(jiǎn)單迭代法快一倍左右，數(shù)據(jù)容量也小。

3.超松弛迭代法

為了進(jìn)一步加快收斂速度，采用超松弛迭代法。這里引入收斂因子α，將某節(jié)點(diǎn)的新舊電位值之差乘以收斂因子，在加到該節(jié)點(diǎn)的舊電位值上，以之作為該節(jié)點(diǎn)的新電位值。表達(dá)式如下所示：

研究可知，迭代收斂的速度與收斂因子有明顯關(guān)系，如下表所示：

表1 收斂因子與迭代次數(shù)關(guān)系表

由于收斂因子的選取一般依經(jīng)驗(yàn)而定，對(duì)于正方形磁場(chǎng)場(chǎng)域，可由下面的經(jīng)驗(yàn)公式計(jì)算最佳收斂因子：

除此之外，迭代收斂的速度還與電位初始值的給定、網(wǎng)格剖分精細(xì)和工程精度誤差有關(guān)。

三、實(shí)例分析

下面用實(shí)例闡述通過MATLAB使用超松弛迭代法求解場(chǎng)域電位的方法。例：在一個(gè)長(zhǎng)直接地金屬槽中，其橫截面如圖2所示，其側(cè)壁與底面電位均為零，底蓋電位的相對(duì)值為10，試求槽中間區(qū)段的電位分布。

圖2 長(zhǎng)直接地金屬槽

對(duì)于金屬槽中間區(qū)段的電場(chǎng)分析，可理想為二維平行平面場(chǎng)的問題，選定直角坐標(biāo)系，槽內(nèi)電位函數(shù)滿足拉普拉斯方程，構(gòu)成如下所示的第一類邊值問題：

為了全面了解超松弛迭代法的應(yīng)用，本例先在粗略的網(wǎng)格剖分下，計(jì)算槽內(nèi)九點(diǎn)區(qū)域電位的近似解，再推算區(qū)域內(nèi)任意點(diǎn)的電位情況。解題過程如下：

（1）場(chǎng)域離散化。用正方形網(wǎng)格對(duì)場(chǎng)域進(jìn)行粗略剖分，使得步距，方向的等分?jǐn)?shù)均為，使場(chǎng)域內(nèi)存在九個(gè)點(diǎn)。

（4）給定初值：進(jìn)取零值為初始值。

（5）給定迭代解收斂的指標(biāo)。規(guī)定，當(dāng)各網(wǎng)格內(nèi)點(diǎn)相鄰二次迭代的近似值的絕對(duì)誤差的絕對(duì)值均小于10-5時(shí)，迭代終止。

經(jīng)計(jì)算可知，迭代次數(shù)為13，用MATLAB求解出九點(diǎn)的電位值和電位分布圖像如下所示：

圖3 用MATLAB計(jì)算得到的九點(diǎn)電位值

圖4 用MATLAB繪制得到的電位分布圖

下面研究如何求出任意點(diǎn)的電位。該問題的關(guān)鍵在于如何設(shè)置合適的程序終止條件。為此提出如下思路：引入一個(gè)判斷矩陣P和兩個(gè)計(jì)算矩陣Va、Vb，其中Vb作為初始矩陣，Va作為迭代矩陣，判斷矩陣則作為循環(huán)執(zhí)行與否的條件，而迭代矩陣與初始矩陣之差作為是否滿足誤差限的判斷條件。當(dāng)某一次迭代使得迭代矩陣與初始矩陣之差的絕對(duì)值小于誤差限，則令判斷矩陣P=0，以此跳出循環(huán)。給出MATLAB代碼實(shí)現(xiàn)如下：

N=100;

Va=zeros(N);

Vb=zeros(N);

for i=2:N-1

Va(1,i)=10;

end

e=1e-5;

p=ones(N);

p(:,N)=0;

p(N,:)=0;

p(1,:)=0;

p(:,1)=0;

while sum(abs(p(:)))～=0

Vb=Va;

for i=N-1:-1:2

for j=2:N-1

Va(i,j)=Vb(i,j)+1.17/4*(Vb(i-1,j)+Va(i+1,j)+Vb(i,j+1)+Va(i,j-1)-4*Vb(i,j));

if abs(Va(i,j)-Vb(i,j))

p(i,j)=0;

end

end

end

end

contour(Va);

grid on

通過更改程序中的N值，可以求出場(chǎng)域點(diǎn)陣中任意點(diǎn)的電位值。假令N=100，繪制出電位分布曲線如圖5所示：

圖5 用MATLAB繪制得到100×100點(diǎn)陣的電位分布圖

[1]謝處方,饒克謹(jǐn)．電磁場(chǎng)與電磁波(第四版)[M]．北京:高等教育出版社,2006．

[2]倪光正．工程電磁場(chǎng)原理(第二版)[M]．北京:高等教育出版社,2009．

[3]葛哲學(xué)．精通MATLAB[M]．北京:電子工業(yè)出版社,2008．

[4]王飛,裴永祥．有限差分方法的MATLAB編程[J]．新疆師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,22(4):22-27．