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

李曹杰，張海湘，楊雪花

（湖南工業大學 理學院，湖南 株洲 412007）

1 研究背景

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

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

現討論如下一類橢圓型Dirichlet 邊值問題：

其邊界條件為

式（1）（2）中：Ω為矩形區域[0,L1]×[0,L2]內部；Γ為邊界。

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

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

對結點表示如下：

在結點處考慮方程（1），有

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

2 理論基礎

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

證明由帶微分余項的Taylor 公式，有

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

對g的二階導數應用Taylor 公式求解，可得

由式（10）（11）可得：

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

引理1 證畢。

引理2設

為Ω上的網格函數，記

當步長h1、h2滿足

則有

證明采用反證法，若

記

則一定存在內部結點(i0,j0)使得vi0,j0=M，且其周圍結點至少有一個的值嚴格小于M，因此當步長滿足上述關系時，有

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

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

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

記

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

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

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

定義網格函數

則有

因此有

由引理2，有

于是

引理3 證畢。

3 數值離散格式

3.1 緊致差分格式

現設u(x,y)有八階連續偏導，在式（3）兩端同時作用緊算子AB，因為AB=BA，可得

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

為使算式簡潔，記誤差項

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

省去誤差項，用數值解ui,j代替真解Ui,j，得到如下緊致差分格式：

3.2 高階Richardson 外推差分格式

若p(h)可近似表示為

用h/2 代替h，得

結合式（39）（40），有

定理1設定解問題

存在光滑解，則有

證明記

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

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

記

記

由引理3，且φ=0，有

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

結合式（53）與（54），再根據式（41），可得

定理1 證畢。

利用文獻[7]，可得一個具有六階截斷誤差的差分格式：

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

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

4 數值算例

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

表1 算例1 數值結果Table 1 Numerical results of example 1

表2 算例2 的數值結果Table 2 Numerical results of examples 2

算例1

算例2

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

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

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

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

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

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

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

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

表5 顯示，兩個算例在不同步長下的L2收斂率都約為7 階收斂，即數值方法得到的收斂率略小于理論的8 階收斂。可能的原因是：每一次步長的計算結果需要由3 種不同步長的數值解通過線性關系得出。而3 種不同步長的數值解本身就會產生舍入誤差，在高精度下三者疊加會影響到截斷誤差的收斂率。由此多次外推法有一定的局限性。

5 結語

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