非局部守恒Allen-Cahn方程的高效算子分裂格式

崔晨， 吳哲， 翟術(shù)英

(華僑大學(xué) 數(shù)學(xué)科學(xué)學(xué)院， 福建 泉州 362021)

1 預(yù)備知識

近年來，帶有非局部算子的Allen-Cahn方程等相場模型引起了廣泛的關(guān)注，且在數(shù)學(xué)、物理、力學(xué)、圖像處理和材料科學(xué)等領(lǐng)域得到了成功的應(yīng)用. 例如，相變[1]、近場動力學(xué)理論[2]、圖像修復(fù)[3]和非局部熱傳導(dǎo)[4]. 由于經(jīng)典的Allen-Cahn方程不保持其初始質(zhì)量，Rubinstein等[5]提出了一個守恒型的Allen-Cahn方程，這個問題與相變的擴(kuò)散界面模型有密切的關(guān)系；隨后他們的模型在文獻(xiàn)[6]中得到了廣泛的分析和數(shù)值研究.

文中，考慮如下帶有周期邊界條件的非局部Allen-Cahn方程，即

(1)

(2)

式(2)中：卷積核J滿足ⅰ)J(x)≥0，?x∈Ω；ⅱ)J(x-y)=J(y-x)；ⅲ) J是Ω周期的.

非局部Allen-Cahn方程(1)可以看作能量泛函的L2梯度流[8]，即

(3)

對能量泛函E(u)關(guān)于時間t求導(dǎo)，可得

(4)

由上述推導(dǎo)可知，能量泛函E(u)是非增的.

與經(jīng)典Allen-Cahn方程[9-12]的大量研究相比， 非局部Allen-Cahn方程的數(shù)值結(jié)果較少.額外的非局部項給開發(fā)快速有效的算法帶來了巨大的數(shù)值挑戰(zhàn)，因此許多學(xué)者致力于這一方面的研究.Du等[13]研究了非局部Allen-Cahn問題的傅里葉譜逼近；Zhai等[14-15]結(jié)合譜方法和解析法，研究了分?jǐn)?shù)階非局部Allen-Cahn模型的快速顯式算子分裂方法；Weng等[16]對此方法進(jìn)行進(jìn)一步推廣，并給出了誤差分析；Liu等[17]研究了幾類非局部和分?jǐn)?shù)階模型的數(shù)值分析及快速算法；Guan等[18]研究了非局部Allen-Cahn和Cahn-Hilliard方程的二階凸分裂格式，且利用非線性多重網(wǎng)格方法求解離散產(chǎn)生的非線性方程；Tian等[19]對非局部擴(kuò)散方程的不同逼近形式做了比較并給出理論分析.

鑒于此，本文給出求解守恒型非局部Allen-Cahn方程的數(shù)值格式及質(zhì)量守恒定理，并通過兩個數(shù)值算例驗證算法的有效性.

2 數(shù)值格式

2.1 非局部算子的離散格式

首先，從文獻(xiàn)[20]關(guān)于非局部算子Lδ的分析，可將算子Lδ轉(zhuǎn)化為卷積形式.對于任意函數(shù)u(x，t)，有

Lδu(x，t)=(J*1)u(x，t)-J*u(x，t).

(5)

于是，可將Lδu在點(diǎn)(xi，tm)，(0≤i≤N，0≤m≤M)處離散為

(6)

(7)

(8)

因此，針對任意的N+1階向量u=(u0,u1,…,uN-1,uN)T， 有

Lδu=-h×

(9)

式(9)中：*i的值為所在行其余值和的相反數(shù).

2.2 算子分裂法求解非局部Allen-Cahn方程

u(t+τ)=SA(t)SB(t)SC(t)u(t)+O(τ).

(10)

式(10)中：非線性方程解析求解，非局部方程利用C-N格式求解，拉格朗日乘子方程利用數(shù)值積分求解，從而得到簡單高效的數(shù)值格式.

非線性方程SA可通過解析式得到， 即

(11)

式(11)中：um表示第m層數(shù)值解.

(12)

對上述等式進(jìn)行整理，并寫成矩陣形式有

(13)

式(13)中：I表示單位矩陣.

拉格朗日乘子方程SC可通過數(shù)值積分離散為

(14)

(15)

(16)

將式(16)代入式(14)，可得

(17)

結(jié)合等式(11)，(13)，(17)， 可獲得算子分裂格式為

(18)

證明：方程(18)的第3個式子與e=(1,1,…，1)T作內(nèi)積，并對i求和，可得

(19)

3 數(shù)值算例

通過數(shù)值算例來驗證理論分析的正確性，包括數(shù)值格式的收斂性、能量穩(wěn)定性、質(zhì)量守恒等. 所有的算例均選擇空間區(qū)域Ω=[-1，1]，高斯核J由下列等式給出[20]，即

(20)

3.1 算例1

取初值u0(x)=0.1cos(2πx)，x∈[-1，1].為驗證空間收斂階，將時間剖分固定為M=3000.計算當(dāng)ε=0.1，δ=0.5，T=1時，不同N值的空間收斂階， 如表1所示.由表1可知： 隨著網(wǎng)格剖分變細(xì)，E2及E∞變得越來越小， 且收斂精度逐漸接近預(yù)期的二階收斂.

表1 空間收斂階(ε=0.1，δ=0.5，T=1，M=3 000)Tab.1 Spatial convergence order of different schemes (ε=0.1，δ=0.5，T=1，M=3 000)

將空間剖分固定為N=3000， 計算當(dāng)ε=0.1，δ=0.50，T=1時，不同M值的時間收斂階， 如表2所示.由表2可知：隨著網(wǎng)格剖分變細(xì)，E2與E∞均變得越來越小， 且此方法在時間上達(dá)到一階精度.

表2 時間收斂階(ε=0.1，δ=0.50，T=1，N=3 000)Tab.2 Spatial convergence order of different schemes (ε=0.1，δ=0.5，T=1，N=3 000)

3.2 算例2

將能量函數(shù)E(u)進(jìn)行離散，可得

(21)

取初值u0(x)=0.2sin(πx)cos(πx)，x∈[-1，1]， 其相應(yīng)的參數(shù)為ε=0.1，T=200，N=500，M=1 000.分別取δ=0.10，0.20，0.21，0.30，得到不同δ值的數(shù)值解、能量變化及質(zhì)量圖像，分別如圖1～2所示.

(a) 數(shù)值解(δ=0.10) (b) 能量變化(δ=0.10)

(c) 數(shù)值解(δ=0.20) (d) 能量變化(δ=0.20)

(e) 數(shù)值解(δ=0.21) (f) 能量變化(δ=0.21)

(g) 數(shù)值解(δ=0.30) (h) 能量變化(δ=0.30)圖1 算例2不同δ值的數(shù)值解和能量變化圖Fig.1 Numerical solution and energy change diagram of different δ values of example 2

(a) δ=0.10 (b) δ=0.20 (c) δ=0.30圖2 算例2不同δ值的質(zhì)量圖Fig.2 Quality map of different δ values of example 2

由圖1可知:能量函數(shù)E(u)隨著時間t的增加而減小，滿足能量遞減規(guī)律，且能量耗散與δ的大小有關(guān)，δ越大， 方程達(dá)到穩(wěn)態(tài)所需時間越長.從圖2可知：對于不同δ值，所給出格式均滿足質(zhì)量守恒的性質(zhì).

4 結(jié)束語

文中研究了守恒型非局部Allen-Cahn方程. 首先利用算子分裂方法將原方程分解為3個子問題，并構(gòu)造相應(yīng)數(shù)值求解格式.理論分析表明，格式滿足質(zhì)量守恒. 最后，通過數(shù)值算例驗證了所給格式的有效性.