汪亞楠, 蔡耀雄

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

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

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

(1)

非局部項(xiàng)G*u2的定義為

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

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

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

(2)

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

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

1 數(shù)值格式

1.1 預(yù)備知識(shí)

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

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

(3)

(4)

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

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

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

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

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

積分因子(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)兩邊同時(shí)乘以e-Lt,可得

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

(7)

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

vt=e-Ltf(eLtv).

(8)

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

u(0)=un,

(9)

(10)

un+1=u(s).

(11)

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

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

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

.

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

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

2 數(shù)值算例

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

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

2.1 算例1

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

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

表1 不同時(shí)間步長(zhǎng)下的時(shí)間誤差和時(shí)間收斂階Tab.1 Time errors and time convergence rates at different time

由表1可知：時(shí)間收斂階都達(dá)到了預(yù)期精度.

2.2 算例2

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

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

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

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

采用四步四階顯式穩(wěn)定性積分因子龍格庫(kù)塔格式可得不同時(shí)刻數(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,固定空間步長(zhǎng)h=1,時(shí)間步長(zhǎng)τ=0.1.

采用四步四階顯式穩(wěn)定性積分因子龍格庫(kù)塔格式可得不同時(shí)刻數(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時(shí)液體晶體的生長(zhǎng),初始值為

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

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

3 結(jié)束語(yǔ)

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