鄧楊芳， 黃蓉， 翁智峰

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

1958年，Cahn等[1]在研究熱力學(xué)兩相物質(zhì)(如合金、聚合物等)之間的相互擴(kuò)散現(xiàn)象時，提出Cahn-Hilliard方程.Cahn-Hilliard方程是一類重要的四階非線性擴(kuò)散方程，用來描述固體表面上微滴的擴(kuò)散、河床的遷移、生物種群的競爭與排斥等擴(kuò)散現(xiàn)象，也廣泛應(yīng)用于工程流體力學(xué)[2-3]、材料科學(xué)[4]、圖像分析[5-6]和生命科學(xué)[7-8]等領(lǐng)域.國內(nèi)外學(xué)者對Cahn-Hilliard方程提出了許多數(shù)值解法，如有限差分[9-11]、有限元[12-13]、有限體積[14]、譜方法[15-16]等.Zhang等[17]構(gòu)造一個無條件能量穩(wěn)定的有限差分格式，并基于該格式提出自適應(yīng)時間步長策略.Cheng等[18]采用快速穩(wěn)定的顯式算子分裂方法求解Cahn-Hilliard方程.Li等[19]考慮求解二維Cahn-Hilliard方程的半隱穩(wěn)定化傅里葉譜方法，并對幾種經(jīng)典的穩(wěn)定化技術(shù)進(jìn)行比較研究，確定相應(yīng)的穩(wěn)定區(qū)域.Cheng等[20]提出一個能量穩(wěn)定的四階差分格式求解帶周期邊界條件的Cahn-Hilliard方程，并證明該格式解的唯一性和能量穩(wěn)定性.Wang 等[21]提出一個能量穩(wěn)定的線性擴(kuò)散Crank-Nicolson格式求解Cahn-Hilliard方程.Zhang等[22]構(gòu)造一類無條件能量穩(wěn)定的線性二階預(yù)估-校正時間步進(jìn)格式求解Cahn-Hilliard方程.

上述方法都是基于網(wǎng)格剖分求解微分方程問題.李淑萍等[23]利用重心插值配點法求解微分方程初邊值問題，重心插值配點法包括重心Lagrange插值和重心有理插值.對于等距插值節(jié)點，高階Lagrange插值公式構(gòu)造的近似函數(shù)值容易出現(xiàn)Runge現(xiàn)象，具有極大的數(shù)值不穩(wěn)定性.為了避免這一現(xiàn)象的出現(xiàn)，Berrut等[24]分別將Lagrange插值公式改為重心Lagrange插值和重心有理插值[25]，克服數(shù)值格式的振蕩現(xiàn)象.重心插值配點法作為一種新型的無網(wǎng)格計算方法，具有計算格式簡單、精度高、程序?qū)嵤┓奖?、?jié)點適應(yīng)性好等特點，使用Chebyshev節(jié)點有效克服了Runge現(xiàn)象.重心插值配點法被推廣到求解各類積分和微分方程，如高維Fredholm積分方程[26]、非線性拋物方程[27]、粘彈性波方程[28]、Allen-Cahn方程[29-30]、Black-Scholes方程[31]等.本文將重心插值配點法推廣到求解Cahn-Hilliard方程.

1 Cahn-Hilliard方程的重心插值配點法

1.1 Cahn-Hilliard方程

Cahn-Hilliard方程是相場模型方程之一，用來模擬二元合金中的相分離現(xiàn)象，即

φt(x，t)=Δ(F′(φ(x，t))-ε2Δφ(x，t))，x∈[a，b]，t∈[0，T].

(1)

邊界條件為

φ(a，t)=α1(t)，φ(b，t)=α2(t)，φ′(a，t)=β1(t)，φ′(b，t)=β2(t).t∈[0，T]，

(2)

初始條件為

φ(x，0)=ψ(x)，x∈[a，b].

(3)

Cahn-Hilliard方程也可看作是H-1內(nèi)積下的一個梯度流，它的自由能泛函為

(4)

式(4)中：?是梯度算子.

令μ=F′(φ)-ε2Δφ，方程(1)與μ作內(nèi)積可得

(5)

因此，Cahn-Hilliard方程滿足能量遞減規(guī)律，在采用數(shù)值方法求解Cahn-Hilliard方程時，不僅要保證數(shù)值精度，而且要保證能量的穩(wěn)定性.

1.2 重心型插值

1.2.1 重心Lagrange插值 設(shè)有插值節(jié)點tj和一組對應(yīng)的實數(shù)fj(j=0，1，…，m)，采用多項式插值，則在次數(shù)不超過m的多項式空間可以找到唯一的插值多項式p(t)，滿足p(tj)=fj，j=0，1，…，m，Lagrange插值公式為

(6)

(7)

(8)

(9)

將式(9)代入式(6)，有

(10)

由式(8)～(10)，可得重心Lagrange 插值公式，即

(11)

1.2.2 重心有理插值 設(shè)有插值節(jié)點tj(j=0，1，…，m)和一組對應(yīng)的實數(shù)fj，選擇一個整數(shù)d且滿足0≤d≤m，令pj(t)為插值d+1個點對(tj，fj)，(tj+1，fj+1)，…，(tj+d，fj+d)(j=0，1，…，m-d)的次數(shù)最多為d的多項式，則

(12)

式(12)中:λj(t)=(-1)j/((t-tj)…(t-tj+d)).

將pj(t)寫成Lagrange插值形式，整理得

(13)

(14)

將式(13)，(14)代入式(12)，得到高階重心有理插值公式，即

(15)

重心有理插值基于混合函數(shù)的思想，可以有效克服等距插值的不穩(wěn)定問題.

1.2.3 重心插值的微分矩陣 由式(15)，函數(shù)r(t)在節(jié)點t0，t1，…，tm處的v階導(dǎo)數(shù)為

(16)

由文獻(xiàn)[32]，一階和二階微分矩陣的計算公式為

(17)

利用數(shù)學(xué)歸納法，可以得到v階微分矩陣的遞推公式

(18)

1.3 一般迭代格式下的Cahn-Hilliard方程

(19)

式(19)中：k為迭代次數(shù).

1.4 Cahn-Hilliard方程的重心插值配點格式

φt+ε2φxxxx-s(x，t)φxx-g(x，t)φx=0.

(20)

首先，對空間域[a，b]和時間域[0，T]進(jìn)行網(wǎng)格剖分：a=x0

(21)

將式(21)代入式(20)，使其在點x0，x1，…，xM上成立，有

(22)

(23)

記φi(tj)=φ(xi，tj):=φi，j，得到φi(t)在點t0，t1，…，tN上的重心插值函數(shù)，即

(24)

將式(24)代入式(23)，使其在點t0，t1，…，tN上成立，有

(25)

φi=[φi，0，φi，1，…，φi，N]T，

si=diag(si，0，si，1，…，si，N)=diag(si(t0)，si(t1)，…，si(tN))，

gi=diag(gi，0，gi，1，…，gi，N)=diag(gi(t0)，gi(t1)，…，gi(tN))，i=0，1，…，M.

方程(25)可以簡記為矩陣形式，即

(26)

式(26)中:符號?代表矩陣的Kronecker積；IM，IN分別為M+1，N+1階單位矩陣；C(v)，D(v)(v=1，2，3，4)分別表示節(jié)點x0，x1，…，xM和t0，t1，…，tN上的重心插值v階微分矩陣.設(shè)

則式(26)可以表示為

[(IM?D(1))+ε2(C(4)?IN)-s(C(2)?IN)-g(C(1)?IN)]φ=0.

(27)

類似地，采用重心插值格式逼近s(x，t)和g(x，t)，得到Cahn-Hilliard方程在重心插值配點法下的迭代計算格式，即

[IM?D(1)+ε2(C(4)?IN)-diag(3(φ(k))2-1)(C(2)?IN)-
diag(6φ(k)((C(1)?IN)φ(k)))(C(1)?IN)]φ(k+1)=0.

(28)

1.5 插值誤差估計

由插值逼近誤差理論估計最后的誤差范圍，有助于檢驗算法在計算過程中的整體誤差結(jié)果.由文獻(xiàn)[33]的定理3.1，得到定理1.

定理1如果φ(x，t)∈Cn+1(Ω)，Ω=[a，b]×(0，T]∈R2是一個帶Lipschitz連續(xù)邊界的非空開集，φh(x，t)為φ(x，t)的重心Lagrange插值函數(shù)，則有如下誤差估計成立，即

(29)

(30)

(31)

證明：設(shè)xi(i=0，1，…M)為區(qū)間[-1，1]上的Chebyshev節(jié)點，x0=1，xM=-1，φ(x)∈Cm+1(-1，1)，由文獻(xiàn) [33]，有

(32)

由拉格朗日插值余項定理，有

(33)

式(33)中：v表示求導(dǎo)的次數(shù)，v=1，2，3，4.

當(dāng)M?1時，由Stirling公式知

(34)

結(jié)合式(32)～(34)，有

2.2.4 聚3-吲哚基環(huán)氧丙烷(Poly-IDPOs)、聚3-咔唑基環(huán)氧丙烷(Poly-CDPOs)

(35)

(36)

式(36)中：當(dāng)n較大時，c0，c1是與M無關(guān)的常數(shù).

(37)

(38)

對于二元函數(shù)φ(x，t)，有

|φ(x，t)-φh(xm，t)+φh(xm，t)-φh(xm，tn)|≤|φ(x，t)-φh(xm，t)|+|φh(xm，t)-φh(xm，tn)|≤

(39)

即有式(29)成立.取v=4，由式(38)即可得到式(30)，同理可得式(31).

從逼近誤差結(jié)果可知，該算法是指數(shù)收斂的，且微分算子的階數(shù)決定了代數(shù)方程的收斂階.

2 數(shù)值算例

給出兩個數(shù)值算例來驗證重心插值配點格式的高精度和有效性.記uh，u分別為數(shù)值解和真解向量，定義相對誤差和絕對誤差分別為

(40)

式(40)中：‖·‖∞表示無窮范數(shù).

2.1 算例1

為對比兩種重心插值配點格式的精度，選取三角函數(shù)精確解的方程作為測試實例.考慮Cahn-Hilliard方程，即

(41)

當(dāng)ε=0.3時，求解問題(41)，重心Lagrange插值配點格式和重心有理插值配點格式的誤差比較，如表1所示.由表1可知：節(jié)點數(shù)量在一定范圍內(nèi)增加時，兩種重心插值配點格式的計算誤差與節(jié)點數(shù)量之間呈指數(shù)下降，取較少節(jié)點數(shù)時，如M=20，25，N=10，10，重心Lagrange插值配點格式的精度高于重心有理插值配點格式的精度；隨著計算節(jié)點的進(jìn)一步增加，兩種重心插值配點格式的絕對誤差Ea和相對誤差Er都出現(xiàn)微小振蕩，重心有理插值配點的表現(xiàn)更平穩(wěn)，但仍都維持在10-11和10-10量級；兩種重心插值配點格式都具有很高的計算精度和較好的數(shù)值穩(wěn)定性；綜合比較可知，重心Lagrange插值配點格式的精度更高，重心有理插值配點格式的數(shù)值穩(wěn)定性更好.

表1 兩種重心插值配點格式的誤差比較Tab.1 Error comparison of two barycentric interpolation collocation formats

取M=25，N=15，重心Lagrange插值配點格式和重心有理插值配點格式的數(shù)值解圖及誤差分布圖，分別如圖1，2所示.圖1，2中:φh為數(shù)值解，δ為數(shù)值解和真解之間的誤差.

(a) 重心Lagrange插值配點 (b) 重心有理插值配點 圖1 兩種重心插值配點格式的數(shù)值解圖(M=25，N=15)Fig.1 Numerical solution graph of two barycentric interpolation collocation formats (M=25，N=15)

(a) 重心Lagrange插值配點 (b) 重心有理插值配點 圖2 兩種重心插值配點格式的誤差分布圖(M=25，N=15)Fig.2 Error distribution diagram of two barycentric interpolation collocation point formats (M=25，N=15)

由圖1，2可知：兩種重心插值配點格式的數(shù)值圖像幾乎一致，均逼近于真實解，但重心Lagrange插值配點格式的精度更高.

2.2 算例2

為了進(jìn)一步驗證數(shù)值格式的有效性，考慮 Cahn-Hilliard方程的離散能量函數(shù)，即

(42)

φ(-1，t)=-1，φ(1，t)=1，φ′(-1，t)=φ′(1，t)=0.

取M=30，N=30，在不同ε2取值下，兩種重心插值配點格式求解Cahn-Hilliard方程對應(yīng)的能量遞減圖及數(shù)值解圖，分別如圖3，4所示.

(a) 重心Lagrange插值配點(ε2=0.09) (b) 重心有理插值配點(ε2=0.09)

(c) 重心Lagrange插值配點(ε2=0.04) (d) 重心有理插值配點(ε2=0.04)圖3 不同ε2取值下兩種重心插值配點格式的能量遞減圖(M=30，N=30)Fig.3 Energy declining graphs of two barycentric interpolation collocation point formats under different values of ε2 (M=30， N=30)

(a) 重心Lagrange插值配點(ε2=0.09) (b) 重心有理插值配點(ε2=0.09)

(c) 重心Lagrange插值配點(ε2=0.04) (d) 重心有理插值配點(ε2=0.04) 圖4 不同ε2取值下兩種重心插值配點格式的數(shù)值解圖(M=30，N=30)Fig.4 Numerical solution graph of two barycentric interpolation collocation schemes under different ε2 values (M=30， N=30)

由圖3可知：兩種重心插值配點格式均滿足能量遞減規(guī)律，即隨著時間t遞增，能量函數(shù)E(φ)不斷遞減，最后達(dá)到穩(wěn)定狀態(tài)；當(dāng)ε2變小時，趨于穩(wěn)定狀態(tài)的時間會變長，能量值會變小.由圖4可知：不管ε2取何值，重心Lagrange插值配點插值效果都更好，數(shù)值解圖像更光滑，兩種重心插值配點格式圖像最終都趨于穩(wěn)態(tài).

3 結(jié)束語

利用重心Lagrange插值配點格式和重心有理插值配點的一般迭代格式求解Cahn-Hilliard方程，并給出了重心lagrange插值的逼近誤差.數(shù)值算例比較了兩種重心插值配點格式的精度和穩(wěn)定性.當(dāng)剖分節(jié)點數(shù)較少時，重心lagrange插值精度更高；當(dāng)剖分節(jié)點數(shù)較多時，重心有理插值配點穩(wěn)定性更好.最后，驗證了兩種格式均滿足能量遞減規(guī)律.接下來將進(jìn)一步研究重心插值配點推廣到高階非線性的問題.