非局部Swift-Hohenberg方程的積分因子龍格庫塔格式

汪亞楠, 蔡耀雄

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

Swift-Hohenberg(SH)模型[1]最初是由Swift和Hohenberg在研究Rayleigh-Bénard對流時引入的,并且已經(jīng)成為導致復雜模式形成的非線性動力系統(tǒng)的范例之一,其研究成果被廣泛應用于復雜流體和生物組織[2-3].近年來,隨著流體力學、反應擴散化學和生物學等學科的發(fā)展,非局部SH模型的研究吸引了眾多學者的關(guān)注.Roberts[4]指出,將二維局部SH模型用作特定物理系統(tǒng)中空間模式演化的可靠模型是不夠的,應使用非局部SH模型.2014年,Morgan等[5]首次提出帶有非局部非線性項的SH模型,即

ut=-(Δ+1)2u+εu-uG*u2, (x,t)∈Ω×(0,T].

(1)

非局部項G*u2的定義為

卷積核G滿足以下3個條件：1)G(x)≥0,對于任意的x∈Ω;2)G是Ω-周期性;3)G(x)=G(|x|)

是一個給定的徑向?qū)ΨQ函數(shù).

非局部SH方程(式(1))可看作能量函數(shù)E(u)的L2梯度流,有

(2)

關(guān)于局部SH方程的理論和數(shù)值研究已有大量的優(yōu)秀成果[6-7],而關(guān)于非局部SH方程的研究工作還十分有限.非局部SH方程中的非局部和非線性項給研究帶來巨大挑戰(zhàn).Firth等[8]提出非線性光學系統(tǒng)的非局部模型.Purwins等[9]提出介質(zhì)氣體放電動力學的非局部模式.Du等[10]提出非局部算子的向量演算.Zhang等[11]提出一個帶有拉格朗日乘子的守恒型非局部 SH方程.Weng等[12]利用傅里葉譜方法求解帶有非局部非線性項的SH方程.

近年來,作為求解偏微分方程的有效數(shù)值方法,積分因子龍格庫塔法受到廣泛關(guān)注.Ju等[13]提出關(guān)于半線性拋物方程的保極值原理的積分因子龍格庫塔方法.Ahmed等[14]提出關(guān)于非齊次邊界條件系統(tǒng)的高階積分因子法.Zhang等[15]提出保守Allen-Cahn方程的顯式三階保結(jié)構(gòu)格式.Nan等[16]提出求解非局部Allen-Cahn方程的高階極值原理積分因子龍格庫塔方法.基于此,本文結(jié)合積分因子龍格庫塔法和譜方法[17]對式(1)進行有效求解,并提出4種快速有效求解非局部Swift-Hohenberg方程的數(shù)值格式.

1 數(shù)值格式

1.1 預備知識

空間區(qū)域Ω=[-a,a]2上的網(wǎng)格剖分為

為了求解周期邊界條件下的非局部SH方程(式(1)),基于譜方法的相關(guān)理論,采用傅里葉級數(shù)逼近{ui,j}進行空間離散求解.

(3)

(4)

引理1[12]對于任意的網(wǎng)絡函數(shù)u∈Mh,擴散算子(Δ+1)2u及非局部卷積算子G*u的譜離散形式分別為

1.2 顯式穩(wěn)定性積分因子龍格庫塔求解非局部SH方程

相較于傳統(tǒng)的強穩(wěn)定性龍格庫塔方法,強穩(wěn)定性積分因子龍格庫塔可避免線性算子的時間步長限制,但由于非線性算子帶來的時間步長限制,需要引入一種無時間步長限制的保界格式.

顯式穩(wěn)定性積分因子龍格庫塔法(eSIFRK)[13,18-19]是一種顯性格式且保極值原理的離散方法,其時間步長的約束僅取決于非線性項,故時間步長大小的選擇與空間網(wǎng)格大小無關(guān).

采用顯式穩(wěn)定性積分因子龍格庫塔法,分別建立求解式(1)的一步一階、二步二階、三步三階、四步四階的無時間步長限制的保界格式.

積分因子(IF)法[20]的步驟如下.

首先,將式(1)改寫為

ut=Lu+f(u).

(5)

式中：Lu=-(Δ+1)2u;f(u)=εu-uG*u2.

其次,引入新變量v,定義為

v=e-Ltu.

(6)

式(6)中：e-Lt為積分因子.

然后,式(5)兩邊同時乘以e-Lt,可得

e-Ltut-Le-Ltu=e-Ltf(u).

(7)

注意到vt=e-Ltut-Le-Ltu,則式(7)可簡化為

vt=e-Ltf(eLtv).

(8)

根據(jù)文獻[13],可得式(8)的s步穩(wěn)定性積分因子龍格庫塔格式為

u(0)=un,

(9)

(10)

un+1=u(s).

(11)

給定s=1,2,3,4,可得一步一階顯式穩(wěn)定性積分因子龍格庫塔格式eSIFRK(1,1)為

un+1=eτL[un+τf(un)].

二步二階顯式穩(wěn)定性積分因子龍格庫塔格式eSIFRK(2,2)為

.

三步三階顯式穩(wěn)定性積分因子龍格庫塔格式eSIFRK(3,3)為

四步四階顯式穩(wěn)定性積分因子龍格庫塔格式eSIFRK(4,4)為

2 數(shù)值算例

數(shù)值解的L2范數(shù)誤差ErrL2(τ,h)為

數(shù)值解的最大誤差Errmax(τ,h)為

2.1 算例1

為了驗證時間收斂階,算例1的空間區(qū)域Ω=[-20,20]2,初始值為

u0(x,y)=0.01×(cos πx+cos πy+2cos 0.25πy).

表1 不同時間步長下的時間誤差和時間收斂階Tab.1 Time errors and time convergence rates at different time

由表1可知：時間收斂階都達到了預期精度.

2.2 算例2

算例2考慮帶有隨機初始擾動的非局部SH方程,其空間區(qū)域Ω=[-64,64]2,初始值為

u0(x,y)=0.07+0.001×rand(x,y).

式中：rand(x,y)表示取值在[-1,1]的隨機數(shù).

設(shè)置系數(shù)ε=0.035,δ=0.5,固定空間步長h=1,時間步長τ=1.

采用四步四階顯式穩(wěn)定性積分因子龍格庫塔格式可得不同時刻數(shù)值解圖像和能量圖(T=2 000),如圖1所示.由圖1可知：能量是遞減的.

(a) T=0 (b) T=300 (c) T=500 (d) T=700

2.3 算例3

為了研究數(shù)值解與能量的變化,算例3空間區(qū)域Ω=[-64,64]2,初始值為

設(shè)置系數(shù)ε=0.035,δ=0.5,固定空間步長h=1,時間步長τ=0.1.

采用四步四階顯式穩(wěn)定性積分因子龍格庫塔格式可得不同時刻數(shù)值解圖像及能量圖(T=1 000),如圖2所示.由圖2可知：能量是遞減的.

(a) T=0 (b) T=5 (c) T=20 (d) T=50

2.4 算例4

考慮空間區(qū)域Ω=[-150,150]2時液體晶體的生長,初始值為

采用四步四階顯式穩(wěn)定性積分因子龍格庫塔格式可得不同時刻數(shù)值解圖像及能量圖(T=1 000),如圖3所示.由圖3可知：隨著時間的推移,微晶相互碰撞并開始形成晶界;當T=1 000時,能量隨著時間逐漸減小達到穩(wěn)態(tài).

(a) T=0 (b) T=50 (c) T=200 (d) T=400

3 結(jié)束語

結(jié)合積分因子龍格庫塔法和傅里葉譜方法求解非局部SH方程,得到一種快速有效的數(shù)值格式,數(shù)值算例表明文中算法具有較好的穩(wěn)定性及較高的準確性,且滿足能量遞減性質(zhì).