仇佳棟,梁琪琪,胡志濤,韓丹夫
(杭州師范大學(xué)數(shù)學(xué)學(xué)院,浙江 杭州 311121)
非線性偏微分方程在量子力學(xué)、固體物理學(xué)、流體力學(xué)和等離子體物理學(xué)等許多學(xué)科中發(fā)揮著重要作用[1].為了描述以離子聲速運動的坐標(biāo)系中一維朗繆爾(Langmuir)和離子聲波的非線性動力學(xué),耦合Schr?dinger-KdV(CSK)方程組[2-3]被提出.
考慮具有初始條件和周期性邊界條件的一維CSK方程組[4]:
iεut+puxx-qvu-s|u|2u=0, (x,t)∈×(0,T],
vt+αvxxx+(βvm+ρ|u|2)x=0, (x,t)∈×(0,T],
u(x,t)=u(x+l,t),v(x,t)=v(x+l,t), (x,t)∈×(0,T],
u(x,0)=φ(x),v(x,0)=φ(x),x∈.
在此物理模型中通常有以下幾類守恒量.
這幾類守恒量很可能與時間方向上的精度密切相關(guān).近十年來,人們提出了大量的數(shù)值方法來解決CSK方程,如有限元方法(finite element method,FEM)[5]、徑向基函數(shù)配點法[6]、分解法[7]、變分迭代[8]、Padé有理逼近的指數(shù)時間差分三層隱式格式(exponential time differencing three-layer implicit scheme,ETDT-P)[9]、同倫攝動[10]及四階緊致差分格式(fourth-order conservative compact finite difference scheme,FCS)[4]等.
本文首先介紹六階緊致差分格式并將其用于CSK方程的數(shù)值求解,隨后提出一種六階緊致差分格式的通用矩陣形式并分析該格式的離散守恒性,最后給出數(shù)值實驗的結(jié)果.
注意到空間一階導(dǎo)數(shù)在每點上的六階緊致差分[11]近似
(1)
通過計算及以上算子的定義,式(1)可以寫成如下算子形式:
(2)
類似地,空間二階導(dǎo)數(shù)的六階緊致差分[11]近似為
等價于
(3)
結(jié)合式(2)和(3),可得空間三階導(dǎo)數(shù)的六階差分近似為
(4)
(5)
(9)
記Hp(Ωh)={u|u={uj},j=0,1,…,J,uj=uj+J}為定義在Ωh上以J為周期的實值或復(fù)值網(wǎng)格函數(shù)所組成的空間.網(wǎng)格函數(shù)空間Hp(Ωh)上的離散內(nèi)積及其相應(yīng)離散L2-范數(shù)和L∞-范數(shù)定義如下:
引理1對于任意網(wǎng)格函數(shù)u,w∈Hp(Ωh),以下等式成立:
引理2對于任意網(wǎng)格函數(shù)u,w∈Hp(Ωh),以下等式成立:
由引理3,易得如下結(jié)果:
引理4[13]對于任意網(wǎng)格函數(shù)u,w∈Hp(Ωh),以下不等式成立:
關(guān)于CSK方程, Xie等[4]給出了一個四階緊致格式的守恒性,本文構(gòu)造的格式(6)、(7)是一個六階緊致格式,其守恒性定理有:
定理1格式(6)、(7)保持如下離散守恒性:
‖Un‖2=‖U0‖2,
(10)
(11)
證明令
(12)
結(jié)合式(6)和(12)得
(13)
(16)
將式(16)兩邊同時乘以iε并利用內(nèi)積的共軛對稱性可得
(17)
由引理1、式(15)和(17)得
(18)
(19)
(20)
將引理1用于式(20),即得‖Un+1‖2=‖Un‖2.將式(7)與向量1∶=(1,1,…,1)T∈Hp(Ωh)進行內(nèi)積運算,可得
〈B1B2Vn+1,1〉=〈B1B2Vn,1〉.
再利用周期性邊界條件,即得式(11)成立.
以下定義
易知格式(6)、(7)等價于以下矩陣形式:
其中循環(huán)對稱矩陣A1,A2,B1,B2是分別對應(yīng)于算子A1,A2,B1,B2的J階方陣,記為
方程(6)和(7)可寫成如下等價形式:
(22)
定理2格式(6)、(7)保持如下離散守恒性:
n=0,1,…,N.
(23)
證明將式(21)與δtUn進行內(nèi)積運算并利用循環(huán)矩陣的可交換性,可得
(24)
取式(24)的實部得
(25)
將式(25)兩邊同時乘以2τ并將上標(biāo)從0至n求和,可得
(26)
(27)
(28)
由ψ(Vn,Vn+1)的定義,易得
以及
將式(28)兩邊同時乘以2τ并將上標(biāo)從0至n求和,可得
(29)
又因為
方程(29)等價于
(30)
將式(26)和(30)分別乘以ρ和q/2并相減,即得式(23)成立.
因格式(6)和(7)是非線性的,計算較為復(fù)雜,為了保持格式原有的計算精度,本文將其轉(zhuǎn)化為如下線性迭代格式:
(34)
其中U*0=(U0*+U0)/2,U*n=(3Un-Un-1)/2,Un*=2Un-Un-1,n≥1.
可以證明線性格式(33)、(34)在時間和空間上分別是二階和六階收斂的.
本節(jié)給出一些例子來說明所提出格式的守恒性和六階精度.記守恒量I1,I2,I3和I4的離散形式為
守恒量的誤差定義為
此外,本文用L2-范數(shù)(‖u-U‖+‖v-V‖)及L∞-范數(shù)(‖u-U‖L∞+‖v-V‖L∞)來檢驗所提出格式的精度.
例1[10]考慮如下柯西問題:
iut+uxx-vu=0, (x,t)∈×(0,T],
vt+vxxx+(3v2+|u|2)x=0, (x,t)∈×(0,T],
u(x,0)=φ(x),v(x,0)=φ(x),x∈.
其解析解為u(x,t)=exp(i(x+t/4)),v(x,t)=3/4,在h=π/20、τ=0.001、空間域取[0,2π]的情況下數(shù)值求解CSK方程,表1列出了不同時刻的各個守恒量誤差,佐證了所提出格式的守恒性.
表1 不同時刻的各個守恒量誤差Tab.1 Errors of invariants at different time
取h=π/10,τ=0.1,表2列出了不同時刻不同步長下的兩種范數(shù)誤差,驗證了所提出格式在空間上六階收斂.
表2 不同時刻的收斂速率Tab.2 Convergence rates at different time
例2[5]考慮如下耦合問題:
其解析解為
其中c是任意正常數(shù),初始條件為
此外,為滿足CSK方程的物理意義即當(dāng)|x|趨于無窮時|u|和v趨于0,設(shè)置邊界條件為u(a,t)=u(b,t)=0,v(a,t)=v(b,t)=0.取參數(shù)ε=1,c=0.45,在h=0.25、τ=0.000 01、空間域取[-30,30]的情況下數(shù)值求解CSK方程.表3列出了t=0.001時,U(x,t)虛部數(shù)值解(ImU)、實部數(shù)值解(ReU)、V(x,t)的數(shù)值解及其誤差范數(shù)(‖ImEu‖, ‖ReEu‖, ‖Ev‖),以及與其他方法和精確解的比較.由表3可見,本文方法與其他方法具有相近的誤差范數(shù)(實際上這是邊界誤差造成的,因為此時空間和時間上的內(nèi)部誤差都遠小于邊界誤差),而在離邊界較遠的內(nèi)點(如0)上,本文方法的精度遠遠大于其他方法.圖1模擬了|U|、V的數(shù)值解和解析解從t=0到t=30的演變過程,可以發(fā)現(xiàn)數(shù)值解與精確解高度一致.
表3 本文方法數(shù)值解與其他方法及精確解的比較 Tab.3 Comparison of numerical solutions with exact solutions and other methods
續(xù)表3
注:t∈[0,30],[a,b]=[-70,30],h=0.5,τ=0.001.
通過構(gòu)造數(shù)個循環(huán)對稱陣得到一種通用的六階緊致差分格式并將其用于數(shù)值求解耦合Schr?dinger-KdV(CSK)方程.利用CSK方程的幾類守恒量對該數(shù)值方法的性能進行了衡量.數(shù)值實驗說明了該方法的精度和穩(wěn)定性,特別是在內(nèi)點取得了較其他方法更高的精度.因為六階緊致差分格式中構(gòu)造的數(shù)個矩陣具有與四階緊致差分格式中矩陣類似的性質(zhì),該格式還可推廣到其他用四階緊致差分格式求解的方程如Klein-Gordon-Schr?dinger方程[15]、耦合Gross-Pitaevskii方程[16]等.