基于MQ-RBF-FD求解對流擴(kuò)散方程的二階算法

姚林，唐泉

(新疆師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，新疆 烏魯木齊 830017)

對流擴(kuò)散方程能夠模擬許多物理、化學(xué)、經(jīng)濟(jì)和金融的現(xiàn)象從而進(jìn)行預(yù)判，對其數(shù)值解的研究由來已久，已有的數(shù)值方法包括有限差分(FD)[1]、有限元(FE)[2]和邊界元方法[3]等。一般采用局部插值格式[4]進(jìn)行推導(dǎo)，需要一個(gè)網(wǎng)格來支持操作。對流擴(kuò)散方程的特點(diǎn)是包括一階導(dǎo)數(shù)的對流項(xiàng)和二階導(dǎo)數(shù)的擴(kuò)散項(xiàng)，對于網(wǎng)格的選擇要求很嚴(yán)格，如果選擇的網(wǎng)格結(jié)構(gòu)不合適，數(shù)值逼近將會產(chǎn)生迎風(fēng)效應(yīng)和數(shù)值擴(kuò)散。徑向基函數(shù)(RBF-FD)方法[5]是一種局部有效的無網(wǎng)格方法，廣泛用于求解多尺度和偏微分方程等，解決了選擇合適網(wǎng)格的困難，具有數(shù)值精度高、構(gòu)造數(shù)值格式容易和程序操作簡單等優(yōu)點(diǎn)。因此，本文將基于RBF-FD方法的優(yōu)點(diǎn)，構(gòu)造出穩(wěn)定和有效的數(shù)值格式，避免出現(xiàn)數(shù)值振蕩和數(shù)值擴(kuò)散。

本文主要研究對流擴(kuò)散方程的數(shù)值逼近問題，考慮這類方程：

邊界條件和初始條件為

其中，Ω是二維區(qū)域，t∈(0,T]，T是終止時(shí)間，κ=[κx,κy]T是擴(kuò)散系數(shù)，ν=[νx,νy]T是對流系數(shù)，u=u(x,t)表示依賴時(shí)間t的函數(shù)，g(x,t)表示邊界函數(shù)，u0(x)表示初始條件。

本文提出一種新穎的二階離散格式求解對流擴(kuò)散方程[6]。采用二次元局部的徑向基函數(shù)(MQ-RBF-FD)方法進(jìn)行空間離散[7]，時(shí)間分裂使用交替迭代算子分裂方法[8-9]，BDF2方法用于時(shí)間離散。首先，使用維數(shù)分裂格式把二維方程轉(zhuǎn)化成沿X和Y方向的算子，梯度算子和拉普拉斯算子使用MQ-RBF-FD公式進(jìn)行離散[10]，推導(dǎo)出空間二階格式；其次，時(shí)間分裂采用交替迭代格式，可以有效減少分裂誤差，迭代格式使用BDF2方法進(jìn)行離散，目的是使時(shí)間格式上也達(dá)到二階精度；最后，進(jìn)一步測試出適當(dāng)?shù)牡綌?shù)，目的是平衡時(shí)間和空間步長，再選取適合的形狀參數(shù)c，最終使得數(shù)值解達(dá)到高階精度。數(shù)值實(shí)驗(yàn)?zāi)軌蝌?yàn)證本文提出的二階方法要大大優(yōu)于FD方法。

1 空間二階算法的推導(dǎo)過程

我們應(yīng)用五點(diǎn)RBF-FD公式離散對流擴(kuò)散方程中的梯度算子和拉普拉斯算子。采用維數(shù)分裂格式分裂二維對流擴(kuò)散方程：

其中

(1)

同理

這里使用φi(x)替換上式中的u(x)

φi(x1)=ω1φ1(x1-h)+ω2φ2(x1)+ω3φ3(x1+h)，

(2)

Δφi(x1)=γ1φ1(x1-h)+γ2φ2(x1)+γ3φ3(x1+h)，

(3)

(4)

公式(4)表示MQ徑向基函數(shù)，從式(2)和(3)中求出ωi,γi(i=0,1,2)，同理求出ωi,γi(i=0,3,4)[11-12]。Ax,Ay離散格式如下：

其中，h是空間步長，通過維數(shù)分裂，MQ-RBF-FD方法，簡單技巧處理，即可在空間格式上獲得二階格式，相比于FD方法，精度更高。通過空間離散，二維對流擴(kuò)散方程轉(zhuǎn)化成常微分方程組，下面討論時(shí)間分裂和時(shí)間離散方法。

2 時(shí)間分裂和離散的推導(dǎo)形式

時(shí)間分裂使用交替迭代格式[13]，時(shí)間離散采用BDF2方法。按照第1節(jié)處理方法，可以得到全離散格式：

為達(dá)到二階精度，我們需要找到最優(yōu)迭代步數(shù)去平衡空間步長和時(shí)間步長。下面通過兩個(gè)數(shù)值算例驗(yàn)證算法的合理性。

3 數(shù)值實(shí)驗(yàn)

3.1 數(shù)值實(shí)例1

算例1，我們考慮熱傳導(dǎo)方程：

這里u(x,y,t)屬于溫度函數(shù)，Ω=[0,1]×[0,1]，精確解：

uexact(x,y,t)=exp(-(κx+κy)π2t)sin(π(x+y))。

首先，通過數(shù)值實(shí)驗(yàn)選出最優(yōu)形狀參數(shù)c，我們選擇擴(kuò)散系數(shù)κx=κy=0.01，τ=h，T=1，迭代步數(shù)i=5。圖1說明最優(yōu)形狀參數(shù)c在1附近，我們統(tǒng)一選擇c=1。表1表示收斂階和誤差，說明MQ-RBF-FD方法優(yōu)于FD方法。其次，選擇擴(kuò)散系數(shù)κx=0.1,κy=0.001,τ=h,T=1，形狀參數(shù)c=1，迭代步數(shù)i=5。表2同樣表示收斂階和誤差，說明擴(kuò)散系數(shù)選取不同影響收斂階和誤差。

表1 RBF-FD方法和FD方法收斂階和誤差比較(κx=κy=0.01)Table 1 Convergence order and error comparison for RBF-FD and FD method(κx=κy=0.01)

圖1 不同的形狀參數(shù)c之間的誤差比較 Fig.1 Error comparison among different shape parameters

表2 RBF-FD方法和FD方法收斂階和誤差比較(κx=0.1，κy=0.001)Table 2 Convergence order and error comparison for RBF-FD and FD method(κx=0.1，κy=0.001)

3.2 數(shù)值實(shí)例2

算例2，我們考慮含有對流項(xiàng)的對流擴(kuò)散方程：

解析解：

首先，我們選取擴(kuò)散系數(shù)κx=κy=1，對流系數(shù)νx=νy=1，a=1,b=0.1,τ=0.5h,T=0.1，選取合適的迭代步數(shù)i=121，最優(yōu)形狀參數(shù)c=2，表3表示收斂階和誤差。其次，選取對流系數(shù)νx=νy=10，最優(yōu)形狀參數(shù)c=0.2，其他參數(shù)不變，表4表示收斂階和誤差。表3和表4比較說明對流系數(shù)不同，結(jié)果影響很大。最后，選取空間步長h=1/30，迭代步數(shù)i=261，其他參數(shù)與表3系數(shù)相同，得到數(shù)值解。圖2表示數(shù)值解圖和精確解圖之間的比較。

表3 RBF-FD方法和FD方法收斂階和誤差比較(υx=υy=1，c=2)Table 3 Convergence order and error comparison for RBF-FD and FD method(υx=υy=1，c=2)

表4 RBF-FD方法和FD方法收斂階和誤差比較(υx=υy=10，c=0.2)Table 4 Convergence order and error comparison for RBF-FD and FD method(υx=υy=10，c=0.2)

4 結(jié)論

本文提出一種新穎的二階算法求解對流擴(kuò)散方程，著重描述了二階算法的推導(dǎo)過程以及時(shí)間離散的推導(dǎo)形式，闡述了維數(shù)分裂和交替迭代格式在處理對流擴(kuò)散方程中的優(yōu)勢，說明了迭代步數(shù)在平衡時(shí)間和空間步長中的重要作用。數(shù)值算例印證了本文提出的二階算法的合理性和可行性，且精度遠(yuǎn)遠(yuǎn)高于FD方法。將來我們會考慮給出二階算法的理論證明。