潘晗霜，羅智文，王淑芬

(復旦大學 數學科學學院，上海 200433)

1 基礎知識

本文考慮無斜率選擇的分子束外延模型[1-3].設Ω是2上的有界區(qū)域，u和w滿足如下方程：

?tu-ε2Δw+

(1)

w+Δu=0.

(2)

并考慮邊值條件：

?nu=?nw=0 在區(qū)域?Ω×[0,T]上，

(3)

和初值條件：

u(x,0)=u0(x),

(4)

容易驗證u，w具有守恒性質，這里(·)表示L2內積：

(u,1)=(w,1)=0.

(5)

對于上述方程的能量泛函定義為

(6)

將能量中的第一項稱為Ehrlich-Schwoebel能量項，該項具有非常重要的物理意義.

文獻[4]中通過擾動分析的方法對方程進行分析，得出了對應方程初邊值問題的適定性.除此之外，還得出了很多有趣的粗化性質，例如Li等[5]證明了當ε→0時，對于長時間t，能量下界被O(-lnt)所限制.由于梯度流粗化過程的長時間性，對外延方程提出的數值方法需要具有長時間穩(wěn)定以及數值精度高的性質.

關于能量穩(wěn)定的數值格式，受Eyre[6]所做工作的影響，在研究中涌現出了一類重要的數值處理方法——凹凸分離的方法，文獻[7]中首次對外延生長模型提出了時間一階無條件能量穩(wěn)定的凹凸分離格式，但格式的非線性導致數值求解比較具有挑戰(zhàn).針對這個問題，文獻[8]提出了一種基于凹凸分離的線性格式.通過對能量巧妙地分離，使得非線性項對應的能量位于凸部，并且通過顯式處理，使得得到的格式為線性格式，極大地簡化了數值計算難度.

為了實現兩階時間精度，文獻[9]中基于同樣的凹凸分離方法，通過使用修正的Crank-Nicolson格式得到兩階數值格式，并且格式保持無條件能量穩(wěn)定.類似地，還可以使用后向差分(Backward Differentiation Formula, BDF)格式構造兩階精度凹凸分離格式.但是，這兩個格式都是非線性的，使得數值實現變得比較困難.為了處理非線性的問題，文獻[10]中提出了一種線性迭代的能量穩(wěn)定格式，通過迭代算法來有效地處理高度非線性的格式.事實上，對MBE的能量穩(wěn)定格式還有很多研究[11-14].

文獻[15]中提出了一種無條件能量穩(wěn)定的線性格式，其格式設計為：

Aτ(vh)=0,

(7)

(8)

使用后向差分格式處理對時間的導數項，而對非線性項使用顯式外推格式，同時，為了保證格式能量穩(wěn)定，還人為添加了一個兩階時間精度O(τ2)的正則項AτΔ2(un+1-un).

本文在文獻[15]提出的格式基礎上，考慮了新的正則項設計以及非線性項外推，從而得到了一種新的兩階數值格式，這個數值格式同樣可以線性求解且能量穩(wěn)定，本文將給出相應的穩(wěn)定性和收斂性分析，最后使用數值實驗驗證了理論分析中得出的各項性質.

2 新的兩階數值格式

首先，使用有限元方法對方程進行空間上的離散： 用Th={K}表示Ω上的擬一致三角剖分，這里h表示網格劃分大小的離散參數.記P=Pq(K)為在單元網格K上定義的不高于q階的多項式函數集合，并定義分片多項式函數空間Xh≡{v∈X∩C0(Ω)|v|K∈P(K)，?K∈Th}?X.這里X表示在Ω上積分為零的Sobolev函數空間X={v∈Hq(Ω)|(v,1)=0}.

((Rhφ-φ),χ)=0, (Rhφ-φ,1)=0 ?χ∈Xh.

(9)

(10)

(11)

對于n≥1時的中間步格式設置，我們將文獻[15]的格式中的正則項由Aτ(vh)修改為(vh)，對于非線性項外推，我們把其由形式修改為從而得到新的中間步格式： 給定和尋找使得

(12)

(13)

注1對正則項修改是非常自然的，使得兩個時間方向離散的處理統(tǒng)一起來，對于非線性項，現在的處理方式也是非常自然的，它是對非線性項函數做外推，來取代對原來解的外推.

3 穩(wěn)定性分析

首先，定義離散的能量泛函：

(14)

在穩(wěn)定性證明中，我們會用到文獻[15]中的如下引理.

引理1給定函數φ1,φ2∈X和v∈X，我們定義函數Φφ1,φ2,v(s): [0,1]→，

(15)

則Φφ1,φ2,v滿足

(16)

對數值格式(12)，(13)有如下能量穩(wěn)定結果.

定理1給定A≥1，數值格式(12)，(13)具有能量下降的性質，即

(17)

證 首先從式(13)可以推出

(18)

(19)

接下來我們依次估計式(19)右邊的項.對于第1項，由Cauchy-Schwarz不等式可得

(20)

對于式(19)右邊第2項，由式(13)可得

(21)

對式(19)右邊第3項存在估計：

(22)

直接引用文獻[15]中的結果，對Ehrlich-Schwoebel能量項存在估計：

(23)

由此成立下面的不等式：

(24)

下面我們依次估計式(24)右邊積分中的3項S1，S2和S3.首先對于S1，利用引理1的結論以及Cauchy-Schwarz不等式：

-|

(25)

類似于對S1的估計，我們可以得到：

(26)

而對于S3，注意到對任意正實數a，b，成立a(a+b)<(1+a2)(1+b2)，因此，

(27)

綜合上面3項的估計，可得

(28)

對于式(28)中的梯度項，可以利用Cauchy-Schwarz不等式做放縮：

(29)

(30)

綜合式(19),(20),(21),(22),(28)可得到下式：

(31)

(32)

從而對修正的能量格式成立能量下降的性質.

4 收斂性分析

首先，通過引用文獻[16]的結論，對于任意φ∈Hq+1(Ω)∩X的Ritz投影Rh，成立如下估計：

(33)

這里我們將證明格式是兩階收斂的.用(u,w)表示方程(1)，(2)的真實解，為了使后面分析成立，需要對(u,w)作出如下正則性假設：

u∈L∞(0,T;Hq+1)∩H1(0,T;Hq+1)∩H2(0,T;H1)∩W2,∞(0,T;L2)∩H3(0,T;L2),
w∈L∞(0,t;Hq+1)∩H1(0,T;H1).

(34)

下面將證明格式的誤差收斂階.

這里C2為任意大于0的常數，則有如下的誤差估計：

證 首先，定義如下的誤差函數：

(35)

(36)

對于中間步,將方程和數值格式(12),(13)相減可得誤差方程:

(37)

其中：

根據映射Rh的定義可知： (χ)=0,?χ∈Xh，由此式(37)可改寫為： 當n≥1時，對任意vh∈Xh,φh∈Xh，成立

(38)

(39)

(40)

這里用到了一個變換：

(41)

運用這個記號，式(40)中左邊第1項可以寫為：

(42)

(43)

接下來估計式(40)的右邊項.對于前兩項，運用Ritz投影的性質(33)以及Cauchy-Schwarz不等式可得

(44)

以及

(45)

對于第3項和第4項的分析，同樣可參考文獻[15]的結論：

(46)

(47)

為了估計式(40)右邊第5項，注意到：

(48)

以及

(49)

根據Taylor展開和Cauchy-Schwarz不等式，結合式(48)和(49)，對于式(40)右邊第5 項成立：

(50)

(Nn+1，

(Nn+1,

(51)

將式(42)～(51)帶回誤差方程(40)，得到:

(52)

對式(52)做從n=1到n=m求和，再在兩邊乘上2τ，可得

(53)

容易驗證

根據式(33)和(36)的結果可得

(54)

(55)

結合式(55)以及式(14)的結果即可得到定理2的結論.

5 數值結果

這一節(jié)中，我們將通過幾個數值實驗來驗證前面提到的一些性質，首先，我們驗證格式的收斂階，在下面的例子中我們都是使用文章第1節(jié)提出的數值格式，這里設定區(qū)域Ω=[0,1]2和ε2=0.05，方程的真值解為

ue(x,y,t)=cos(πx)cos(πy)e-t.

(56)

為了滿足方程以及邊界條件，我們在方程的右邊添加上外力項：

?tue+

這里，

表1 空間收斂性

表2 時間收斂性

下面，我們將進行數值模擬來驗證結果是否具有能量下降等性質.對于方程本身的參數，我們取定ε2=0.005，而將方程的初值選取為每個點均為-0.05至0.05之間的隨機數，設定空間區(qū)域Ω為[0，12.8]×[0，12.8]的正方形區(qū)域.對于模型的參數設定，首先就空間的離散，考慮到計算的效率我們設定N=128，為了觀察方程解的長時間性質，取定T=20000，而對時間步長，為了快速計算，我們在不同的時間區(qū)間選取了不同的時間步長： 當t≤200時，時間步長τ=0.004；當2002000時，取時間步長為0.15.

我們首先觀察該數值格式解出的解在不同時刻的形態(tài)，圖1中所示為解在t=1,500,1000,5000,10000,20000時刻的形態(tài).

接下來，我們來驗證解是否具有能量下降的性質.圖2展示了能量隨時間的變化情況，從該圖中我們可以明顯觀察到能量是隨時間下降的，也就從數值上驗證了方法的能量穩(wěn)定性.

定義表面粗度R(t)以及平均梯度W(t)分別為:

(57)

從文獻[1，5] 中我們已經知道，對于沒有斜率選擇的分子束外延方程，E(t),R(t)和W(t)分別滿足如下性質：

(58)

以下，我們將對數值結果分別驗證其是否滿足上面提到的3條性質.為此，圖3(a)顯示了能量E(t)對lnt的變化情況，圖中紅色虛線表示線性擬合的結果，這里我們擬合出的結果為E(t)=-38.1890 lnt-60.1026，而對于表面粗度R(t)和時間的關系，圖3(b)顯示了ln(R(t))相對于lnt的變化情況，同樣這里給出了線性擬合的結果，即R(t)=0.4454t0.482 4，最后對于平均梯度W(t)，可以查看圖3(c)的結果，這里擬合出的結果為W(t)=2.3245t0.243 8.從這些結果我們可以發(fā)現，該格式計算出的數值結果與文獻[1，5]中的結論相符.