A stabilized second order exponential time differencing multistep method for thin film growth model without slope selection
Wenbin Chen, Weijia Li, Zhiwen Luo, Cheng Wang, Xiaoming Wang
AA STABILIZED SECOND ORDER EXPONENTIAL TIMEDIFFERENCING MULTISTEP METHOD FOR THIN FILMGROWTH MODEL WITHOUT SLOPE SELECTION
WENBIN CHEN ∗ , WEIJIA LI † , ZHIWEN LUO ‡ , CHENG WANG § , AND
XIAOMINGWANG ¶ Abstract.
In this paper, a stabilized second order in time accurate linear exponential timedifferencing (ETD) scheme for the no-slope-selection thin film growth model is presented. An artificialstabilizing term Aτ ∂ ∆ u∂t is added to the physical model to achieve energy stability, with ETD-basedmulti-step approximations and Fourier collocation spectral method applied in the time integral andspatial discretization of the evolution equation, respectively. Long time energy stability and detailed (cid:96) ∞ (0 , T ; (cid:96) ) error analysis are provided based on the energy method, with a careful estimate of thealiasing error. In addition, numerical experiments are presented to demonstrate the energy decayand convergence rate. Key words. epitaxial thin film growth, exponential time differencing, long time energy stability,convergence analysis, second order scheme
AMS subject classifications.
1. Introduction.
Consider the continuum model of no-slope-selection thin filmepitaxial growth, which is the L gradient flow of the following energy functional: E ( u ) = (cid:90) Ω ε | ∆ u | −
12 ln(1 + |∇ u | ) d x , (1.1)where Ω = [0 , L ] , u : Ω × [0 , T ] → R is a scaled height function of thin film withperiodic boundary condition, and ε > E ( u ) rep-resents the surface diffusion, and the logarithm term models the Ehrlich-Schowoebel(ES) effect which describes the effect of kinetic asymmetry in the adatom attachment-detachment; see [13, 26, 27, 37]. Consequently, the growth equation is the L gradientflow of (1.1): ∂u∂t = − ε ∆ u − ∇ · (cid:18) ∇ u |∇ u | (cid:19) , x ∈ Ω , t ∈ (0 , T ] . (1.2)Taking the inner product with (1.2) by ∂u∂t , we obtain the energy decay property forthe continuous system: d E d t = −(cid:107) u t (cid:107) L . ∗ Shanghai Key Laboratory of Mathematics for Nonlinear Sciences, School of Mathematical Sci-ences; Fudan University, Shanghai, China 200433 ( [email protected] ) † School of Mathematical Sciences; Fudan University, Shanghai, China 200433( ) ‡ School of Mathematical Sciences; Fudan University, Shanghai, China 200433( ) § Mathematics Department; University of Massachusetts; North Dartmouth, MA 02747, USA( corresponding author: [email protected] ) ¶ Department of Mathematics, Southern University of Science and Technology, Shenzhen, China518055, and Fudan University, Shanghai, China 200433, and Florida State University, Tallahassee,FL 32306, USA ( [email protected] ) 1 a r X i v : . [ m a t h . NA ] J u l W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
The equation (1.2) is referred to as the no-slope-selection (NSS) equation, since (1.2)predicts an unbounded mound slope m ( t ) = O ( t ) [16, 27]. On the other hand, theslope-selection (SS) case refers to the case when the logarithm term −∇ · (cid:16) ∇ u |∇ u | (cid:17) in (1.2) is replaced by ∇ · ( |∇ u | ∇ u ), which predicts that mound-like or pyramidstructures in the surface profile tend to have a uniform, constant mound slope [27].The solution u ( x , t ) is mass conservative, i.e., (cid:90) Ω u ( x , t ) d x = (cid:90) Ω u ( x ,
0) d x , t > , (1.3)due to the periodic boundary condition. For simplicity of presentation, herein weassume u to have zero mean-value. Scaling laws that characterize the surface coars-ening during the film growth have been a physically interesting problem; see [16, 25,26, 28, 33]. As for the continuum model, the well-posedness of the initial-boundary-value problem for both SS and NSS equation has been given by [26, 27] through theperturbation analysis and Galerkin spectral approximations. Li and Liu [28] provedthat the minimum energy scales as log ε for small ε >
0, and the average energy isbounded below by O ( − log t ) for large t . This implies that long time energy stabilityis required to simulate the coarsening process.One popular way to construct energy stable numerical scheme is to split theenergy functional into convex and concave parts, then apply implicit and explicittime stepping algorithms to the corresponding terms [14], respectively; see the firstsuch numerical scheme in [41] for the molecular epitaxy growth model. Since then,there have been various works dedicated to deriving high order and energy stableschemes under this idea, such as [8, 9, 10, 32, 35, 36, 39]. In particular, a first orderin time linear scheme was proposed in [8] to solve the NSS equation, with the energystability established. On the other hand, the concave term turns out to be nonlinearin this approach, and there has been a well-known difficulty of the convex-splittingapproach to construct unconditionally stable higher-order schemes for a nonlinearconcave term. Many efforts have been devoted to overcome this subtle difficulty, suchas introducing an artificial stabilizing term in the growth equation to balance theexplicit treatment of the nonlinear term; see [15, 30, 31, 43]. In addition, there havebeen some other interesting approaches for the stability of a numerically modifiedenergy functional, such as the invariant energy quadratization (IEQ) [44] and thescalar auxiliary variable (SAV) methods [40].Other than these direct approaches to establish the energy stability, some otherideas have been reported in recent years to obtain higher order temporal accuracywith explicit treatment of nonlinear terms, such as the time exponential time differ-encing (ETD) approach. In general, an exact integration of the linear part of theNSS equation is involved in the ETD-based scheme, followed by multi-step explicitapproximation of the temporal integral of the nonlinear term [3, 4, 12, 19, 20]. Anapplication of such an idea to various gradient models has been reported in recentworks [11, 21, 22, 23, 24, 42, 45], with the high order accuracy and preservation ofthe exponential behavior observed in the numerical experiments. At the theoreticalside, some related convergence analyses have also been reported, while a rigorous en-ergy stability analysis has only been provided for the first order scheme [21]. For thesecond and higher order numerical schemes, a theoretical justification of the energybound has not been available.In this paper, we work on a second order in time accurate ETD-based scheme, withFourier pseudo-spectral approximation in space, and the energy stability analysis is ETDMs2 for no-slope-selection thin film model Aτ ∂ ∆ u∂t is added in the growthequation, where A is a positive constant and τ is the time step size. Also, we apply alinear Lagrange approximation to the nonlinear term as in [21]. This approach enablesus to perform a careful energy estimate, so that a decay property for a modifieddiscrete energy functional is proved. This in turn leads to a uniform in time boundfor the original energy functional. In our knowledge, this is the first such result for asecond order ETD-based scheme, with the primary difficulty focused on the explicitextrapolation for the nonlinear terms. Moreover, we provide a novel convergenceanalysis for the proposed scheme. Instead of analyzing the operator form of thenumerical error function, here we start from the continuous in time ODE systemsatisfied by the error function, which is derived from the corresponding equations ofthe numerical solution. With a careful treatment of the aliasing error and H estimateof the numerical solution, we are able to derive an (cid:96) ∞ (0 , T ; (cid:96) ) error estimate for theproposed scheme.There have been quite a few physically interesting quantities that may be obtainedfrom the solutions of the NSS equation (1.2), such as the energy, average surfaceroughness and the average slope. A theoretical analysis in [28] has implied a lowerbound for the energy dissipation as the order of O ( − ln t ), and an upper bound for theaverage roughness as the order of O ( t / ), which in turn implies an O ( t / ) order forthe average slope. In particular, the O ( − ln t ) energy dissipation scale leads to a longtime scale nature of the NSS equation (1.2), in comparison with the O ( t − / scale forthe slope-selection version and the standard Cahn-Hilliard model. In addition, a lowerbound for the energy, estimated as γ := L (cid:0) ln(4 ε L − ) − ε π + 1 (cid:1) as derived in [8,9, 41], indicates an intuitive O ( ε − ) law for the saturation time scale for ε (cid:28) min L ,because of the O ( − ln t ) energy dissipation scale. All these facts have demonstratedthe necessity of energy stable numerical approach to accurately capture the physicallyinteresting quantities, since the energy stability turns out to be a global in time naturefor the numerical scheme. In the extensive numerical experiments, the proposedstabilized ETD scheme has produced highly accurate solutions, which obtain the longtime asymptotic growth rate of the surface roughness and average slope with relativeaccuracy within 4%.The rest of this article is organized as follows. In Section 2 we present the fullydiscrete numerical scheme. The numerical energy stability is proved in Section 3,followed by (cid:96) ∞ (0 , T ; H h ) and (cid:96) ∞ (0 , T ; H h ) bounds of the numerical solution. Subse-quently, an (cid:96) ∞ (0 , T ; (cid:96) ) error analysis is given by Section 4, consisting of two lemmasconcerning the error of the nonlinear term. In addition, numerical experiments areprovided in Section 5, including temporal convergence test and simulation results ofthe scaling laws for energy, average surface roughness and average slope. Finally,some concluding remarks are given by Section 6.
2. The stabilized second order in time ETD mulstistep scheme (sEDTMs2).
Some definitions in [1] are recalled. Consider Ω = [0 , L ] d . We define W m,p (Ω) := { v ∈ L p (Ω) | D α v ∈ L p (Ω) for 0 ≤ | α | ≤ m } , where α = ( α , . . . , α d ) is a d -tuple of non-negative integers with | α | = (cid:80) di =1 α i , and D α = (cid:81) di =1 D α i x i . For simplicity, we denote (cid:107) · (cid:107) W m,p (Ω) as (cid:107) · (cid:107) m,p , and the related semi-norm as | · | m,p . In particular, if p = 2,we set W m,p (Ω) as H m (Ω), and (cid:107) · (cid:107) m, as (cid:107) · (cid:107) H m .In this paper, we focus on the case when d = 2 and follow the notations used incollocation spectral method; see [5, 6, 17, 18, 21, 38]. Let N be a positive integer, W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG Ω N be a uniform 2 N × N mesh on Ω, with (2 N + 1) grid points ( x i , y j ), where x i = ih , y j = jh with h := L N , 0 ≤ i, j ≤ N . Denote u e as the exact solution of(1.2), u ( t ) = u e ( · , t ) | Ω N be the restriction of the exact solution on Ω N , τ be the timestep size TN t , and t i = iτ for 0 ≤ i ≤ N t . Let H mper (Ω) = { v ∈ H m (Ω) | v is periodic } , M N be the space of 2D periodic grid functions on Ω N , and B N be the space oftrigonometric polynomials in x and y of degree up to N . In this paper, we denote C one generic constant which may depend on ε , the solution u , the initial value u andtime T , but is independent of the mesh size h and time step size τ .Let us firstly recall the following Berstein inequality introduced in [38, p. 33,Lemma 2.5]: Lemma 2.1.
For any u ∈ B N and ≤ p ≤ ∞ , we have (cid:107) ∂ mx u (cid:107) L p (Ω) ≤ CN m (cid:107) u (cid:107) L p (Ω) , m ≥ . (2.1)Now we introduce an interpolation operator I N onto B N that reserves the functionvalue on (2 N + 1) grid points, i.e., ( I N f )( x i , y j ) = f ( x i , y j ) for 0 ≤ i, j ≤ N :( I N f )( x, y ) = N (cid:88) k,l = − N +1 ( ˆ f c ) k,l exp (cid:18) π i L ( kx + ly ) (cid:19) , with i = √− , (2.2)where the coefficients { ( ˆ f c ) k,l } are given by the discrete Fourier transform of the 4 N grid points: ( ˆ f c ) k,l = 14 N N (cid:88) i,j =1 f ( x i , y j ) exp (cid:18) − π i L ( kx i + ly j ) (cid:19) . (2.3)For any f ∈ M N , denote ˜ f = I N f as the continuous extension of f . When f and ∂ α f ( | α | ≤ m ) are continuous and periodic on Ω, I N has the following approximationproperty ([7, Theorem 1.2, p. 72]): (cid:107) ∂ k f − ∂ k I N f (cid:107) L ≤ Ch m − k (cid:107) f (cid:107) H m , for 0 ≤ k ≤ m, m > d , (2.4)for dimension d . Also, the following H m bound of I N is excerpted in [18, Lemma 1]: Lemma 2.2.
For any ϕ ∈ B N in dimension d , we have (cid:107)I N ϕ (cid:107) H k ≤ ( √ d (cid:107) ϕ (cid:107) H k , k ≥ . (2.5)Given I N f , the discrete spatial partial derivatives can be defined as:( D x f ) i,j = N (cid:88) k,l = − N +1 kπ i L ( ˆ f c ) k,l exp (cid:18) π i L ( kx i + ly j ) (cid:19) , ( D x f ) i,j = N (cid:88) k,l = − N +1 − k π L ( ˆ f c ) k,l exp (cid:18) π i L ( kx i + ly j ) (cid:19) . Similarly we have D y and D y . For any f, g ∈ M N , and f = ( f , f ) T , g = ( g , g ) T ∈M N × M N , we can define the discrete gradient, divergence and Laplace operators: ∇ N f = (cid:18) D x fD y f (cid:19) , ∇ N · f = D x f + D y f , ∆ N f = D x f + D y f. ETDMs2 for no-slope-selection thin film model L (denoted as (cid:96) ) inner product ( · , · ) N and norm (cid:107) · (cid:107) N :( f, g ) N = h N (cid:88) i,j =1 f ij g ij , (cid:107) f (cid:107) N = (cid:112) ( f, f ) N , ( f , g ) N = h N (cid:88) i,j =1 ( f ij g ij + f ij g ij ) , (cid:107) g (cid:107) N = (cid:112) ( f , f ) N . Similarly, we can define the discrete Sobolev norm (cid:107) · (cid:107) H h and the discrete Sobolevsemi-norm | · | H h : (cid:107) f (cid:107) H h = (cid:16) (cid:88) | α |≤ (cid:107) D α f (cid:107) N (cid:17) , | f | H h = (cid:16) (cid:88) | α | =2 (cid:107) D α f (cid:107) N (cid:17) , where as above α = ( α , α ) is a 2-tuple of nonnegative integers with | α | = α + α ,and D α = D α x D α y . Furthermore, the following summation by parts formulas areavailable in [21, proposition 2.2]:
Lemma 2.3.
For any f, g ∈ M N , and f = ( f , f ) T , g = ( g , g ) T ∈ M N ×M N ,we have the following summation by parts formula: ( f, ∇ N · g ) N = − ( ∇ N f, g ) N , ( f, ∆ N g ) N = − ( ∇ N f, ∇ N g ) N = (∆ N f, g ) N . Recall that the exact solution u in (1.2) is assumed to be mean value free, thuswe consider the subset of zero-mean grid functions: M N = { v ∈ M N | ( v, N = 0 } = { v ∈ M N | ˆ v = 0 } . Define L N = ε ∆ N | M N , which is symmetric positive definite on M N .In turn, the spatial discretization of (1.2) becomes: Given u ∈ M N , find ˆ u :[0 , T ] → M N such thatdˆ u d t = − L N ˆ u − f N (ˆ u ) , L N = ε ∆ N , t ∈ (0 , T ] , ˆ u (0) = u (0) , (2.6)where f N (ˆ u ) = ∇ N · (cid:16) ∇ N ˆ u |∇ N ˆ u | (cid:17) . Multiplying both sides of (2.6) by e L N t yieldsd e L N t ˆ u d t = − e L N t f N (ˆ u ) . (2.7)Integrating (2.7) from t n to t n +1 givesˆ u ( t n +1 ) = e − L N τ ˆ u ( t n ) − (cid:90) t n +1 t n e − L N ( t n +1 − t ) f N (ˆ u ( t ))d t. (2.8)By [21, Lemma 4.1], the error between the exact solution u ( t ) and the solution ˆ u ( t )of (2.8) is of order O ( N − m ), given that u e ( t ) is smooth enough. W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
Two exponential time differencingschemes (ETD1 and ETDMs2) have been proposed in [21], using the energy convexsplitting method. Ju et al used a similar spatial discretization as (2.8) is derived, withmodified L N and f N (ˆ u ( t )): L c = L N − κ ∆ N , κ > , (2.9) f e (ˆ u ( t )) = f N (ˆ u ( t )) + κ ∆ N ˆ u ( t ) . (2.10)For the ETD1, the term f e (ˆ u ( t )) in [ t n , t n +1 ] is simply approximated by f e (ˆ u ( t n )).For the ETDMs2, f e (ˆ u ( t )) is approximated by the linear Lagrange interpolation: f e (ˆ u ( t n )) + t − t n τ [ f e (ˆ u ( t n )) − f e (ˆ u ( t n − ))] , t ∈ [ t n , t n +1 ] . Let u h ( t ) be the numerical solution of ETD1 and ETDMs2. Denote u h ( t n ) as u nh for n ≥
0. Integrating from t n to t n +1 , the following expressions are obtained in [21]: • ETD1: u n +1 h = e − τL c u nh − φ ( L c ) f e ( u nh ) . • ETDMs2: u n +1 h = e − τL c u nh − φ ( L c ) f e ( u nh ) − φ ( L c )[ f e ( u nh ) − f e ( u n − h )] , in which φ ( x ) = x − ( I − e − xτ ) , φ ( x ) = x − [1 − ( xτ ) − ( I − e − xτ )] . The convergence analysis for both ETD1 and ETDMs2 was given by [21]. However,a theoretical proof of the long time energy stability of ETDMs2 was not provided.
In this section we propose a second order in timestabilized exponential time differencing multistep scheme (sETDMs2). In order toguarantee the long time energy stability, a stabilizing term Aτ d ∆ N ˆ u ( t ) d t is introduced,in which A is a positive constant. Also, we apply the same Lagrange approximationfor f N as in [21]. Define a continuous in time function for s ∈ [0 , τ ]: f N ,L ( s, ˆ u ( t n ) , ˆ u ( t n − )) := f N (ˆ u ( t n )) + sτ [ f N (ˆ u ( t n )) − f N (ˆ u ( t n − ))] . We present the following sETDMs2 scheme: • sETDMs2: For n ≥
1, find u n +1 s : [0 , t n +1 ] → M N such that for any t ∈ ( t n , t n +1 ],d u n +1 s ( t )d t + Aτ d∆ N u n +1 s ( t )d t = − L N u n +1 s ( t ) − f N ,L ( t − t n , u ns , u n − s ) , (2.11)with u n +1 s ( t ) := u ns ( t ) for t ∈ [0 , t n ].The initial step is obtained as follow: Find u s : [0 , t ] → M N such thatd u s ( t )d t + Aτ d∆ N u s ( t )d t = − L N u s ( t ) − f N ( u s ) , t ∈ (0 , t ] , (2.12)where u s = u (0). In fact, as for the initial step, the stabilizing term can also bereplaced by − Aτ d ∆ N u s ( t ) d t . Integrating (2.11) and (2.12) from t n to t n +1 , one obtainsexplicit expressions of u s and u n +1 s for sETDMs2: ETDMs2 for no-slope-selection thin film model • sETDMs2 (matrix form): For n ≥ u n +1 s = e −K N τ u ns − φ ( K N ) G N ( u ns ) − φ ( K N )[ G N ( u ns ) − G N ( u n − s )] , (2.13)and the intial step is computed by u s = e −K N τ u s − φ ( K N ) G N ( u s ) , (2.14)in which K N = ε ( I + Aτ ∆ N ) − ∆ N , G N ( v ) = ( I + Aτ ∆ N ) − f N ( v ) ,φ ( K N ) = K − N ( I − e −K N τ ) , φ ( K N ) = K − N [ I − ( K N τ ) − ( I − e −K N τ )] . Note that sETDMs2 and ETDMs2 have similar matrix forms for solutions at { t i } N t i =1 ,and they share similar computational complexity.
3. Long time energy stability of the sETDMs2 scheme.
In later analy-sis we shall make frequent use of the Gagliardo-Nirenberg inequality and Agmon’sinequality (see [1, 2, 34]): for function v in R d , | v | W j,p (Ω) ≤ C (cid:107) v (cid:107) − θL q (Ω) (cid:107) v (cid:107) θW m,r (Ω) , p = jd + θ (cid:18) r − md (cid:19) + (1 − θ ) 1 q , (3.1) (cid:107) v (cid:107) L ∞ (Ω) ≤ C (cid:107) v (cid:107) − n m L (Ω) (cid:107) v (cid:107) n m H m (Ω) , m > d , v ∈ H m (Ω) , (3.2)with 1 ≤ q, r ≤ ∞ , ≤ j < m, jm ≤ θ ≤ E ( u n +1 s ( t )) = E N ( u n +1 s ( t )) + √ − (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + ( √ τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) , (3.3)where E N ( v ) is the discrete version of the continuous energy functional E ( v ) in (1.1): E N ( v ) = (cid:18) −
12 ln(1 + |∇ N v | ) , (cid:19) N + ε (cid:107) ∆ N v (cid:107) N , ∀ v ∈ M N . (3.4)The main theoretical result of this section is the following theorem. Theorem 3.1.
Assume that A ≥ √ . The numerical system (2.11) is energystable in the sense that for any ≤ n ≤ N t − , ˜ E ( u n +1 s ) + (cid:32) A − √ (cid:33) τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) ≤ ˜ E ( u ns ) . Proof . We define β ( v ) := v | v | and denote ˆf N ,L ( t − t n , u ns , u n − s ) := β ( ∇ N u ns ) + t − t n τ [ β ( ∇ N u ns ) − β ( ∇ N u n − s )] . By definition we have: For any v, w ∈ M N and ζ ≥ f N ( v ) = ∇ N · β ( ∇ N v ) , f N ,L ( ζ, v, w ) = ∇ N · ˆf N ,L ( ζ, v, w ) . W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
Taking an inner product with (2.11) by d u n +1 s ( t ) d t gives that: For n ≥ (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N + Aτ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N + ε (cid:107) ∆ N u n +1 s ( t ) (cid:107) N d t = (cid:18) ˆf N ,L ( t − t n , u ns , u n − s ) , d ∇ N u n +1 s ( t )d t (cid:19) N . (3.5)Recall that E N ( u ) = (cid:0) − ln(1 + |∇ N u | ) , (cid:1) N + ε (cid:107) ∆ N u (cid:107) N , thusdd t E N ( u n +1 s ( t )) = (cid:18) − β ( ∇ N u n +1 s ( t )) , d ∇ N u n +1 s ( t )d t (cid:19) N + ε (cid:107) ∆ N u n +1 s ( t ) (cid:107) N d t . Therefore, we can rewrite (3.5) as below: (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N + Aτ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N + dd t E N ( u n +1 s ( t ))= (cid:18) ˆf N ,L ( t − t n , u ns , u n − s ) − β ( ∇ N u n +1 s ( t )) , d ∇ N u n +1 s ( t )d t (cid:19) N . (3.6)Integrating (3.6) from t n to t n +1 , we have (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + Aτ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + E N ( u n +1 s ) − E N ( u ns )= (cid:90) t n +1 t n (cid:18) β ( ∇ N u ns ) − β ( ∇ N u n +1 s ( t )) , d ∇ N u n +1 s ( t )d t (cid:19) N dt + (cid:90) t n +1 t n t − t n τ (cid:18) β ( ∇ N u ns ) − β ( ∇ N u n − s ) , d ∇ N u n +1 s ( t )d t (cid:19) N dt := ( I ) + ( II ) . (3.7)As for ( I ), an application of H¨older’s inequality gives( I ) ≤ (cid:90) t n +1 t n (cid:107) β ( ∇ N u n +1 s ( t )) − β ( ∇ N u ns ) (cid:107) N (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N dt. (3.8)Using the inequality | β ( v ) − β ( w ) | ≤ | v − w | in [21, Lemma 3.5] and reversing theorder of integration and summation, we obtain (cid:107) β ( ∇ N u n +1 s ( t )) − β ( ∇ N u ns ) (cid:107) N ≤ (cid:107)∇ N u n +1 s ( t ) − ∇ N u ns (cid:107) N = (cid:13)(cid:13)(cid:13)(cid:13)(cid:90) tt n d ∇ N u n +1 s ( l )d t dl (cid:13)(cid:13)(cid:13)(cid:13) N ≤ τ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) . (3.9)Substituting (3.9) into (3.8) and applying H¨older’s inequality again, we get( I ) ≤ τ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) . (3.10) ETDMs2 for no-slope-selection thin film model (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N ≤ τ λ (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N + τ λ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N , where λ is a positive constant to be decided later. Thus we arrive:( I ) ≤ λ (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + τ λ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) . (3.11)The term ( II ) could be estimated in a similar way:( II ) ≤ (cid:90) t n +1 t n t − t n τ (cid:107)∇ N u ns − ∇ N u n − s (cid:107) N (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N dt ≤ τ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) (cid:90) t n +1 t n t − t n τ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) N dt ≤ τ √ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) ≤ τ √ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + τ √ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) , (3.12)where we have applied the Cauchy-Schwarz inequality to obtain the last inequality.Again, we make use of the interpolation inequality and get( II ) ≤ √ λ (cid:13)(cid:13)(cid:13)(cid:13) d u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + τ λ √ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + 14 √ λ (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + τ λ √ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) . (3.13)Substituting (3.11) and (3.13) into (3.7), we arrive at0 ≥ (cid:32) − √ √ λ (cid:33) (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) − √ λ (cid:13)(cid:13)(cid:13)(cid:13) d u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + (cid:32) A − (2 √ λ √ (cid:33) τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) − λτ √ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + E N ( u n +1 s ) − E N ( u ns ) , (3.14)from which we can obtain the following restraints on λ and A :1 − √ √ λ ≥ √ λ ⇒ λ ≥ √ , (3.15) A − (2 √ λ √ ≥ λ √ ⇒ A ≥ √ λ ≥ (cid:32) √ (cid:33) = 2 + √ . (3.16)0 W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
Setting λ = √ , (3.14) is equal to − (cid:32) A − √ (cid:33) τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) ≥ √ − (cid:13)(cid:13)(cid:13)(cid:13) d u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) − √ − (cid:13)(cid:13)(cid:13)(cid:13) d u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + ( √ τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) − ( √ τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + E N ( u n +1 s ) − E N ( u ns ) ≥ ˜ E ( u n +1 s ) − ˜ E ( u ns ) . (3.17)Thus we have proved the energy stability of (2.11). Remark 3.2.
With a slight modification of the proof, we can obtain the energydecay property for the following numerical energy functional: ˆ E ( u n +1 s ) := E N ( u n +1 s ) + τ √ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) . (3.18) Applying the interpolation inequality to the LHS of (3.7) , substituting (3.10) and (3.12) into the RHS of (3.7) , we see that √ Aτ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + E N ( u n +1 s ) − E N ( u ns ) ≤ (6 + √ τ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + τ √ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) . Rearranging the terms, we also have (cid:32) √ A − √ (cid:33) τ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) + E N ( u n +1 s ) ≤ τ √ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u ns ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n − ,t n ; (cid:96) ) + E N ( u ns ) = ˆ E ( u ns ) . Therefore, ˆ E ( u ns ) + 2 (cid:32) √ A − √ (cid:33) τ (cid:13)(cid:13)(cid:13)(cid:13) d ∇ N u n +1 s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ( t n ,t n +1 ; (cid:96) ) ≤ ˆ E ( u ns ) , so that the energy estimate ˆ E ( u n +1 s ) ≤ ˆ E ( u ns ) is valid under the same assumption of A as in Theorem 3.1, if we notice that A ≥ (cid:16) √ (cid:17) = √ . Given the above energy stability and the fact that u n +1 s ( t ) ∈ M N , now we areable to derive a uniform in time bound of (cid:107) u n +1 s ( t ) (cid:107) H h . Proposition 3.3.
Assume that the initial solution u (0) has H h (Ω) -regularityand A ≥ √ . Then we have global in time bounds for the numerical solution E N ( u n +1 s ) ≤ E N ( u s ) , (cid:107) u n +1 s ( t ) (cid:107) H h ≤ C , for ≤ n ≤ N t − , (3.19) ETDMs2 for no-slope-selection thin film model where C only depends on ε , Ω and (cid:107) u (0) (cid:107) H h .Proof . Theorem 3.1 implies that E N ( u n +1 s ) ≤ ˜ E ( u n +1 s ) ≤ ˜ E ( u s ) . (3.20)Taking an inner product with d u s ( t )d t on both sides of (2.12) and performing a similaranalysis as in (3.5)-(3.11) yields that (cid:13)(cid:13)(cid:13)(cid:13) d u s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L (0 ,τ ; (cid:96) ) + Aτ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L (0 ,τ ; (cid:96) ) + E N ( u s ) − E N ( u s ) ≤ (cid:13)(cid:13)(cid:13)(cid:13) d u s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L (0 ,τ ; (cid:96) ) + τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L (0 ,τ ; (cid:96) ) . (3.21)In addition, the fact A ≥ √ > indicates that˜ E ( u s ) ≤ (cid:13)(cid:13)(cid:13)(cid:13) d u s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L (0 ,τ ; (cid:96) ) + 2 A − τ (cid:13)(cid:13)(cid:13)(cid:13) d∆ N u s ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L (0 ,τ ; (cid:96) ) + E N ( u s ) ≤ E N ( u s ) . (3.22)Therefore, we have E N ( u n +1 s ) ≤ E N ( u s ) for any 0 ≤ n ≤ N t −
1. Now we can applyLemma 3.2 and Remark 3.3 in [9, p. 586] to derive the desired conclusion.Next we provide a finite time uniform H h bound for the numerical solution u n +1 s ( t ); this subtle estimate will be used in the optimal rate convergence analysis. Proposition 3.4.
Assume that the initial solution u (0) has H h -regularity and A ≥ √ . Then we have the finite time H h bound for the numerical solution (cid:107) u n +1 s ( t ) (cid:107) H h ≤ C , for ≤ n ≤ N t − , (3.23) where C only depends on ε , Ω and (cid:107) u (cid:107) H h .Proof . Taking an inner product with − ∆ N u n +1 s ( t ) on both sides of (2.11) leadsto 12 d (cid:107)∇ N ∆ N u n +1 s ( t ) (cid:107) N d t + Aτ (cid:107)∇ N ∆ N u n +1 s ( t ) (cid:107) N d t + ε (cid:107)∇ N ∆ N u n +1 s ( t ) (cid:107) N = − t − t n − τ (cid:16) ∇ N f N ( u ns ) , ∇ N ∆ N u n +1 s ( t ) (cid:17) N + t − t n τ (cid:16) ∇ N f N ( u n − s ) , ∇ N ∆ N u n +1 s ( t ) (cid:17) N . (3.24)For any v ∈ H h (Ω) with periodic boundary conditions, recall that ˜ v = I N v is thecontinuous extension of v . By (2.4) we have (cid:107)∇ N f N ( v ) (cid:107) N = (cid:107)∇ N ∇ N · β ( v ) (cid:107) N = (cid:107)∇∇ · I N ( β ( ∇ ˜ v )) (cid:107) ≤ C (cid:107) β ( ∇ ˜ v ) (cid:107) H . (3.25)With an application of elliptic regularity, we get (cid:107) β ( ∇ ˜ v ) (cid:107) H ≤ C ( (cid:107) ∆ β ( ∇ ˜ v ) (cid:107) + (cid:107) β ( ∇ ˜ v ) (cid:107) ) , for ˜ v = ˜ u ns , ˜ u n +1 s . (3.26)2 W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
We also notice that (cid:107) β ( ∇ ˜ v ) (cid:107) ≤ (cid:107)∇ ˜ v (cid:107) ≤ C (cid:107) ∆˜ v (cid:107) , and | ∆ β ( ∇ ˜ v ) | ≤ |∇ ∆˜ v | |∇ ˜ v | + 6 | ∆˜ v | |∇ ˜ v | (1 + |∇ ˜ v | ) + 2 |∇ ˜ v | |∇ ∆˜ v | (1 + |∇ ˜ v | ) + 8 |∇ ˜ v | | ∆˜ v | (1 + |∇ ˜ v | ) ≤ |∇ ∆˜ v | + 14 | ∆˜ v | . (3.27)Also, since ∆˜ v has zero mean because of the periodic boundary conditions, thePoincar´e inequality (cid:107) ∆˜ v (cid:107) ≤ C (cid:107)∇ ∆˜ v (cid:107) is available. Combining this estimate with(3.1), we obtain (cid:107) β ( ∇ ˜ v ) (cid:107) H ≤ C ( (cid:107)∇ ∆˜ v (cid:107) + (cid:107) ∆˜ v (cid:107) L ) ≤ C (cid:107)∇ ∆˜ v (cid:107) + C (cid:107) ∆˜ v (cid:107)(cid:107) ∆˜ v (cid:107) H ≤ C (cid:107)∇ ∆˜ v (cid:107) , (3.28)in which we have used the H bound of ˜ v = ˜ u ns , ˜ u n +1 s . Substituting the aboveestimates into (3.24) and applying Lemma 3.3, we see that12 d (cid:107)∇ N ∆ N u n +1 s ( t ) (cid:107) N d t + Aτ (cid:107)∇ N ∆ N u n +1 s ( t ) (cid:107) N d t + ε (cid:107)∇ N ∆ N u n +1 s ( t ) (cid:107) N ≤ ε (cid:107)∇ N ∆ N u n +1 s ( t ) (cid:107) N + Cε − ( (cid:107)∇ N ∆ N u ns (cid:107) N + (cid:107)∇ N ∆ N u n − s (cid:107) N ) . (3.29)In turn, an integration from t n to t n +1 implies that12 (cid:0) (cid:107)∇ N ∆ N u n +1 s (cid:107) N − (cid:107)∇ N ∆ N u ns (cid:107) N (cid:1) + Aτ (cid:0) (cid:107)∇ N ∆ N u n +1 s (cid:107) N − (cid:107)∇ N ∆ N u ns (cid:107) N (cid:1) ≤ Cε − τ ( (cid:107)∇ N ∆ N u ns (cid:107) N + (cid:107)∇ N ∆ N u n − s (cid:107) N ) . (3.30)A summation from n = 1 to n shows that12 (cid:107)∇ N ∆ N u n +1 s (cid:107) N + Aτ (cid:107)∇ N ∆ N u n +1 s (cid:107) N ≤ Cτ ε − n (cid:88) i =0 (cid:107)∇ N ∆ N u is (cid:107) N + 12 (cid:107)∇ N ∆ N u s (cid:107) N + Aτ (cid:107)∇ N ∆ N u s (cid:107) N . An application pf the discrete Gronwall inequality yields (cid:107)∇ N ∆ N u n +1 s (cid:107) N ≤ C, (3.31)where C depends only on ε, u (0) and T .
4. Error analysis of the sETDMs2 scheme.
The goal of this section is toprovide an optimal rate convergence analysis for the sETDMs2 scheme (2.13). Recallthat h = L N , τ = TN t , β ( v ) = v | v | .The following inverse estimate is needed in the error analysis. Lemma 4.1.
For any u ∈ B N and ≤ p ≤ ∞ , | u | W m,p (Ω) ≤ Ch − k − ( − p ) d (cid:107) u (cid:107) H m − k (Ω) , ≤ k ≤ m. (4.1) ETDMs2 for no-slope-selection thin film model Proof . By (3.1) and (3.2) we have (cid:107) ∂ mx u (cid:107) L p (Ω) ≤ (cid:107) ∂ mx u (cid:107) − θL (Ω) (cid:107) ∂ mx u (cid:107) θH k (Ω) , θ = (cid:18) − p (cid:19) dk , where the integer k ≥ ( − p ) d for 2 ≤ p < ∞ and k > d for p = ∞ . Hence, | u | W m,p (Ω) ≤ | u | − θH m (Ω) (cid:107) u (cid:107) θH m + k (Ω) . (4.2)Meanwhile, an application of (2.1) (in Lemma 2.1) implies that | u | H m + k (Ω) ≤ Ch − k | u | H m (Ω) , k ≥ . (4.3)Combining (4.2) and (4.3) leads to | u | W m,p (Ω) ≤ Ch − ( − p ) d (cid:107) u (cid:107) H m (Ω) . (4.4)For k ≥
1, we simply apply (2.1) again: | u | W m,p (Ω) ≤ Ch − k − ( − p ) d (cid:107) u (cid:107) H m − k (Ω) . (4.5)We get the conclusion.The following estimate is needed in the analysis of the nonlinear term. Lemma 4.2.
For any v, w ∈ M N ∩ H h , g ∈ M N ∩ H h , we have (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N g |∇ N v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N ≤ C (1 + h )( (cid:107) ∆ N g (cid:107) N + (cid:107) g (cid:107) N ) , (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N w · ∇ N g ( ∇ N w + ∇ N v )(1 + |∇ N v | )(1 + |∇ N w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N ≤ C (1 + h )( (cid:107) ∆ N g (cid:107) N + (cid:107) g (cid:107) N ) , in which C depends only on Ω , (cid:107) v (cid:107) H h and (cid:107) w (cid:107) H h .Proof . Define ˜ v = I N v, ˜ w = I N w, ˜ g = I N g as the continuous extension of v, w, g , respectively. Thanks to (2.4), we have (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N g |∇ N v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N = (cid:13)(cid:13)(cid:13)(cid:13) ∇ · I N (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ (cid:13)(cid:13)(cid:13)(cid:13) ∇ · (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) + Ch (cid:13)(cid:13)(cid:13)(cid:13) ∇ ˜ g |∇ ˜ v | (cid:13)(cid:13)(cid:13)(cid:13) H . (4.6)For the first term on the RHS, using the Gagliardo-Nirenberg inequality (3.1), theinterpolation inequality, the H bound of (cid:107) ˜ v (cid:107) and the Young inequality, we get (cid:13)(cid:13)(cid:13)(cid:13) ∇ · (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ C ( | ˜ g | H + | ˜ g | , | ˜ v | , ) ≤ C | ˜ g | H + C (cid:107) ˜ g (cid:107) H (cid:107) ˜ g (cid:107) H (cid:107) ˜ v (cid:107) H ≤ C | ˜ g | H + C (cid:107) ˜ g (cid:107) H ≤ C (cid:107) ˜ g (cid:107) + C (cid:107) ∆˜ g (cid:107) . (4.7)For the second term, the elliptic regularity and the interpolation inequality shows that (cid:13)(cid:13)(cid:13)(cid:13) ∇ ˜ g |∇ ˜ v | (cid:13)(cid:13)(cid:13)(cid:13) H ≤ C (cid:18)(cid:13)(cid:13)(cid:13)(cid:13) ∆ (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13) ∇ ˜ g |∇ ˜ v | (cid:13)(cid:13)(cid:13)(cid:13)(cid:19) ≤ C (cid:18)(cid:13)(cid:13)(cid:13)(cid:13) ∆ (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) + (cid:107) ˜ g (cid:107) + (cid:107) ∆˜ g (cid:107) (cid:19) . W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
On the other hand, we see that (cid:12)(cid:12)(cid:12)(cid:12) ∆ (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ≤ C ( |∇ ∆˜ g | + |∇∇ ˜ g ||∇∇ ˜ v | + |∇ ˜ g ||∇∇ ˜ v | + |∇ ˜ g || ∆ ∇ ˜ v | ) . (4.8)Then we obtain (cid:13)(cid:13)(cid:13)(cid:13) ∆ (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ C ( | ˜ g | H + | ˜ g | , | ˜ v | , + | ˜ g | , | ˜ v | , + | ˜ g | , ∞ | ˜ v | H ) . (4.9)Applying Lemma 2.1, Lemma 4.1 and the Sobolev inequality, we get (cid:13)(cid:13)(cid:13)(cid:13) ∆ (cid:18) ∇ ˜ g |∇ ˜ v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ Ch − | ˜ g | H + Ch − (cid:107) ˜ g (cid:107) H (cid:107) ˜ v (cid:107) H + C (cid:107) ˜ g (cid:107) H (cid:107) ˜ v (cid:107) H + Ch − | ˜ g | L ∞ | ˜ v | H ≤ Ch − | ˜ g | H + Ch − (cid:107) ˜ g (cid:107) H (cid:107) ˜ v (cid:107) H + C (cid:107) ˜ g (cid:107) H (cid:107) ˜ v (cid:107) H + Ch − (cid:107) ˜ g (cid:107) H (cid:107) ˜ v (cid:107) H ≤ C ( h − + 1) (cid:107) ˜ g (cid:107) H . (4.10)Substituting the above estimate into (4.6), we arrive at (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N g |∇ N v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N ≤ C (cid:107) ∆˜ g (cid:107) + C (cid:107) ˜ g (cid:107) + C (1 + h ) (cid:107) ˜ g (cid:107) H ≤ C (1 + h ) (cid:107) ∆˜ g (cid:107) + C (1 + h ) (cid:107) ˜ g (cid:107) = C (1 + h )( (cid:107) ∆ N g (cid:107) N + (cid:107) g (cid:107) N ) . (4.11)Similar estimates can be applied on the second part of the conclusion: (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N w · ∇ N g ( ∇ N w + ∇ N v )(1 + |∇ N v | )(1 + |∇ N w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N = (cid:13)(cid:13)(cid:13)(cid:13) ∇ · I N (cid:18) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ (cid:13)(cid:13)(cid:13)(cid:13) ∇ · (cid:18) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) + Ch (cid:13)(cid:13)(cid:13)(cid:13) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:13)(cid:13)(cid:13)(cid:13) H . (4.12)Analyzing the first term on the RHS in a similar fashion as in (4.7), we have (cid:13)(cid:13)(cid:13)(cid:13) ∇ · (cid:18) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ C ( | ˜ w | , | ˜ g | , + | ˜ g | H + | ˜ g | , ( | ˜ w | , + | ˜ v | , )) ≤ C (cid:107) ∆˜ g (cid:107) + C (cid:107) ˜ g (cid:107) . (4.13)The second term can also be analyzed by the elliptic regularity and the interpolationinequality: (cid:13)(cid:13)(cid:13)(cid:13) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:13)(cid:13)(cid:13)(cid:13) H ≤ C (cid:13)(cid:13)(cid:13)(cid:13) ∆ (cid:18) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) + C (cid:13)(cid:13)(cid:13)(cid:13) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:13)(cid:13)(cid:13)(cid:13) ≤ C (cid:13)(cid:13)(cid:13)(cid:13) ∆ (cid:18) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) + C ( (cid:107) ˜ g (cid:107) + | ˜ g | H ) . ETDMs2 for no-slope-selection thin film model (cid:13)(cid:13)(cid:13)(cid:13) ∆ (cid:18) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) ≤ C | ˜ w | H | ˜ g | , ∞ + C | ˜ g | H + C | ˜ g | , ( | ˜ w | , + C | ˜ v | , ) + C | ˜ g | , | ˜ w | , + C | ˜ g | , | ˜ v | , + C | ˜ g | , | ˜ w | , | ˜ v | , + C | ˜ g | , ( | ˜ w | , + | ˜ v | , ) . (4.14)Again we apply Lemma 2.1, (3.1), (3.2) and the Sobolev inequality as in (4.10): (cid:13)(cid:13)(cid:13)(cid:13) ∇ ˜ w · ∇ g ( ∇ ˜ w + ∇ ˜ v )(1 + |∇ ˜ v | )(1 + |∇ ˜ w | ) (cid:13)(cid:13)(cid:13)(cid:13) H ≤ C ( (cid:107) ˜ g (cid:107) + | ˜ g | H ) + Ch − (cid:107) ˜ w (cid:107) H (cid:107) ˜ g (cid:107) L ∞ + Ch − (cid:107) ˜ g (cid:107) H + Ch − (cid:107) ˜ g (cid:107) H ( (cid:107) ˜ w (cid:107) H + (cid:107) ˜ v (cid:107) H )+ C (cid:107) ˜ g (cid:107) H ( (cid:107) ˜ w (cid:107) H + (cid:107) ˜ v (cid:107) H + (cid:107) ˜ w (cid:107) H (cid:107) ˜ v (cid:107) H )+ Ch − (cid:107) ˜ g (cid:107) H ( (cid:107) ˜ w (cid:107) H + (cid:107) ˜ v (cid:107) H ) ≤ C (1 + h − ) (cid:107) ˜ g (cid:107) H ≤ C (1 + h − )( (cid:107) ˜ g (cid:107) + (cid:107) ∆˜ g (cid:107) ) . (4.15)Substituting (4.13) and (4.15) into (4.12) yields that (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N w · ∇ N g ( ∇ N w + ∇ N v )(1 + |∇ N v | )(1 + |∇ N w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N ≤ C (1 + h ) (cid:107) ∆˜ g (cid:107) + C (1 + h ) (cid:107) ˜ g (cid:107) = C (1 + h )( (cid:107) ∆ N g (cid:107) N + (cid:107) g (cid:107) N ) , (4.16)which is the desired conclusion.Now we present an estimate of the nonlinear term. Lemma 4.3.
For any v, w, M N ∩ H h and g ∈ M N ∩ H h , we have (cid:16) β ( ∇ N v ) − β ( ∇ N w ) , ∇ N g (cid:17) N ≤ C v,w (2 + ε )(1 + h )2 ε (cid:107) v − w (cid:107) N + ε (cid:107) ∆ N g (cid:107) N + 12 (cid:107) g (cid:107) N , where C v,w is a constant depending on Ω , (cid:107) w (cid:107) H h and (cid:107) v (cid:107) H h .Proof . We begin with the following identity: For any v, w ∈ M N , (cid:16) β ( ∇ N v ) − β ( ∇ N w ) , ∇ N g (cid:17) N = (cid:18) β ( ∇ N v ) − ∇ N w |∇ N v | , ∇ N g (cid:19) N + (cid:18) ∇ N w |∇ N v | − β ( ∇ N w ) , ∇ N g (cid:19) N = (cid:18) ∇ N ( v − w ) , ∇ N g |∇ N v | (cid:19) N + (cid:18) ∇ N w |∇ N w | − |∇ N v | (1 + |∇ N v | )(1 + |∇ N w | ) , ∇ N g (cid:19) N = (cid:18) ∇ N ( v − w ) , ∇ N g |∇ N v | (cid:19) N + (cid:18) ∇ N ( w − v ) , ∇ N w · ∇ N g ( ∇ N w + ∇ N v )(1 + |∇ N v | )(1 + |∇ N w | ) (cid:19) N . W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
An application of integration by parts gives (cid:18) ∇ N v |∇ N v | − ∇ N w |∇ N w | , ∇ N g (cid:19) N = − (cid:18) v − w, ∇ N · (cid:18) ∇ N g |∇ N v | (cid:19)(cid:19) N − (cid:18) w − v, ∇ N · (cid:18) ∇ N w · ∇ N g ( ∇ N w + ∇ N v )(1 + |∇ N v | )(1 + |∇ N w | ) (cid:19)(cid:19) N ≤ (cid:107) v − w (cid:107) N (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N g |∇ N v | (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N + (cid:107) v − w (cid:107) N (cid:13)(cid:13)(cid:13)(cid:13) ∇ N · (cid:18) ∇ N w · ∇ N g ( ∇ N w + ∇ N v )(1 + |∇ N v | )(1 + |∇ N w | ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13) N . (4.17)Applying Lemma 4.2, we have the following estimate: (cid:18) ∇ N v |∇ N v | − ∇ N w |∇ N w | , ∇ N g (cid:19) N ≤ C (1 + h ) (cid:107) v − w (cid:107) N ( (cid:107) ∆ N g (cid:107) N + (cid:107) g (cid:107) N ) ≤ ε (cid:107) ∆ N g (cid:107) N + 12 (cid:107) g (cid:107) N + C (1 + h ) (cid:18) ε + 12 (cid:19) (cid:107) v − w (cid:107) N , (4.18)where C is a constant dependent on (cid:107) w (cid:107) H h , (cid:107) v (cid:107) H h , Ω N .Below is the main result of this section. Theorem 4.4.
Assume that the exact solution satisfies u e ∈ H (0 , T ; H m +4 per (Ω)) ∩ H (0 , T ; H (Ω)) with m ≥ . Let u ( t ) = u e ( t ) | Ω N and u s = u (0) ∈ H h . Denote by { u ns } N t n =1 the solution of the sETDMs2 scheme (2.13) . If the time step satisfies τ ≤ , (4.19) then we have (cid:107) u ( t n ) − u ns (cid:107) N ≤ C ( τ + N − m ) , ≤ n ≤ N t , (4.20) in which C > is independent of τ and the spatial discretization parameter N .Proof . Define error function e ( t ) = u ( t ) − u n +1 s ( t ) ∈ M N . Below we denote u e ( t n ), u ( t n ) as u ne and u n , respectively. Subtracting (2.11) from (1.2) givesd e ( t )d t + Aτ d∆ N e ( t )d t + ε ∆ N e ( t )= − (cid:2) f N ,L ( t − t n , u n , u n − ) − f N ,L ( t − t n , u ns , u n − s ) (cid:3) + R ( t ) , (4.21)for t ≤ t n +1 , where R ( t ) is the truncation error: R ( t ) = − ε (∆ − ∆ N ) u ( t ) + Aτ d∆ N u ( t )d t + (cid:104) ∇ N · β ( ∇ N u ( t )) − ∇ · β ( ∇ u e ( t )) (cid:12)(cid:12)(cid:12) Ω N (cid:105) + (cid:16) t − t n τ ∇ N · (cid:2) β ( ∇ N u ( t )) − β ( ∇ N u n − ) (cid:3) − t − t n − τ ∇ N · [ β ( ∇ N u ( t )) − β ( ∇ N u n )] (cid:17) := R ( t ) + R ( t ) + R ( t ) + R ( t ) . (4.22) ETDMs2 for no-slope-selection thin film model e ( t ) on both sides of (4.22) yields12 d (cid:107) e ( t ) (cid:107) N d t + Aτ (cid:107) ∆ N e ( t ) (cid:107) N d t + ε (cid:107) ∆ N e ( t ) (cid:107) N = (cid:16) ˆf N ,L ( t − t n , u ( t n ) , u ( t n − )) − ˆf N ,L ( t − t n , u ns , u n − s ) , ∇ N e ( t ) (cid:17) N + ( R ( t ) , e ( t ) N )= (NL) + ( R ( t ) , e ( t )) N . (4.23)The truncation error term could be handled by the Cauchy-Schwarz inequality:( R ( t ) , e ( t )) N ≤ (cid:107)R ( t ) (cid:107) N (cid:107) e ( t ) (cid:107) N . (4.24)The nonlinear term on the RHS has the following form:(NL) = t + τ − t n τ (cid:16) β ( ∇ N u n ) − β ( ∇ N u ns ) , ∇ N e ( t ) (cid:17) N − t + t n τ (cid:16) β ( ∇ N u n − ) − β ( ∇ N u n − s ) , ∇ N e ( t ) (cid:17) N . (4.25)An application of Lemma 3.4 and Lemma 4.3 to (4.25) results in(NL) ≤ C (1 + h )(2 + ε )2 ε (cid:0) (cid:107) e ( t n ) (cid:107) N + (cid:107) e ( t n − ) (cid:107) N (cid:1) + 3 ε (cid:107) ∆ N e ( t ) (cid:107) N + 32 (cid:107) e ( t ) (cid:107) N , (4.26)where C depends on (cid:107) e ( t n ) (cid:107) H h , (cid:107) e ( t n − ) (cid:107) H h and Ω N , and thus is a constant depen-dent only of (cid:107) u (cid:107) H h and Ω N .Substituting (4.24) and (4.26) into (4.23), we obtain12 d (cid:107) e ( t ) (cid:107) N d t + Aτ (cid:107) ∆ N e ( t ) (cid:107) N d t + ε (cid:107) ∆ N e ( t ) (cid:107) N ≤ C (1 + h )(2 + ε )2 ε (cid:0) (cid:107) e ( t n ) (cid:107) N + (cid:107) e ( t n − ) (cid:107) N (cid:1) + 2 (cid:107) e ( t ) (cid:107) N + (cid:107)R ( t ) (cid:107) N . (4.27)Denote ω ( t ) = (cid:107) e ( t ) (cid:107) N + Aτ (cid:107) ∆ N e ( t ) (cid:107) N . Multiplying both sides by e − t givesdd t e − t ω ( t ) ≤ (cid:20) C (1 + h )(2 + ε )2 ε (cid:0) (cid:107) e ( t n ) (cid:107) N + (cid:107) e ( t n − ) (cid:107) N (cid:1) + (cid:107)R ( t ) (cid:107) N (cid:21) e − t . (4.28)Integrating (4.28) from t n to t n +1 and multiplying both sides by e t n , we have e − τ ω ( t n +1 ) − ω ( t n ) ≤ − e − τ C (1 + h )(2 + ε ) ε (cid:0) (cid:107) e ( t n ) (cid:107) N + (cid:107) e ( t n − ) (cid:107) N (cid:1) + (cid:107)R ( t ) (cid:107) L ( t n ,t n +1 ; (cid:96) ) . Since e x ≥ x for x ∈ R , from the above inequality we can further obtain that ω ( t n +1 ) − ω ( t n ) − τ ω ( t n +1 ) ≤ C (1 + h )(2 + ε ) τε (cid:0) (cid:107) e ( t n ) (cid:107) N + (cid:107) e ( t n − ) (cid:107) N (cid:1) + (cid:88) i =1 (cid:107) R i ( t ) (cid:107) L ( t n ,t n +1 ; (cid:96) ) . (4.29)8 W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
Note that R i ( t ) = (cid:101) R i ( t ) (cid:12)(cid:12)(cid:12) Ω N , for 1 ≤ i ≤
3, where (cid:107) (cid:101) R ( t ) (cid:107) L = (cid:107) − ε (cid:0) ∆ u e ( t ) − ∆ I N u ( t ) (cid:1) (cid:107) L ≤ Cε h m (cid:107) u e ( t ) (cid:107) H m +4 , (4.30) (cid:107) (cid:101) R ( t ) (cid:107) L = (cid:13)(cid:13)(cid:13)(cid:13) Aτ d∆ I N u e ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L = Aτ (cid:13)(cid:13)(cid:13)(cid:13) ∆ I N d u e ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) L ≤ Cτ h m (cid:13)(cid:13)(cid:13)(cid:13) d u e ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) H m +4 + Aτ (cid:13)(cid:13)(cid:13)(cid:13) d u e ( t )d t (cid:13)(cid:13)(cid:13)(cid:13) H , (4.31) (cid:107) (cid:101) R ( t ) (cid:107) L = (cid:107)∇ · [ β ( ∇ N I N u e ( t )) − β ( ∇ u e ( t ))] (cid:107) L ≤ C (cid:107)I N u e ( t ) − u e ( t ) (cid:107) H ≤ C (cid:107) u e ( t ) (cid:107) H m +2 h m , (4.32)where we have applied the approximation property of I N in (2.4). As for R ( t ), let α ( t ) := ∇ N · β ( ∇ N u ( t )), so that R ( t ) = − α ( t ) + t − t n − τ α ( t n ) − t − t n τ α ( t n − )= (cid:90) t n t ( t − t n − )( t n − s ) τ α (cid:48)(cid:48) ( s ) d s − (cid:90) t n − t ( t − t n )( t n − − s ) τ α (cid:48)(cid:48) ( s ) d s. (4.33)Applying the Cauchy-Schwarz inequality and the estimate (2.4) gives that (cid:107)R ( t ) (cid:107) N ≤ Cτ (cid:107)∇ N · β ( ∇ N u ( t )) (cid:107) H ( t n − ,t n +1 ; L h ) ≤ Cτ (cid:107) u e ( t ) (cid:107) H ( t n − ,t n +1 ; H ) . (4.34)Then we arrive at (cid:88) i =1 (cid:107) R i ( t ) (cid:107) L ( t n ,t n +1 ; L h ) = (cid:88) i =1 (cid:107) (cid:101) R i ( t ) (cid:107) L ( t n ,t n +1 ; L ) + (cid:107) R ( t ) (cid:107) L ( t n ,t n +1 ; L h ) ≤ C ( h m + τ )( (cid:107) u e ( t ) (cid:107) H ( t n ,t n +1 ; H m +4 ) + (cid:107) u e ( t ) (cid:107) H ( t n − ,t n +1 ; H ) ) . Substituting the above estimate into (4.29) yields ω ( t n +1 ) − ω ( t n ) − τ ω ( t n +1 ) ≤ C (1 + h )(2 + ε ) τε (cid:0) (cid:107) e ( t n ) (cid:107) N + (cid:107) e ( t n − ) (cid:107) N (cid:1) + C ( h m + τ )( (cid:107) u e ( t ) (cid:107) H ( t n ,t n +1 ; H m +4 ) + (cid:107) u e ( t ) (cid:107) H ( t n − ,t n +1 ; H ) ) . (4.35)A summation of the above inequality from 1 to n results in ω ( t n +1 ) − ω ( t ) − τ n +1 (cid:88) i =2 ω ( t i ) ≤ C (1 + h )(2 + ε )2 ε τ n (cid:88) i =0 (cid:107) e ( t i ) (cid:107) N + C ( h m + τ ) , where C is dependent of (cid:107) u e (cid:107) W , (0 ,T ; H m +4 ) , (cid:107) u e ( t ) (cid:107) H (0 ,T ; H ) . Since 4 τ ≤ , wehave12 ω ( t n +1 ) ≤ (cid:20) C (1 + h )(2 + ε ) ε + 4 (cid:21) τ n (cid:88) i =0 ω ( t i ) + C ( h m + τ ) + ω ( t ) . (4.36) ETDMs2 for no-slope-selection thin film model ≤ t ≤ τ ,d e ( t )d t + Aτ d∆ N e ( t )d t + ε ∆ N e ( t ) = −∇ N · (cid:2) β ( ∇ N u ) − β ( ∇ N u s ) (cid:3) + R ( t )= R ( t ) , (4.37)where we have applied the set-up u s = u , and R ( t ) iss the truncation error: R ( t ) = − ε (∆ − ∆ N ) u ( t ) + Aτ d∆ N u ( t )d t − [ ∇ · β ( ∇ u ( t )) − ∇ N · β ( ∇ u ( t ))] − ∇ N · [ β ( ∇ u ( t )) − β ( ∇ N u )] . (4.38)Similarly, we take an inner product with e ( t ) on both sides,12 d (cid:107) e ( t ) (cid:107) N d t + Aτ (cid:107) ∆ N e ( t ) (cid:107) N d t + ε (cid:107) ∆ N e ( t ) (cid:107) N ≤ (cid:107) e ( t ) (cid:107) N + C ( h m + τ ) . (4.39)Repeating the process from (4.28) to (4.36) for ω ( t ) in [0 , τ ], since ω (0) = 0, we get e − τ ω ( t ) ≤ Cτ ( h m + τ ) , which again leads to ω ( t ) ≤ Ce τ ( τ h m + τ ) . (4.40)A substitution of (4.40) into (4.36) implies that12 ω ( t n +1 ) ≤ (cid:20) C (1 + h )(2 + ε ) ε + 4 (cid:21) τ n (cid:88) i =0 ω ( t i ) + C ( h m + τ ) . (4.41)Applying the discrete Gronwall’s inequality and recall that h is bounded above, weget ω ( t n +1 ) ≤ C T,ε ( h m + τ ), i.e., (cid:107) e ( t n +1 ) (cid:107) N + Aτ (cid:107) ∆ N e ( t n +1 ) (cid:107) N ≤ C T,ε,u ( h m + τ ) , (4.42)where C T,ε,u is a constant dependent on (cid:107)∇ N ∆ N u (cid:107) N , T and ε . Remark 4.5.
Note that in Theorem 4.4, the parameter ε is not involved in the re-quirement of time step τ , in comparison with a previous work [31]. This improvementcomes from the treatment of the nonlinear term in Lemma 4.2, where we have boundedthe nonlinear term by the sum (instead of the product) of (cid:107) ∆ N e ( t ) (cid:107) and (cid:107) e ( t ) (cid:107) , andthus avoided an ε − coefficient of (cid:107) e ( t ) (cid:107) after applying Cauchy-Schwarz inequality toobtain an ε coefficient of (cid:107) ∆ N e ( t ) (cid:107) . In fact, the restriction can be further relaxed,since we only need τ < in (4.36) for the coefficient of ω ( t n +1 ) to be positive. Remark 4.6.
There have been some recent works for the linearized energy stablenumerical schemes for various 3-D gradient flows, such as a recent paper [29] for theCahn-Hilliard model. In fact, all the results presented in this article could be extendedto the 3-D equation, with careful treatments in the Sobolev analysis, although the NSSequation (1.2) is of physical interests only in the 2-D case. We may consider anassociated extension for other related gradient models in the future works. W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
5. Numerical results.
In this section, we first demonstrate various numericalexperiments to verify the temporal convergence and energy stability of the sETDMs2scheme (2.13). In particular, a comparison with the the ETDMs2 scheme (in [21])are performed, with κ = in ETDMs2, and A = , − , √ in sETDMs2. Andalso, some long time simulation results computed by the sETDMs2 are presented,and some long time characteristics of the numerical solution are studied, such as theenergy decay rate, the average surface roughness and average slope. Throughout this subsection, we consider solving(1.2) with Ω = [0 , π ] , ε = 0 . T = 1, and the initial value u e (0) = sin x cos y on the uniform N × N mesh with N = 256. The time step size is set to be τ =0 . ∗ [2 − , − , − , − , − , − ]. With an additional time-dependent forcing term g ( t ), the exact solution to equation (1.2) is given by u e ( t ) = sin x cos y cos t : g ( t ) = − sin x cos y sin t + 4 ε sin x cos y cos t + − x cos y cos t (cos t ) (1 + cos 2 x cos 2 y )+ (cos t ) (cos x cos y sin 2 x cos 2 y − sin x sin y cos 2 x sin 2 y )[1 + (cos t ) (1 + cos 2 x cos 2 y )] . To solve the molecular epitaxial growth equation with a forcing term, we appliedboth the sETDMs2 and the ETDMs2 schemes. For each algorithm, we computedthe relative error (cid:13)(cid:13) u N t − u N t τ (cid:13)(cid:13) N / (cid:107) u N t (cid:107) N between the exact solution and numericalsolutions. Results are displayed in Figure 5.1, from which the second order temporalconvergence is clear observed. Also, different choices of the stabilization coefficient A in the sETDMs2 scheme don’t have an evident impact on the reported error values. -5 -4 -3 -2 Time step: -6 -5 -4 -3 -2 -1 L e rr o r Discrete L2 relative error of sETDMs2 and ETDMs2
Fig. 5.1 . Temporal convergence rates of the sETDMs2, ETDMs2 schemes.
In this subsection, we set Ω = [0 , . , ε = 0 . T = 25000 and use a random initial data. The reason for such a choice comes from ETDMs2 for no-slope-selection thin film model
21a subtle fact that, a random initial data leads to a solution with various wave fre-quencies, so that more detailed and interesting phenomenon associated with differentwave frequencies are presented in the coarsening process, while a smooth initial valuemay yield a solution with trivial structure in a short time.We use a coarser uniform mesh with N = 128, and set time step size τ = 0 .
001 for t < τ = 0 .
01 for 200 ≤ t < τ = 0 .
02 for 1000 ≤ t < τ = 0 .
04 for t ≥ u at time t = 1, 1500, 5000, 15000, respectively.It can be observed that the solution has saturated to a one-hill-one-valley structureat the final time. Specially, when A in sETDMs2 is lower than the restriction inTheorem 3.1, the numerical solution still converges to the stable state. (k = 1/8) ETDMs2 at t = 1500 (k = 1/8) ETDMs2 at t = 5000 (k = 1/8) ETDMs2 at t = 15000 (A1) sETDMs2 at t = 1500 (A1) sETDMs2 at t = 5000 (A1) sETDMs2 at t = 15000 (A2) sETDMs2 at t = 1500 (A2) sETDMs2 at t = 5000 (A2) sETDMs2 at t = 15000 (A3) sETDMs2 at t = 1500 (A3) sETDMs2 at t = 5000 (A3) sETDMs2 at t = 15000 Fig. 5.2 . Snapshots of the numerical solutions. The first row is for ETDMs2 and the next 3rows are for sETDMs2 with A = , − , √ , respectively. Recall the discrete energy functional in (3.4), for convenience we repeat it here: E N ( u ) = (cid:18) −
12 ln(1 + |∇ N u | ) , (cid:19) N + ε (cid:107) ∆ N u (cid:107) N , ∀ u ∈ M N . (5.1)2 W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG
Also, consider the average surface roughness h N ( u ) and the average slope m N ( u ): h N ( u, t ) = (cid:115) h | Ω | (cid:88) M N | u ( x i,j , t ) − ¯ u ( t ) | , with ¯ u ( t ) := h | Ω | (cid:88) M N u ( x i,j , t ) . (5.2) m N ( u, t ) = (cid:115) h | Ω | (cid:88) M N |∇ u ( x i,j , t ) | . (5.3)For the no-slope-selection growth model (1.2), it has been proved that E N ∼ O ( − ln( t )), h N ∼ O ( t ) and m N ∼ O ( t ) as t → ∞ . (See [16, 27, 28] and references therein.) Theevolution of E N is demonstrated in Figure 5.3, while h N and m N in Figure 5.4. Thelinear fitting results are also presented in Figure 5.3 and Figure 5.4 for the solutionof sETDMs2 in time interval [1 , A have little impact onthe long time characteristics of the numerical solutions. Table 5.1
Fitting coefficients of long time characteristics for sETDMs2. E N ( t ) ≈ a ln( t ) + b h N ( t ) ≈ at b m N ( t ) ≈ at b a b a b a bA = 1 / A = 1 /
100 -39.38 -54.27 0.4377 0.4865 2.132 0.2614 A = √ -39.26 -54.77 0.4346 0.4879 2.137 0.2609 -1 Time -500-450-400-350-300-250-200-150-100-500 E ne r g y Energy evolution of sETDMs2 solution
Fig. 5.3 . Semi-log plot of the energy E N with ε = 0 . of sETDMs2 with A = , − , √ . Fitted line has the form a ln( t ) + b , with coefficients shown in Table 5.1. The numerical results displayed in Table 5.1 have demonstrated the robustnessof the energy stable numerical solvers in the long time scale simulations. With the
ETDMs2 for no-slope-selection thin film model -1 Time -2 -1 A v e r age he i gh t Evolution of average height -1 Time -1 M ound w i d t h g r o w t h r a t e Evolution of average mound width growth rate
Fig. 5.4 . The log-log plot of (1) the average surface roughness h N and (2) the average slope m N using sETDMs2 with ε = 0 . , A = , − , √ . Fitted lines have the form at b , withcoefficients shown in Table 5.1. time step size of order 10 − to 10 − , and such a small interface width parameter ε ,the numerical solver is able to capture the long time asymptotic growth rate of thesurface roughness and average slope, with relative accuracy within 4%. Therefore, asecond order accurate, energy stable numerical algorithm is a worthwhile approachfor a gradient flow with a long time scale coarsening process.
6. Concluding remarks.
In this paper we have presented a linear, second orderin time accurate, energy stable ETD-based scheme for a thin film model without slopeselection. An unconditional long time energy stability is justified at a theoretical level.In addition, we provide an O ( τ )-order convergence analysis in the (cid:96) ∞ (0 , T ; (cid:96) ) norm,when τ ≤ . Moreover, various numerical experiments have demonstrated that theproposed second-order scheme is able to produce accurate long time numerical resultswith a reasonable computational cost. Acknowledgment.
This work is supported in part by the grants NSFC 11671098,91630309, a 111 project B08018 (W. Chen), NSF DMS-1418689 (C. Wang), NSFDMS-1715504 and Fudan University start-up (X. Wang). C. Wang also thanks theKey Laboratory of Mathematics for Nonlinear Sciences, Fudan University, for supportduring his visit.
REFERENCES[1] R. A. Adams and J. J. F. Fournier,
Sobolev spaces , 2nd ed., Academic press, Singapore, 2003.[2] S. Agmon,
Lectures on elliptic boundary value problems , vol. 369, American Mathematical Soc.,2010.[3] B. Benesova, C. Melcher, and E. Suli,
An implicit midpoint spectral approximation of nonlocalCahn-Hilliard equations , SIAM J. Numer. Anal., (2014), 1466–1496.[4] G. Beylkin, J. M. Keiser, and L. Vozovoi, A new class of time discretization schemes for thesolution of nonlinear PDEs , J. Comput. Phys., (1998), 362–387.[5] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang,
Spectral methods: fundamentalsin single domains , Springer, 2006.[6] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang,
Spectral methods: evolution tocomplex geometries and applications to fluid dynamics , Springer Science & Business Media,2007.[7] C. Canuto and A. Quarteroni,
Approximation results for orthogonal polynomials in sobolevspaces , Math. Comput., (1982), no. 157, 67–86.[8] W. Chen, S. Conde, C. Wang, X. Wang, and S. M. Wise, A linear energy stable scheme for athin film model without slope selection , J. Sci. Comput., (2012), no. 3, 546–562.[9] W. Chen, C. Wang, X. Wang, and S. M. Wise, A linear iteration algorithm for a second-order W. CHEN, W. LI, Z. LUO, C. WANG AND X. WANG energy stable scheme for a thin film model without slope selection , J. Sci. Comput., (2014), no. 3, 574–601.[10] W. Chen and Y. Wang, A mixed finite element method for thin film epitaxy , Numer. Math., (2012), no. 4, 771–793.[11] K. Cheng, Z. Qiao, and C. Wang,
A third order exponential time differencing numeri- calscheme for no-slope-selection epitaxial thin film model with energy stability , J. Sci. Com-put., submitted and in review, 2018.[12] S. M. Cox and P. C. Matthews,
Exponential time differencing for stiff systems , J. Comput.Phys., (2002), 430–455.[13] G. Ehrlich and F. G. Hudda,
Atomic view of surface self-diffusion: Tungsten on tungsten , J.Chem. Phys., (1966), no. 3, 1039–1049.[14] D. J. Eyre, Unconditionally gradient stable time marching the Cahn-Hilliard equation , MRSOnline Proceedings Library Archive, (1998).[15] W. Feng, C. Wang, S. M. Wise, and Z. Zhang,
A second-order energy stable backward differ-entiation formula method for the epitaxial thin film equation with slope selection , arXivpreprint, arXiv:1706.01943 (2017).[16] L. Golubovi´c,
Interfacial coarsening in epitaxial growth models without slope selection , Phys.Rev. Lett., (1997), 90–93.[17] D. Gottlieb and S. A. Orszag, Numerical analysis of spectral methods: theory and applications ,vol. 26, SIAM, 1977.[18] S. Gottlieb and C. Wang,
Stability and convergence analysis of fully discrete fourier collocationspectral method for 3-D viscous burgers’ equation , J. Sci. Comput., (2012), no. 1, 102–128.[19] M. Hochbruck and A. Ostermann, Exponential integrators , Acta Numer., (2010), 209–286.[20] M. Hochbruck and A. Ostermann, Exponential multistep methods of Adams-type , BIT Numer.Math., (2011), 889–908.[21] L. Ju, X. Li, Z. Qiao, and H. Zhang, Energy stability and convergence of exponential timedifferencing schemes for the epitaxial growth model without slope selection , Math. Comput., (2018), no. 312, 1859–1885.[22] L. Ju, X. Liu, and W. Leng, Compact implicit integration factor methods for a family ofsemilinear fourth-order parabolic equations , Discrete Cont. Dyn-B., (2014), 1667–1687.[23] L. Ju, J. Zhang, and Q. Du, Fast and accurate algorithms for simulating coarsening dynamicsof Cahn-Hilliard equations , Comp. Mater. Sci., (2015), 272–282.[24] L. Ju, J. Zhang, L. Zhu, and Q. Du,
Fast explicit integration factor methods for semilinearparabolic equations , J. Sci. Comput., (2015), 431–455.[25] R. V Kohn and X. Yan, Upper bound on the coarsening rate for an epitaxial growth model ,Commun. Pur. Appl. Math., (2003), no. 11, 1549–1564.[26] B. Li, High-order surface relaxation versus the Ehrlich-Schwoebel effect , Nonlinearity, (2006), no. 11, 2581–2603.[27] B. Li and J. Liu, Thin film epitaxy with or without slope selection , Eur. J. Appl. Math., (2003), no. 6, 713–743.[28] B. Li and J. Liu, Epitaxial growth without slope selection: energetics, coarsening, and dynamicscaling , J. Nonlinear Sci., (2004), no. 5, 429–451.[29] D. Li and Z. Qiao, On the stabilization size of semi-implicit Fourier-spectral methods for 3DCahn-Hilliard equations , Commun. Math. Sci., (2017), 1489–1506.[30] D. Li, Z. Qiao, and T. Tang, Characterizing the stabilization size for semi-implicit Fourier-spectral method to phase field equations , SIAM J. Numer. Anal., (2016), no. 3, 1653–1681.[31] W. Li, W. Chen, C. Wang, Y. Yan, and R. He, A second order energy stable linear scheme fora thin film model without slope selection , J. Sci. Comput., (2018), no. 3, 1905–1937.[32] X. Li, Z. Qiao, and H. Zhang, Convergence of a fast explicit operator splitting method forthe epitaxial growth model with slope selection , SIAM J. Numer. Anal., (2017), no. 1,265–285.[33] D. Moldovan and L. Golubovic, Interfacial coarsening dynamics in epitaxial growth with slopeselection , Phys. Rev. E, (2000), no. 6, 6190-6214.[34] L. Nirenberg, On elliptic partial differential equations , IL Principio Di Minimo E Sue Appli-cazioni Alle Equazioni Funzionali, (1959), no. 1, 1–48.[35] Z. Qiao, T. Tang, and H. Xie, Error analysis of a mixed finite element method for the molecularbeam epitaxy model , SIAM J. Numer. Anal., (2015), no. 1, 184–205.[36] Z. Qiao, C. Wang, S. M. Wise, and Z. Zhang, Error analysis of a finite difference scheme forthe epitaxial thin film model with slope selection with an improved convergence constant ,Inter. J. Numer. Anal. Mod., (2017), 283–305.ETDMs2 for no-slope-selection thin film model [37] R. L. Schwoebel, Step motion on crystal surfaces. II , J. Appl. Phys., (1969), no. 2, 614–618.[38] J. Shen, T. Tang, and L. Wang, Spectral methods: algorithms, analysis and applications , vol. 41,Springer Science & Business Media, 2011.[39] J. Shen, C. Wang, X. Wang, and S. M. Wise,
Second-order convex splitting schemes for gradientflows with Ehrlich–Schwoebel type energy: application to thin film epitaxy , SIAM J. Numer.Anal., (2012), no. 1, 105–125.[40] J. Shen, J. Xu, and J. Yang, The scalar auxiliary variable (sav) approach for gradient flows ,J. Comput. Phys., (2018), 407–416.[41] C. Wang, X. Wang, and S. M. Wise,
Unconditionally stable schemes for equations of thin filmepitaxy , Discrete Contin. Dyn. Syst., (2010), no. 1, 405–423.[42] X. Wang, L. Ju, and Q. Du, Efficient and stable exponential time differencing Runge-Kuttamethods for phase field elastic bending energy models , J. Comput. Phys., (2016), 21–38.[43] C. Xu and T. Tang,
Stability analysis of large time-stepping methods for epitaxial growthmodels , SIAM J. Numer. Anal., (2006), no. 4, 1759–1779.[44] X. Yang, J. Zhao, and Q. Wang, Numerical approximations for the molecular beam epitaxialgrowth model based on the invariant energy quadratization method , J. Comput. Phys., (2017), 104–127.[45] L. Zhu, L. Ju, and W. Zhao,
Fast high-order compact exponential time differencing Runge-Kutta methods for second-order semilinear parabolic equations , J. Sci. Comput.,67