李小珍

(安徽國防科技職業(yè)學(xué)院)

0 引言

拋物型方程屬于一種偏微分方程，這一類方程往往與人類的生產(chǎn)生活密切相關(guān)，可以通過求解偏微分方程來解釋一些自然現(xiàn)象和物理意義.求解偏微分方程的數(shù)值解，目前來說，主要歸結(jié)為直接法和迭代法這兩種常見的方法[1].其中，直接法是通過有限差分方式把微分方程中的微分形式用差商來代替，以此將一個拋物型偏微分方程轉(zhuǎn)化為便于研究的代數(shù)方程或方程組.但這種方法也有著諸多不便.相較于迭代法來說，它占用內(nèi)存大，運算速度慢[2].通過討論兩種較常見的迭代法，即雅可比迭代法、高斯-賽德爾迭代法剖析網(wǎng)格算法的優(yōu)勢.這兩種迭代法均可收斂，但通過計算比較得出，在相同精度下，高斯-賽德爾迭代法所需迭代的次數(shù)更少，可節(jié)省更多時間[3-4].然而，高斯-賽德爾迭代法的計算量依然存在計算量較大和收斂速度慢的問題[5].因此，利用兩重網(wǎng)格算法進(jìn)行拋物方程迭代求解[6].兩層網(wǎng)格方法的主要思想就是在細(xì)網(wǎng)格上迭代，再在粗網(wǎng)格上校正，如果沒有達(dá)到所需誤差精度便一直重復(fù)上述兩步驟.這種方法運算速度快，占用內(nèi)存少，一般來說，僅需很少的迭代次數(shù)便可達(dá)到所需的精度.

1 拋物方程及算法概述

1.1 拋物型偏微分方程

拋物型偏微分方程是偏微分方程的一種.通過判別式的符號來對偏微分方程進(jìn)行分類，在微分方程[1]

(1)

在公式(1)中，x表示在一維平面中的位置，t表示時間.當(dāng)f(x,t)=0的情況，即拋物型偏微分方程.一維拋物型微分方程中的典型例子是一維的熱傳導(dǎo)方程(在此處與擴散方程表達(dá)式相同，因為熱擴散的方式是一樣的).基本形式可表示為[5]：

(2)

有邊值條件為u(0,t)=u(1,t)=0，并有初值條件u(x,0)=u0(x).一維的熱傳導(dǎo)方程可模擬桿中的熱傳導(dǎo)，在這種情況時，u(x,t)表示熱流的密度，它與溫度成正比；K表示熱傳導(dǎo)系數(shù)，它取決于材料和它自身導(dǎo)熱的程度.為了后文研究和計算方便，(2)式中取熱傳導(dǎo)系數(shù)K=1，即化為

(3)

1.2 迭代法的基本概念

(1)迭代法的基本概念

對于任意一個給定的線性方程組Ax=b，其中矩陣A=(aij)n×n，b=(b1,…,bn)T.假設(shè)它有唯一解，可以通過迭代產(chǎn)生序列,并構(gòu)造以下的向量序列

x(k+1)=Mx(k)+g

(4)

其中，k表示迭代次數(shù)(k=0,1,2,…,n)，M為n階矩陣，對于給定的線性方程組，用公式(4)逐步代入求近似解的方法稱為迭代法[5].

(2)雅可比迭代法[3]

雅可比迭代法是以普魯士的著名數(shù)學(xué)家雅可比(Carl Gustav Jacobi)來命名的.假設(shè)一個線性方程組為Ax=b，采用Jacobi的方法，則要將A分解為三部分，分別為下三角矩陣L，對角陣D，上三角矩陣U，A=L+U+D，如式(5)所示：

也就是說，b=(L+D+U)x，移項得:Dx=(L+U)x+b，解得：x=(D)-1(L+U)+(D)-1b.對于公式(4)，經(jīng)過如下一系列迭代：

(6)

得到迭代公式：

(3)高斯-賽德爾迭代法

高斯-賽德爾迭代法是以德國數(shù)學(xué)家卡爾·弗里德里希·高斯和路德維?！べ惖聽杹砻?，也稱為李伯曼方法或連續(xù)位移方法[4].先假設(shè)一個線性方程組為Ax=b，同樣將系數(shù)矩陣A分解為三部分，對角陣D，下三角矩陣L，上三角矩陣U.不同的是，A=D-L-U，故高斯-賽德爾迭代法的迭代公式為

(i=1,2,…,n;k=0,1,2,…,n)

(7)

(4)兩層網(wǎng)格算法

兩層網(wǎng)格方法主要思想是通過細(xì)網(wǎng)格迭代來減少誤差中的高頻部分，粗網(wǎng)格校正來減少誤差的低頻部分[6].通過這種方法，可以滿足“收斂速度與網(wǎng)格步長h無關(guān)”這種要求.

拋物方程用有限差分法將其化作如下的矩陣形式：Au=f，然后用兩層網(wǎng)格的方法進(jìn)一步求解：

第一步：使用Jacobi迭代(也可用其他迭代如Gauss-Seidel迭代)迭代V1次(V1可以自由選取合適值，此處V1選取3)，記得結(jié)果為un.

第二步：計算余量rh,通過等式rh=fn-An·un.并做余量限制r2h=Rh*rh.

第三步：解粗網(wǎng)格方程：A2h·E2h=r2h得到E2h,其中A2h可由之前定義的差值矩陣和約束矩陣得到，為：A2h=Rh·Ah·Ih；

第四步：做粗網(wǎng)格校正，并將校正量加到uh上，即此時的近似解u=uh+Eh.校正量Eh可通過粗網(wǎng)格計算Eh=Ih·E2h得到；

第五步：計算此時的誤差精度是否達(dá)到要求，如達(dá)到，則已完成所有步驟，上一步驟中的u即為所需要的解；如沒有達(dá)到，則返回第一步中的細(xì)網(wǎng)格迭代，迭代方程中的un用最新得到的u來代替.如此循環(huán)幾次，直到達(dá)到接近真實值為止.

2 結(jié)果與討論

2.1 不同算法的數(shù)值解和精確解比較

利用Matlab軟件使用Jacobi迭代法、Gauss-Seidel迭代法和兩層網(wǎng)格算法對一維熱傳導(dǎo)拋物方程進(jìn)行求解.當(dāng)n=10,k=1時的數(shù)值解和精確解結(jié)果如圖1所示.

圖1 一維熱傳導(dǎo)拋物方程不同算法的數(shù)值解與精確解

由圖1可知，Jacobi迭代法、Gauss-Seidel迭代法和兩層網(wǎng)格法的數(shù)值解與精確解走勢圖變化趨勢非常接近，未出現(xiàn)多峰等現(xiàn)象，說明三者計算結(jié)果的準(zhǔn)確度相差不大.

2.2 Jacobi和Gauss-Seidel迭代法中n與k的關(guān)系

(1)當(dāng)k=1時的迭代次數(shù)

為剖析Jacobi迭代法和Gauss-Seidel迭代法中n值和k值對迭代次數(shù)的影響.先固定k=1，取不同的n值，如10,20,40,80,160，320，取誤差值精度為0.01，比較此時Jacobi迭代的次數(shù).由Matlab軟件計算，結(jié)果見表1.

表1 k=1時不同算法的迭代次數(shù)結(jié)果

(2)當(dāng)n=160時的迭代次數(shù)

然后固定n=160，取不同的k值1,2,3,4,5,6，取誤差值為0.01，比較此時Jacobi迭代的次數(shù)，結(jié)果見表2.

表2 n=160時不同算法的迭代次數(shù)結(jié)果

由表1和表2可以看出，當(dāng)k值固定時，n值越大，則達(dá)到誤差精度所需迭代的次數(shù)越多；當(dāng)n值固定，k值越大，達(dá)到誤差精度所需迭代的次數(shù)反而越少.但通過比較兩種迭代法所需迭代次數(shù)發(fā)現(xiàn)，Gauss-Seidel迭代所需次數(shù)遠(yuǎn)遠(yuǎn)少于Jacobi迭代的次數(shù).

2.3 兩層網(wǎng)格算法運算結(jié)果

兩層網(wǎng)格方法與Jacobi迭代及Gauss-Seidel的迭代方法不同，它的收斂速度不會隨著h的減小而變慢，也就是說，它的迭代次數(shù)與h的取值無關(guān).為驗證這一點，對h分別取值10,20,40,80,160,320，再計算其迭代次數(shù)，結(jié)果見表3.

表3 兩層網(wǎng)格算法迭代次數(shù)運算結(jié)果

從表3可以看出，兩層網(wǎng)格方法所需的迭代次數(shù)的確與n值無關(guān)，這就證明了它的收斂速度不會隨著h的減小而變慢.同時，這種方法所需的迭代次數(shù)很少，當(dāng)n值達(dá)到20后只需要迭代一次，由此可看出兩層網(wǎng)格方法在運算速度和占用運行空間方面的優(yōu)越性和高效性遠(yuǎn)遠(yuǎn)高于經(jīng)典的迭代方法.

3 結(jié)論

偏微分方程及拋物型的偏微分方程的研究，對物理科學(xué)和其他自然科學(xué)具有重要意義.通過比較三種迭代法(Gauss-Seidel迭代法、Jacobi迭代法和兩層網(wǎng)格法)求解一維拋物方程的精度與影響因素，得到以下結(jié)論：

(1)在相同誤差精度條件下，Gauss-Seidel迭代法比Jacobi迭代法的迭代次數(shù)少，收斂速度更快.

(2)取不同的n值和k值時所造成的對迭代次數(shù)的影響.發(fā)現(xiàn)了n值越大，即h越小時，迭代次數(shù)越多；k值越大，所需的迭代次數(shù)反而越少.

(3)兩層網(wǎng)格算法可以解決h對方程收斂速度的影響問題，通過細(xì)網(wǎng)格上迭代，再在粗網(wǎng)格上校正，僅需很少的迭代次數(shù)便可達(dá)到所需的精度.