李偉嘉，韓亦童，嚴 閱

（復旦大學 數(shù)學科學學院，上海 200433）

1 介 紹

這里定義所謂的“能量”函數(shù)：

和對應的常微分方程：

容易驗證，這個常微分方程的解滿足

這里常數(shù)C 由初值u0確定.如果在方程（2）兩邊乘以，并從t1到t2積分，即

那么有所謂的“能量”穩(wěn)定性：

如果考慮長時間計算的問題，也就是T 足夠大，可以考慮作尺度變換t→，從而可以考慮在［0，1］上帶小參數(shù)ε的問題：容易驗證uε（t）＝

如圖1所示，u（t）＝0是方程的不穩(wěn)定平衡點，u（t）＝±1是方程的兩個穩(wěn)定平衡點.當t→＋∞時，對于任意給定的u0，當u0＞0時，u（t）→1；u0＜0時，u（t）→－1.而且當ε越小，uε會變得越陡，求解的問題變得剛性和困難.因此穩(wěn)定的數(shù)值方法非常重要.

上述這類方程的特點是有比較弱的（線性）增長項，另外有比較強的（三次）非線性控制項.當初值比較小的時候，即｜u0｜?1，解隨著時間會指數(shù)增長，由于有非線性項的控制作用，增長速度單調下降趨于0，解也趨于穩(wěn)定.相反，當初值比較大的時候，即｜u0｜≥1，解隨著時間會指數(shù)衰減，由于有線性項的增長作用，衰減速度單調下降趨于0，解同樣趨于穩(wěn)定.

在數(shù)值計算上，顯然如果用最簡單的顯式Euler格式（或其他的顯式格式）是不適合的，為了計算穩(wěn)定，時間步長要選取的非常小，特別是對剛性問題（或小ε）的情況；另一方面，如果用隱式格式，會需要面臨代數(shù)方程的根不唯一的問題［1－2］.最近在研究薄膜外延生長（Molecular Beam Epitaxy，MBE）模型的數(shù)值方法過程［3－11］中，Eyre的凹凸分離的方法［12］被進行了深入的研究和分析，并被推廣應用到Darcy－Stokes流［13］、Cahn－Hilliard模型［14］、相場晶體模型［15］和帶指數(shù)非線性項的波方程［16］等.凹凸分離的基本思想是在能量泛函中，把能量分解為凸項和凹項的組合，對凸能量項用隱式格式處理，對凹能量項用顯式處理，這樣構造的格式具有能量穩(wěn)定、求解非線性代數(shù)方程解唯一等優(yōu)點，但是這個方法對于如何構造兩階或更高階的格式到現(xiàn)在為止只有少量的工作.本文回到最基本的常微分方程，結合常微分方程經(jīng)典的Gear格式（BDF格式）和凹凸分離思想來構造計算格式.BDF格式是一類對剛性問題非常有效的典型的數(shù)值格式［1－2］，這里我們修改BDF格式，通過對能量函數(shù)中的凹項作顯式或外推處理，對凸項直接用隱式處理.本文給出了相應的穩(wěn)定性和收斂性分析，并把這個格式用于求解Allen－Cahn方程，數(shù)值結果表明這個格式和想法對于求解非線性偏微分方程（Partial Differential Equation，PDE）也是有效的.

圖1 常微分方程解的示意圖Fig.1 Sketch map of solutions to the ordinary differential equations

2 兩個凹凸分離的數(shù)值格式

按照凹凸分離的思想，考慮方程（2）的數(shù)值格式.設N 是正整數(shù)，0＝t0＜t1＜…＜tN＝T 是［0，T］上的均勻剖分，τi＝ti－ti－1，i＝1，2…，N，為方便起見我們假設步長是等步長的，也就是τ＝τi＝.對于該問題，能量函數(shù)（1）可以分解為兩個凸函數(shù)的組合，即

注意到這種凹凸的分解是不唯一的.根據(jù)經(jīng)典的凹凸分離格式的構造方法，我們將E′c（u）處理為隱式項，處理為顯式項，構造如下一階和二階數(shù)值格式：

1）一階格式

2）二階格式

這里un為tn時刻的計算解.在上述兩個格式中，我們基于Gear格式（BDF 格式）來盡可能得到大的絕對穩(wěn)定區(qū)域，同時結合外推來處理顯式項（Adams格式）使得整體精度一致.這個格式在文獻中也稱為IMEX（IMplicit－EXplicit）隱顯格式.在常見的IMEX 格式中通常是線性項隱式處理，非線性項顯式處理來使得計算方便.這里我們對非線性項用隱式處理，線性項顯式處理是為了穩(wěn)定性考慮，特別是長時間穩(wěn)定性.當然，在每步的計算過程中，還需要求解一個非線性代數(shù)問題，由于隱式部分是嚴格凸函數(shù)的導數(shù)，所以像牛頓方法、非線性共軛梯度方法等都有好的結果.特別是對我們的例子，三次方程可以有解表達式，而且解是唯一的，對整體計算量增加不大.

為方便起見，我們引入幾個差分記號δ，δ2和D：

3 穩(wěn)定性分析

在這一節(jié)我們將證明上述兩個計算格式（9），（10）是適定的，且能量穩(wěn)定的.

定理1 求解方程（2）的格式（9）是適定的，并且保持能量遞減性和數(shù)值解的有界性，即

證 對于格式（9），設

則f（u）的原函數(shù)

由于

因此F（u）為凸函數(shù)，而且嚴格凸的，所以格式（9）唯一可解.然后，在一階格式（9）兩端同乘以δun＋1，則

由于Ec（u）和Ee（u）都是凸函數(shù)，故有下列不等式成立：

因此，

所以能量遞減性成立：

由E（u）定義，易知

從而解的有界性成立.

定理2 求解方程（2）的格式（10）是適定的，并且在數(shù)值能量函數(shù)

下，如果滿足τ≤2，則保持能量遞減性和數(shù)值解的有界性，即

證 仿照定理1證明，易知格式滿足唯一可解性.在格式（10）兩端同乘以δun＋1，則

結合上述兩式，并由數(shù)值能量定義（19）式和條件τ≤2，有

從而得到數(shù)值能量的單調遞減性.

再由遞推關系式

和數(shù)值能量的定義（19）式得到

從而解的有界性成立.

注1 上述定理成立需要滿足條件τ≤2，如果在格式左邊加入穩(wěn)定項τ（un＋1－un），使之變?yōu)樾露A格式

則新二階格式解除了對時間步長τ的限制，無條件滿足能量遞減性.

4 收斂性

這里我們將證明格式（9）是一階收斂的，而格式（10）是兩階收斂的.定義tn時刻的誤差為

在證明過程中，需要下面基本的不等式（引理1）和Gronwall不等式（引理2），證明分別見文獻［6］和［17］.

引理1［6］對任意的向量a，b∈Rn，都成立

引理2［17］對于固定的T，假設是非負序列，滿足≤C1，其中C1不依賴于τ，N，Nτ＝T，如果?τ＞0，滿足

其中C2＞0是不依賴于τ和N 的常數(shù).那么?τ＞0，有

定理3 假設考慮的時間區(qū)間為t∈［0，T］，u（t）為連續(xù)問題（2），（3）的解，uu為離散數(shù)值格式（9）在t＝tn時刻的解，如果滿足τ≤1／4，則我們有如下誤差估計：

證 首先，計算截斷誤差

由Cauchy－Schwarz不等式得到

因此截斷誤差可以有估計

結合方程（2）和數(shù)值格式（9），得到誤差方程

上式兩端同乘以2τen＋1，則

運用Young不等式和引理1，得到如下估計：

將上述估計代回（29）式，并舍棄一些非必要項，有

對k＝0，1，…，n（n≤N）進行累加，有

這里

由于τ≤1／4，故應用離散的Gronwall不等式，得到

兩邊開方即得證.

為方便起見，在二階格式收斂性證明中，我們引入G 范數(shù)記號.回憶一下G 范數(shù)的定義：設w＝（u，u）T∈R2，定義G 范數(shù)

顯然矩陣G 是對稱正定陣，故上述G 范數(shù)可定義.由于

顯然G1是半正定矩陣，故有

此外，任給vi∈R，i＝0，1，2，有

其中w0＝（v0，v1）T，w1＝（v1，v2）T.

定理4 假設考慮的時間區(qū)間為t∈［0，T］，u（t）為連續(xù)問題（2），（3）的解，un為離散數(shù)值格式（10）在t＝tn時刻的解，如果滿足τ≤1／16，則我們有如下誤差估計：

證 類似于一階格式.首先，計算截斷誤差

通過簡單的估計和Schwarz不等式得到

類似地，有

因此截斷誤差為

兩階格式（10）的誤差方程為

上式兩端同乘以2τen＋1，根據(jù)Young不等式、引理1和（34）式，有

其中wn＋1＝（en，en＋1）T.對k＝0，1，…，n（n≤N）進行累加，舍棄不必要項，并運用不等式（33），有

這里

由于τ≤1／16，故應用離散的Gronwall不等式，得到

兩邊開方即得證.

5 數(shù)值結果

首先格式（9）和格式（10）都需要求解三次方程

我們已經(jīng)知道由這兩個格式得到的解是唯一的實根，可以利用如下公式來求解：

在極端情況下，可能會出現(xiàn)大數(shù)相減導致有效位數(shù)損失的情況［18］，但由于這里考慮的還是低階格式，所以舍入誤差不在這里討論.

下面我們先給出初始值的計算結果，從圖2（a）可以看到對于不同的初始值（u0＝±0.01，±0.5，±1.5），計算解和精確解結果一致.要提到的是上面的兩階格式的初始值u1是用最簡單的顯式Euler格式計算得到的，可以證明這樣的取法其整體誤差還是兩階的.另外我們要說明的是，對于長時間計算T 較大的情況，當τ比較大的時候，在初始狀態(tài)，解有可能出現(xiàn)“過沖”現(xiàn)象，例如在圖2（b）中可以看到，當初始值u0＝0.5，當T＝8，N＝8的時候，對于二階格式，計算解會大于1，這個現(xiàn)象當增加等分次數(shù)時就會消失.這也說明定理3和定理4所要求的條件是充分而非必要的.

圖2 不同初始值的精確解和計算解（a）；T＝8，N＝8時候的數(shù)值解（b）Fig.2 （a）Exact and numerical solutions with different initial values；（b）Numerical solutions when T＝8，N＝8

另外以初始值u0＝0.5為例，看看兩個計算格式是否分別為一階和兩階精度.我們計算當T＝1時刻的誤差收斂階，從表1可以看出兩格式分別是一階和兩階收斂的.對于其他的初值可以得到類似的結論.這里要提到的是關于初始值的選取，當步長τ比較小的時候，我們可以用顯式Euler格式來得到二階格式的第二個初始值u1，但是當τ較大或者ε 很小的時候，如果用顯式Euler格式會得到定性上錯誤的解，也就是正的初始值會趨向負的平衡點u＝－1，這是由于u1可能會是個負值.因此我們建議用本文的一階格式來計算二階格式u1的值.另外對于長時間計算（大的T），本文的格式都有非常好的結果，與預期一致，解都最終收斂到平衡點.對非常小的ε，這兩個格式也是非常穩(wěn)定的.這里就不再附上數(shù)值結果.

表1 粗－細網(wǎng)格誤差的收斂階，等分數(shù)從26 到213Tab.1 Convergence order of errors on coarse and fine grids，with Nfrom 26 to 213

6 對Allen－Cahn方程的離散

下面考慮用我們上面構造的例子來求解Allen－Cahn方程：

這里時間方向上采用之前所構造的一階和二階格式，空間方向上采用經(jīng)典的五點差分格式來求解方程：

1）一階格式

2）二階格式

取Ω＝［0，1］×［0，1］，T＝1，ε2＝0.1，方程真解為u（x，t）＝e－tsin（πx1）sin（πx2），g 由u 計算得到.非線性方程采用牛頓法求解，初始迭代值取前兩步外推值，即，大多數(shù)情況下，這種取法會比初始迭代值取前一步值少迭代一次，從而縮短計算時間.停機準則為＜10－9.首先，驗證誤差在時間方向上的收斂階.取τ＝h，τ為時間步長，h為空間步長，err1，err2分別為一階和二階格式的誤差項.計算結果見表2.然后，驗證誤差在空間方向上的收斂階.我們熟知五點差分格式在空間上是二階收斂的，因此在一階格式中取τ＝h2，使得時間方向誤差和空間方向誤差同階，err 為誤差項.計算結果見表3.從表2，表3中可以看出，我們所構造的格式在時間和空間方向上都得到了滿意的收斂結果.

表2 Allen－Cahn方程時間方向的收斂階分析Tab.2 Convergence analysis in time for Allen－Cahn equations

表3 Allen－Cahn方程空間方向的收斂階分析Tab.3 Convergence analysis in space for Allen－Cahn equations

［1］Hairer E，Wanner G.Solving ordinary differential equations Ⅰ：Nonstiff and problems［M］.Berlin：Springer－Verlag，1993.

［2］Hairer E，Wanner G.Solving ordinary differential equations Ⅱ：Stiff and differential－algebraic problems［M］.Berlin：Springer－Verlag，1996.

［3］Li B，Liu J G.Thin film epitaxal with or without slope selection［J］.Euro J Appl Math，2003，14（6）：713－743.

［4］Li B，Liu J G.Epitaxial growth without slope selection：Energetics，coarsening and dynamic scaling［J］.J Nonlinear Sci，2004，14（5）：429－451.

［5］Xu C J，Tang T.Stability analysis of large time－stepping methods for epitaxial growth models［J］.SIAM J Numer Anal，2006，44（4）：1759－1779.

［6］夏 茜，陳文斌，劉建國.關于薄膜外延生長模型隱式全離散格式的誤差分析［J］.高校計算數(shù)學學報，2012，34（1）：30－51.

［7］Wang C，Wang X M，Wise S M.Unconditionally stable schemes for equations of thin film epitaxy［J］.Discrete Contin Dyn Syst，2010，28（1）：405－423.

［8］Chen W B，Wang Y Q.A mixed finite element method for thin film epitaxy［J］.Numer Math，2012，122（4）：771－793.

［9］Qiao Z H，Sun Z Z，Zhang Z R.The stability and convergence of two linearized finite difference schemes for the nonlinear epitaxial growth model［J］.Numer Meth Part Differ Equ，2012，28（6）：1893－1915.

［10］Chen W B，Conde S，Wang C，et al.A linear energy stable scheme for a thin film model without slope selection［J］.J Sci Comput，2012，52（3）：546－562.

［11］Chen W B，Wang C，Wang X M，et al.A linear iteration algorithm for a second order energy stable scheme for a thin film model without slope selection［J］.J Sci Comput，2014，59（3）：574－601.

［12］Eyre D J.Unconditionally gradient stable time marching the Cahn－Hilliard equation［M］∥Computational and Mathematical Models of Microstructhral Evolution，Vol.52.Warrendale：Materials Research Society，1998：1686－1712.

［13］Chen W B，Gunzburger M，Sun D，et al.Efficient and long－time accurate second－order methods for Stokes－Darcy system［J］.SIAM J Numer Anal，2013，51（5）：2563－2584.

［14］Collins C，Shen J，Wise S M.An efficient，energy stable scheme for the Cahn－Hilliard－Brinkman system［J］.Commun Comput Phys，2013，13（4）：929－957.

［15］Baskaran A，Hu Z Z，Lowengrub J S，et al.Energy stable and efficient finite－difference nonlinear multigrid schemes for the modified phase field crystal equation［J］.J Comput Phys，2013，250：270－292.

［16］Wang L D，Chen W B，Wang C.An energy－conserving second order numerical scheme for nonlinear hyperbolic equation with an exponential nonlinear term［J］.J Comput Appl Math，2015，280：347－366.

［17］Layton W.Introduction to the numerical analysis of incompressible viscous flows［M］.Philadelphia：SIAM，2008.

［18］Higham N J.Accuracy and stability of numerical algorithms，Second edition［M］.Philadelphia：SIAM，2002.