一類橢圓型Dirichlet邊值問題的高精度Richardson外推法

李曹杰，張海湘，楊雪花

（湖南工業(yè)大學(xué) 理學(xué)院，湖南 株洲 412007）

1 研究背景

橢圓型方程在流體力學(xué)、彈性力學(xué)、電磁學(xué)、幾何學(xué)和變分法中都得到了廣泛應(yīng)用。隨著橢圓型方程應(yīng)用的不斷深入，對其數(shù)值解的精度要求越來越高。對這類方程，常見的處理方式是，首先對區(qū)間進(jìn)行網(wǎng)格剖分，然后對方程離散、建立差分格式，再解線性方程組[1]。五點(diǎn)差分格式是最傳統(tǒng)的差分格式，但它的精度只有二階，難以得到高精度的解，不再滿足目前的方程應(yīng)用需要。Richardson 外推法是一種通過線性組合縮小誤差的方法，因此，為了尋求高精度的解，采用外推法獲取方程的高階格式是一個(gè)不錯(cuò)的選擇。在文獻(xiàn)[1-2]中，作者基于五點(diǎn)差分格式外推，得到了橢圓型方程的四階精度格式。文獻(xiàn)[1,3-4]發(fā)展了四階緊致差分格式，可以直接求出橢圓型方程四階精度的數(shù)值解。文獻(xiàn)[5]提出了多項(xiàng)復(fù)合型黏彈性波問題的離散奇異卷積方法。文獻(xiàn)[6-8]分別在不同條件下給出了橢圓型方程六階精度的解，文獻(xiàn)[9]通過引入新的變量來構(gòu)建方程的緊致格式，文獻(xiàn)[10]利用Hopf-Cole 變換，將非線性方程變成線性方程。如此可以借助緊算子構(gòu)建高階差分格式。文獻(xiàn)[11]提出了新的時(shí)空平衡的Sinc 配點(diǎn)方法，求解四階帶奇異的積分微分方程。文獻(xiàn)[12]提出了ADI 配點(diǎn)方法求解二階帶奇異的積分微分方程。文獻(xiàn)[13]提出了計(jì)算高維帶弱奇異核發(fā)展型方程的交替方向隱式歐拉方法。文獻(xiàn)[14]構(gòu)造了二維帶弱奇異核拋物型積分微分方程的交替方向隱式有限差分格式。

本文擬基于文獻(xiàn)[1]中的四階緊致差分格式和文獻(xiàn)[7]中的六階緊致差分格式，首先，給出相關(guān)引理，搭建四階差分格式和極值原理，并證明差分格式的收斂階，在此基礎(chǔ)上建立四階緊致差分格式；然后，用Richardson 外推法，建立橢圓型方程的外推格式，分別得到其六階與八階的高階精度格式；最后，以兩個(gè)算例來驗(yàn)證高階格式的有效性。

現(xiàn)討論如下一類橢圓型Dirichlet 邊值問題：

其邊界條件為

式（1）（2）中：Ω為矩形區(qū)域[0,L1]×[0,L2]內(nèi)部；Γ為邊界。

將區(qū)間[0,L1]作m等分，記h1=L1/m；將區(qū)間[0,L2]作n等分，記h2=L2/n，則矩形區(qū)域被剖分成m×n個(gè)矩形小塊。

為表示方便，本文引入如下記號：

對結(jié)點(diǎn)表示如下：

在結(jié)點(diǎn)處考慮方程（1），有

為了引入緊致差分法，定義如下x和y方向緊算子：

2 理論基礎(chǔ)

引理1如果函數(shù)g(x)在區(qū)間[c-h,c+h]有八階連續(xù)偏導(dǎo)，則

證明由帶微分余項(xiàng)的Taylor 公式，有

將式（7）和（8）相加，整理后有

對g的二階導(dǎo)數(shù)應(yīng)用Taylor 公式求解，可得

由式（10）（11）可得：

將式（12）減去式（9），并應(yīng)用積分中值定理及介值定理，可得

引理1 證畢。

引理2設(shè)

為Ω上的網(wǎng)格函數(shù)，記

當(dāng)步長h1、h2滿足

則有

證明采用反證法，若

記

則一定存在內(nèi)部結(jié)點(diǎn)(i0,j0)使得vi0,j0=M，且其周圍結(jié)點(diǎn)至少有一個(gè)的值嚴(yán)格小于M，因此當(dāng)步長滿足上述關(guān)系時(shí)，有

這與題設(shè)矛盾，故假設(shè)不成立。引理2 證畢。

引理3緊算子A、B定義如式（4）（5）所示，有如下差分格式

這里g(xi,yj)≤0，若步長滿足如下關(guān)系：

記

為上述差分格式（21）的解，則有

證明由g(xi,yj)≤0，有

簡記g(xi,yj)為gi,j，再記

定義網(wǎng)格函數(shù)

則有

因此有

由引理2，有

于是

引理3 證畢。

3 數(shù)值離散格式

3.1 緊致差分格式

現(xiàn)設(shè)u(x,y)有八階連續(xù)偏導(dǎo)，在式（3）兩端同時(shí)作用緊算子AB，因?yàn)锳B=BA，可得

由引理1，記Ui，j=u(xi，yj)，整理得：

為使算式簡潔，記誤差項(xiàng)

將式（35）（36）代入方程（34），整理得到

省去誤差項(xiàng)，用數(shù)值解ui,j代替真解Ui,j，得到如下緊致差分格式：

3.2 高階Richardson 外推差分格式

若p(h)可近似表示為

用h/2 代替h，得

結(jié)合式（39）（40），有

定理1設(shè)定解問題

存在光滑解，則有

證明記

差分格式（38）的誤差方程組為

將式（42）（43）離散后有

記

記

由引理3，且φ=0，有

用（h1/2,h2/2）代替（h1,h2），得到

結(jié)合式（53）與（54），再根據(jù)式（41），可得

定理1 證畢。

利用文獻(xiàn)[7]，可得一個(gè)具有六階截?cái)嗾`差的差分格式：

通過定理1 類似證明，容易得到如下八階外推格式：

同樣，通過外推原理，可以對四階格式（38）進(jìn)行兩次外推，得到如下兩次外推八階格式：

4 數(shù)值算例

本節(jié)將給出數(shù)值算例，以驗(yàn)證本研究所建立方法的有效性，以及在數(shù)值求解橢圓型Dirichlet 邊值問題時(shí)的高效性。表1 和表2 中列出了算例1 與算例2 在四階格式（38）、六階格式（55）（56）、八階格式（57）（58）下的數(shù)值結(jié)果。其中，E(h1,h2)表示在步長h1、h2下網(wǎng)格結(jié)點(diǎn)的最大誤差，且Rate=log2(E(2h1,2h2)/E(h1,h2))，CPU time（s）是差分格式算該步長數(shù)值解的時(shí)間，總CPU 是程序連續(xù)運(yùn)行完各步長結(jié)果的總時(shí)間。

表1 算例1 數(shù)值結(jié)果Table 1 Numerical results of example 1

表2 算例2 的數(shù)值結(jié)果Table 2 Numerical results of examples 2

算例1

算例2

兩個(gè)算例中，一次外推六階差分格式與一次外推八階差分格式分別建立在四階、六階緊致差分格式上。從表1 和表2 中，可以看到一次外推能將收斂率提高兩階。在相同步長下這能得到更高精度的數(shù)值解，且解的收斂速度更快。兩次外推八階差分格式建立在四階緊致格式上，其L∞誤差值均低于理論預(yù)期值，因此，繪出其誤差曲面并求出部分結(jié)點(diǎn)值，以進(jìn)一步探究其原因。算例1 在步長為1/32 和1/64 時(shí)的誤差曲面如圖1所示，其不同結(jié)點(diǎn)處的誤差值與最大誤差值具體見表3。

圖1 算例1 在不同步長下的誤差曲面Fig.1 Error surface of example 1 with different step lengths

表3 算例1 不同結(jié)點(diǎn)處的誤差值與最大誤差Table 3 Error values and maximum errors of example 1 at different nodes

算例2 在步長為1/32 和1/64 時(shí)的誤差曲面如圖2所示，其不同結(jié)點(diǎn)處的誤差值與最大誤差值見表4。

圖2 算例2 在不同步長下的誤差曲面Fig.2 Error surface of example 2 with different step lengths

表4 算例2 不同結(jié)點(diǎn)處的誤差值與最大誤差Table 4 Error values and maximum errors of example 2 at different nodes

從圖1 和圖2、表3 和表4 中的數(shù)據(jù)可以看出，對比某一固定結(jié)點(diǎn)時(shí)，例如表4，在結(jié)點(diǎn)(1/2,1/4)處，其收斂率分別約為8 和9，由此可見算例收斂率符合理論預(yù)期。從圖1 和圖2 中可以看到，同一步長下，有些結(jié)點(diǎn)的誤差值相對較大，并且可以發(fā)現(xiàn)當(dāng)步長不同時(shí)，其誤差最大值在不同結(jié)點(diǎn)處取得，Rate=log2(E(2h1,2h2)/E(h1,h2))約為6，因此在前文的表1 與表2 中，二次外推格式的整體收斂率只達(dá)到了6 階。接下來計(jì)算二次外推格式的L2誤差收斂階，兩個(gè)算例的L2誤差及收斂率結(jié)果如表5所示，表中的Rate=log2(L2(2h1,2h2)/L2(h1,h2))。

表5 兩個(gè)算例的L2 誤差及收斂率Table 5 L2 error and convergence rate of two examples

表5 顯示，兩個(gè)算例在不同步長下的L2收斂率都約為7 階收斂，即數(shù)值方法得到的收斂率略小于理論的8 階收斂?？赡艿脑蚴牵好恳淮尾介L的計(jì)算結(jié)果需要由3 種不同步長的數(shù)值解通過線性關(guān)系得出。而3 種不同步長的數(shù)值解本身就會產(chǎn)生舍入誤差，在高精度下三者疊加會影響到截?cái)嗾`差的收斂率。由此多次外推法有一定的局限性。

5 結(jié)語

本文研究了以基于緊致差分格式的高精度Richardson 外推格式來求橢圓型偏微分方程的數(shù)值解，理論顯示一次外推能將精度提升兩階。算例對比了原緊致差分格式、外推高階差分分格式和同階緊致差分格式的數(shù)值結(jié)果。結(jié)果表明，外推高階差分格式雖然比原格式要花費(fèi)更長的時(shí)間進(jìn)行運(yùn)算，但同步長下得到解的精度大大提升。且外推差分格式運(yùn)算時(shí)間不一定比同階緊致差分格式更長，總CPU 時(shí)間甚至可能會更短，外推格式時(shí)間主要取決于原格式程序的運(yùn)算簡易程度。兩次外推八階差分格式理論上能達(dá)到八階，但由于結(jié)果會受到多次舍入誤差的影響，使得計(jì)算結(jié)果受到影響。L∞誤差只達(dá)到了六階收斂；L2誤差收斂率在七階左右。兩個(gè)算例的結(jié)果均顯示，外推八階格式精度符合理論預(yù)期?？偟膩碚f，Richardson 外推法對于提高差分格式的精度有顯著作用，且易于理解與使用，便于在差分格式中廣泛利用。