改進Douglas分裂方法求解反應擴散方程的全離散誤差分析①

李平，陳浩

重慶師范大學 數學科學學院， 重慶 401331

反應擴散方程的應用非常廣泛， 比如斑圖動力學[1]、圖靈結構[2]、非線性波[3]、孤子或螺旋波[4]、時空混沌[5]. 然而， 由于這類系統存在剛性項和非線性項的耦合， 所以想要有效和精確地模擬這些系統是十分困難的. 到目前為止， 關于反應擴散方程數值解法的研究已經有了很多[6-18]， 例如： 隱顯式方法[6]、隱顯式Runge-Kutta方法[7]、線性隱式Runge-Kutta方法[8]、指數時間差分法[9]、混合Runge-Kutta方法[10]、指數積分法[11]等.

本文主要利用改進Douglas分裂方法[12]求解反應擴散方程， 該方法具有良好的穩定性且計算速度很快. 但在該文獻中， 作者只對線性反應擴散方程進行了時間半離散誤差分析得到了二階收斂的結論. 對于非線性反應擴散方程的全離散誤差分析在文獻中尚未提到. 本文的主要目的是對非線性反應擴散方程進行全離散誤差分析.

本文的結構安排如下： 第一節主要介紹二維非線性反應擴散方程空間、時間離散方法； 第二節給出了所構造的數值格式的全離散誤差分析； 第三節給出了幾個空間二維及三維的數值算例； 最后在第四節給出簡短的總結.

1 離散方法

1.1 空間離散

我們考慮配備齊次Dirichlet邊界條件或齊次Neumann邊界條件下的二維非線性反應擴散方程

(1)

從而得到半離散系統

(2)

其中：uj，k(t)為u(xj，yk，t)的近似解，fj，k(t，uj，k(t))： =f(xj，yk，t，uj，k(t)). 結合齊次Dirichlet邊界條件， 半離散化系統(2)可改寫成矩陣形式

U′=AU+F(t，U)

(3)

U=[u1， 1(t)， …，uN， 1(t)，u1， 2(t)， …，uN， 2(t)， …，u1， N(t)， …，uN， N(t)]T

F(t，U)=[f1， 1(t，u1， 1(t))， …，fN， 1(t，uN， 1(t))， …，f1， N(t，u1， N(t))， …，fN， N(t，uN， N(t))]T

若為齊次Neumann邊界條件時，

1.2 時間離散

在進行時間離散時， 利用改進Douglas分裂方法[12]求解半離散系統(3)， 將其改寫為

U′=A1U+A2U+F(t，U)

(4)

2 收斂性分析

由半離散系統(3)進行分析得方程(1)的精確解格式

(5)

(6)

其中， eτA是矩陣指數. 令s=tn+v， 利用中點公式， 則等式(6)可化為

(7)

由指數函數有理逼近得

所以可得到精確解表達為

(8)

消去式(4)中的中間變量得到數值解表達式為

(9)

假設1假設F(t，U)滿足Lipschitz條件， 即存在L>0， 對?U1，U2∈RN2， ?t∈[t0，T]， 有

‖F(t，U1)-F(t，U2)‖≤L‖U1-U2‖

引理1對?τ>0有

證因為A1對稱負定， 則A1的特征值μk<0，k=1，…，N2. 于是，

同理可證A2的情況.

定理1若假設1成立， 則Modified Douglas Splitting方法是二階收斂的， 即

其中c為正常數.

證將精確解表達式(8)與數值解表達式(9)作差后取范數， 由引理1和Lipschitz條件得

(10)

其中c2為正常數.由

得到

因為

所以

則有

(11)

其中：ω=‖A1+A2‖，c3為正常數.將式(11)代入式(10)可得

en+1≤(1+βτ)en+c1τh2+c3τ3+c1τ2h2+c2τ3

依此類推得

3 數值實驗

下面給出反應擴散方程的數值算例， 主要對收斂性分析進行數值驗證， 實驗是通過Matlab實現的.

算例1考慮二維反應擴散方程 (1)， 其中，Ω=[0， 1]2，t∈[0， 1]. 邊界條件為齊次Dirichlet邊界條件， 以及初始條件

u(x，y， 0)=sin(2πx)sin(2πy)

其中， 擴散項為

f(x，y，t，u)=u-u3+e-t(8π2-2)sin(2πx)sin(2πy)+e-3tsin3(2πx)sin3(2πy)

表1 二維反應擴散方程計算結果

通過表1的結果， 我們可以得出誤差關于h，τ是二階精度的， 與收斂性分析是一致的.

算例2考慮三維反應擴散方程

其中：Ω=[0， 1]3，t∈[0， 1]. 邊界條件為齊次Dirichlet邊界條件， 以及初始條件

u(x，y，z， 0)=sin(2πx)sin(2πy)sin(2πz)

其中， 擴散項為

f(x，y，z，t，u)=u-u3+e-t(12π2-2)sin(2πx)sin(2πy)sin(2πz)+e-3tsin3(2πx)sin3(2πy)sin3(2πz)

表2 三維反應擴散方程計算結果

通過上表的結果， 我們可以得到三維反應擴散方程可以利用同樣的方法進行求解， 得到的誤差關于h，τ也是二階精度的.

算例3考慮二維Schnackenberg方程組

其中：Ω=[0，L]2，t∈[0，T]. 邊界條件為齊次Neumann邊界條件， 以及初始條件

該方程組沒有精確解， 利用該方法進行求解.取N=101，M=8 000時， 其中，a=0.130 5，b=0.769 5，γ=100，Ku=0.05，Kv=1， 當L=1，T=0.5， 1， 2， 3時， 得到數值解u(x，y，t)的圖像如圖1所示.

圖1 數值解

4 結論

參考文獻[12]中針對反應擴散方程提出了一類二階改進Douglas分裂方法， 該方法具備良好的穩定性且計算速度快的特點， 但作者僅僅給出了該方法針對線性問題的半離散誤差分析. 而本文主要采用空間二階中心差分方法對空間方向進行離散， 利用改進Douglas方法對時間方向進行離散得到相應的Modified Douglas Splitting全離散格式. 對該格式進行收斂性分析， 證明該全離散格式關于空間和時間步長是二階收斂的結論. 最后借助相關數值實驗算例進行收斂性驗證.