求解二維聲波方程的高精度Runge-Kutta 方法

陳 麗，朱興文，張朝元*

（1.大理大學(xué)工程學(xué)院，云南大理 671003；2.大理大學(xué)數(shù)學(xué)與計(jì)算機(jī)學(xué)院，云南大理 671003）

過去的地震波計(jì)算方法〔1-5〕因存在數(shù)值頻散的不良缺點(diǎn)滿足不了現(xiàn)在大尺度數(shù)值模擬的特性，楊頂輝教授團(tuán)隊(duì)〔6〕在2003 年將近似解析離散化（nearly analytic discrete，NAD）算子引入到模擬地震波中，已發(fā)展了系列空間離散精度為四階的有關(guān)NAD 算子的數(shù)值算法〔7-10〕，這些算法能有效抑制數(shù)值頻散現(xiàn)象的發(fā)生。

為了在地震波模擬時能更好地壓制數(shù)值頻散現(xiàn)象的發(fā)生，基于二維聲波方程，選擇八階精度的NAD 算子〔11〕去離散常微分方程中的高階空間偏導(dǎo)數(shù)，選擇三階Runge-Kutta 方法〔12〕去離散常微分方程中的時間導(dǎo)數(shù)，從而獲得八階NAD-RK 算法。針對該方法，分析了其理論誤差，結(jié)合八階Lax-Wendroff（LWC）格式和八階交錯網(wǎng)格（SG）格式兩種高階方法進(jìn)行數(shù)值誤差比較研究，詳細(xì)推導(dǎo)其穩(wěn)定性條件。

1 算法推導(dǎo)

設(shè)二維聲波方程為：

此處，d 表示位移，c 表示二維聲波速度。

再令U＝（D，P）T，則方程（2）式變?yōu)橄率剑?/p>

式（3）轉(zhuǎn)化為標(biāo)準(zhǔn)的常微分方程模式，這樣就可用常微分方程模式的求解程序來轉(zhuǎn)化為對二維聲波方程式（1）的對應(yīng)式（3）的求解。過程如下：

第一，由于位移d 和粒子速度p 的二階和三階偏導(dǎo)數(shù)存在高階偏微分算子矩陣L 中，先利用高階偏導(dǎo)數(shù)的計(jì)算公式〔11〕對式（3）右邊的L 進(jìn)行空間高階偏導(dǎo)數(shù)離散化。

第二，在離散化高階空間偏導(dǎo)數(shù)后，式（3）成為一個關(guān)于時間t 的半離散化常微分方程模式，選用三階Runge-Kutta 方法〔12〕進(jìn)行求解。三階Runge-Kutta 方法的計(jì)算公式為：

其中，Δt 為時間步長，U（1）、U（2）為中間變量，Un＝U（nΔt）。

消去式（4）中的變量U（1）、U（2），可以加快計(jì)算速度和減少存儲量，得計(jì)算公式：

把U＝（D，P）T代入式（5）中，得到如下分量的計(jì)算公式：

這里A2＝A·A。

上面計(jì)算公式（5）或式（6a）和式（6b）稱為求解二維聲波方程式（1）的八階NAD-RK 算法。

2 理論誤差

對于求解二維聲波方程的八階NAD-RK 算法，由泰勒級數(shù)展開公式知式（5）右邊的二階和三階空間偏導(dǎo)數(shù)計(jì)算公式的截?cái)嗾`差為O（Δx8+Δz8）〔11〕，即八階NAD-RK 算法在空間上是一種八階精度格式。采用三階Runge-Kutta 方法求解常微分方程（5）式時，時間導(dǎo)數(shù)誤差為O（Δt3）〔12〕。因此，求解二維聲波方程的八階NAD-RK 算法的理論誤差為O（Δt3+Δx8+Δz8）。

3 數(shù)值誤差

下面研究八階NAD-RK 算法的數(shù)值誤差。考慮以下二維初值問題：

其中，頻率f0＝25 Hz，t=0 s 時刻聲波的入射角度θ0=π/4。容易求得該初值問題的解析解為：

定義數(shù)值解的相對誤差〔11〕為：

其中d（tn，xi，zj）為精確解，為數(shù)值解。

在進(jìn)行數(shù)值計(jì)算時，計(jì)算參數(shù)選擇為：計(jì)算域0＜x，z≤10 km、頻率f0=20 Hz、波速c=4 km/s、網(wǎng)格尺寸Δx=Δz=50 m、時間步長Δt=0.001 s。圖1 給出了式（9）中相對誤差Er（%）的計(jì)算結(jié)果。圖1 中的3條曲線分別對應(yīng)八階NAD-RK 算法、八階LWC 格式和八階SG 格式。從圖1 可以看出，在3 種數(shù)值方法中，八階NAD-RK 算法的數(shù)值誤差最小。

圖1 八階NAD-RK 算法、八階LWC 格式和八階SG格式在求解初值問題（7）時在空間步長和時間步長分別為Δx=Δz=50 m 和Δt=0.001 s 情況下的相對誤差

4 穩(wěn)定性條件

為了保持?jǐn)?shù)值迭代的穩(wěn)定性，需要考慮如何選擇合適的時間和空間步長，即Δt，Δx 和Δz。在條件h=Δx=Δz 下，所定義的庫朗數(shù)α=cΔt/h 在數(shù)學(xué)上給出了速度c 與時間和空間兩個網(wǎng)格步長h 之間的關(guān)系。需要確定α 的范圍，使提出的方法穩(wěn)定。在本部分，根據(jù)傅里葉分析和Yang 等〔13〕提出的分析過程，推導(dǎo)二維聲波情況下八階NAD-RK 算法的穩(wěn)定性條件。

設(shè)

然后將諧波解：

帶入式（5），結(jié)合式（6a）和（6b）兩個高階偏導(dǎo)數(shù)的離散公式，可得到如下方程：

其中H 為增長矩陣。

由于增長矩陣H 中部分元素hij的表達(dá)式太復(fù)雜，在這里只展示了H 的前兩行元素的表達(dá)式，如下：

這里，αmax是最大庫朗數(shù)。

5 結(jié)論

為有效抑制地震波數(shù)值模擬過程中的數(shù)值頻散而獲得更高的模擬精度和更快的計(jì)算速度，本文在楊頂輝教授團(tuán)隊(duì)的研究基礎(chǔ)上，基于二維聲波方程，結(jié)合八階NAD 算子和三階Runge-Kutta 方法，推導(dǎo)出八階NAD-RK 算法。針對該方法，分析了其理論誤差，結(jié)合八階LWC 格式和八階SG 格式兩種高階方法進(jìn)行了數(shù)值誤差比較研究，詳細(xì)推導(dǎo)了其穩(wěn)定性條件。研究結(jié)果揭示，八階NAD-RK 算法的理論誤差和穩(wěn)定性條件分別為O（Δt3+Δx8+Δz8）和α≤αmax≈0.541 6；同八階LWC 格式和八階SG 格式相比，八階NAD-RK 算法具有最小的數(shù)值誤差。