Energy stable arbitrary order ETD-MS method for gradient flows with Lipschitz nonlinearity
aa r X i v : . [ m a t h . NA ] F e b Energy stable arbitrary order ETD-MS methodfor gradient flows with Lipschitz nonlinearity
Wenbin Chen ∗ Shufen Wang † Xiaoming Wang ‡ February 23, 2021
Abstract
We present a methodology to construct efficient high-order in timeaccurate numerical schemes for a class of gradient flows with appropriateLipschitz continuous nonlinearity. There are several ingredients to thestrategy: the exponential time differencing (ETD), the multi-step (MS)methods, the idea of stabilization, and the technique of interpolation.They are synthesized to develop a generic k th order in time efficient lin-ear numerical scheme with the help of an artificial regularization term ofthe form Aτ k ∂∂t L p ( k ) u where L is the positive definite linear part of theflow, τ is the uniform time step-size. The exponent p ( k ) is determinedexplicitly by the strength of the Lipschitz nonlinear term in relation to L together with the desired temporal order of accuracy k . To validateour theoretical analysis, the thin film epitaxial growth without slope se-lection model is examined with a fourth-order ETD-MS discretization intime and Fourier pseudo-spectral in space discretization. Our numericalresults on convergence and energy stability are in accordance with ourtheoretical results. AMS subject classifications:
Key words:
Gradient flow, epitaxial thin film growth, exponential time differ-encing, long time energy stability, arbitrary order scheme, multi-step method.
Many natural and engineering processes are gradient flows in the sense thatthe time evolution of the system is in the direction of decreasing certain energy ∗ Shanghai Key Laboratory for Contemporary Applied Mathematics, School of Mathemat-ical Sciences; Fudan University, Shanghai, China 200433 ( [email protected] ) † School of Mathematical Sciences; Fudan University, Shanghai, China 200433( ) ‡ SUSTech International Center for Mathematics and Department of Mathematics and andGuangdong Provincial Key Laboratory of Computational Science and Material Design andNational Center for Applied Mathematics Shenzhen, Southern University of Science and Tech-nology, Shenzhen, China 518055, PRC ( [email protected] ), corresponding author. Aτ ∆ N ( u n +1 − u n ) was added with A = O ( ǫ − ) to guar-antee energy stability. The purpose of this manuscript is to present a systematicapproach to construct ETD-MS based energy stable schemes of arbitrary orderin time with the help of an appropriate stabilizing term for a class of gradientflows on a Hilbert space with a positive linear part and a mild nonlinear partsatisfying Lipschitz condition in some suitable sense. More importantly, we canmake the stabilized coefficient A independent of the small parameter ǫ or time2tep-size τ (i.e., of order O (1)) by a proper choice of p ( k ) (the spatial order ofstabilized term). So far as we know, this is the first result on unconditionallyenergy stable, arbitrary-order ETD-MS based efficient algorithm for this typeof gradient flow.Let E ( u ) be a suitable energy functional on a Hilbert space H (with thedomain being a subspace of H ), and let δEδu be the variational derivative of E with respect to u in H . The associated gradient flow with mobility M ≥ ∂u∂t = − M δEδu . (1.1)Taking inner product with (1.1) and δEδu on space H we formally arrive at thefollowing energy equality: dE ( u ) dt = − M (cid:13)(cid:13)(cid:13)(cid:13) δEδu (cid:13)(cid:13)(cid:13)(cid:13) H ≤ . (1.2)This implies that the flow always evolves in the direction of decreasing energy.Notice that both the classical Allen-Cahn model and the Cahn-Hilliard modelfall into this framework with the same Ginzburg-Landau free energy but differentHilbert space H .The variational derivative of E can be split into two parts, a linear part L u ,and a nonlinear one N ( u ). In additional, large system size for models arising inmaterial sciences often corresponds the existence of a small parameter in front ofthe linear term once we non-dimensionalize the system. Formulating the systemas an abstract ODE in the Hilbert space H or an ODE system after performingappropriate spatial discretization (say Fourier pseudo-spectral in the spatiallyperiodic case) we have dudt = − ǫ L u + N ( u ) . (1.3)We will impose the following two assumptions on (1.3) :1. The operator L is non-negative.Thus we can define operators L α/ for any α ≥
0. The domain of operator L α/ is denoted by V α . For α = 1 and α = 0, it is abbreviated as D ( L ) = V and D ( L ) = H , respectively;2. The nonlinear term is Lipschitz continuous in the sense that: ∃ β > , γ > , C L > kN ( u ) − N ( v ) k V − β ≤ C L k u − v k V γ , ∀ u, v ∈ V γ , (1.4)where V − β is the dual space to V β with the duality induced by the innerproduct on H . 3he case when L is only semi-positive definite or positive definite minus a con-stant multiply of the identity operator can be treated after we shift the L sothat it is positive definite, and modify the “nonlinear” term accordingly (nowit may contain a linear part).A typical gradient flow that fits into this abstract framework is the thinfilm epitaxy growth equation without slope selection (see [21, 34, 36, 37, 40] andreferences therein): ∂u∂t = − ǫ ∆ u − ∇ · (cid:18) ∇ u |∇ u | (cid:19) (1.5)with the energy functional given by E ( u ) = Z Ω (cid:18) ǫ | ∆ u | −
12 ln (cid:0) |∇ u | (cid:1)(cid:19) d x , (1.6)and the Hilbert space H = L , the linear and nonlinear operators are L = ∆ , N ( u ) = −∇ · (cid:16) ∇ u |∇ u | (cid:17) .Most gradient flows do not fit into this framework directly. However, many ofthem enjoy invariant regions in the phase space that attract all solutions such asthe Cahn-Hilliard equation, the Allen-Cahn equation and many other reaction-diffusion equations, and the two dimensional Navier-Stokes equations [49] amongothers. For these models, we could modify the nonlinear term with the help ofan appropriate truncation to derive a “prepared” equation as in the theory ofinertial manifolds so that certain Lipschitz continuity condition is satisfied [49].Hence the framework is quite general.The main contribution of this manuscript is the construction and energy sta-bility analysis of ETD-MS based numerical schemes of arbitrary order for (1.3)under the two assumptions postulated above, with an appropriate stabilizationterm of the form Aτ k ddt L p ( k ) u ( t ) for k th order scheme where A is independentof the small parameter ǫ or the spatial/temporal discretization. The exponent p ( k ) is determined explicitly independent of ǫ or the temporal/spatial grids.These algorithms are highly efficient since only one fixed positive definite con-stant coefficient elliptic problem needs to be solved at each time step. So far aswe know, this is the first time such an efficient systematic approach is presentedalthough it was alluded to in [9].The rest of the manuscript is organized as follows. We present the numericalscheme together with the energy stability analysis in section 2. Numerical resultsare presented in section 3. Concluding remarks are offered in section 4. In this section, we propose a temporal semi-discrete scheme for (1.3), in whicha generic k th order approximation in time is constructed by applying the ex-ponential time differencing and multi-step method. The nonlinear term N ( u )is treated explicitly and the Lagrange approximation is adopted, leading to a4inear system. However, higher order multi-step treatment to nonlinear termgives rise to strong instability, thus to keep the energy decaying property, a k th order regularization term of the form of Aτ k ddt L p ( k ) u ( t ) is added with a carefulchoice of exponent p ( k ). Denote the discrete numerical solution at t = t n by u n . The differential form ofthe numerical scheme for (1.3) is to find u n +1 ( t ) : [ t n , t n +1 ] → an appropriatesubspace of H such that du n +1 ( t ) dt + ǫ L u n +1 ( t ) + Aτ k ddt L p ( k ) u n +1 ( t ) = k − X i =0 ℓ i ( t − t n ) N ( u n − i ) , (2.1)where ℓ i ( s ) are the shifted (to the negative range) Lagrange basis polynomialof degree k with the form of ℓ i ( s ) = Y ≤ m ≤ k − m = i mτ + s ( m − i ) τ = Y ≤ m ≤ k − m = i m + s/τm − i = k − X j =0 ξ i,j s j , (2.2)where { ξ i,j } k − j =0 are coefficients of the polynomial ℓ i ( s ). Obviously, ξ i,j ∼O ( τ − j ). This property will be used later in the stability analysis.Introducing an integrating factor e K t with K = ǫ (cid:16) I + Aτ k L p ( k ) (cid:17) − L andintegrating the equation (2.1) from t n to t n +1 leads to the following k th orderETD-MS scheme u n +1 = e − K τ u n + k − X i =0 Z τ e − K ( τ − s ) ℓ i ( s ) ds N ( u n − i ) . (2.3)Since the Lagrange interpolation polynomials and the linear operator L areknown a priori, the algorithm can be simplified. Let φ j := R τ e − K ( τ − s ) s j ds which can be calculated (explicitly) before hand, it can be shown via recurrenceformula φ = K − (cid:16) I − e − K τ (cid:17) ,φ j = K − (cid:0) τ j − jφ j − (cid:1) , ≤ j ≤ k − . (2.4)Thus the k th order ETD-MS scheme (2.3) can be rewritten as u n +1 = e − K τ u n + k − X i =0 k − X j =0 ξ i,j φ j N ( u n − i ) . (2.5)This is an extremely efficient algorithm, especially since we could pre-calculatethe operators involved. 5 .2 Energy stability In this subsection, we establish the energy stability for the scheme (2.1). Firstwe present an interpolation estimation that will be used later.
Lemma 2.1.
For any u ∈ V β and v ∈ V γ , p ( k ) > max { β, γ } , q ∈ (0 , , andany constants ˆ C, ˜ C to be determined for specific problem, the following inequalityholds τ (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) V β (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) V γ ≤ C (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) H + C τ qp ( k ) β (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) V p ( k ) + C (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) H + C τ − q ) p ( k ) γ (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) V p ( k ) , (2.6) where C ( ˆ C ) = 1 − β/p ( k )2 ˆ C / (1 − β/p ( k )) , C ( ˆ C ) = β p ( k ) ˆ C − p ( k ) /β ,C ( ˜ C ) = 1 − γ/p ( k )2 ˜ C / (1 − γ/p ( k )) , C ( ˜ C ) = γ p ( k ) ˜ C − p ( k ) /γ . (2.7) Proof.
Using the interpolation inequality to control k·k V β by k·k H and k·k V p ( k ) : τ (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) V β (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) V γ ≤ τ (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) − β/p ( k ) H (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) β/p ( k ) V p ( k ) · (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) − γ/p ( k ) H (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) γ/p ( k ) V p ( k ) . (2.8)Denote the term in the right hand side (RHS) of (2.8) by I , then we use theYoung’s inequality twice to show I ≤ τ q (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) − β/p ( k ) H (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) β/p ( k ) V p ( k ) + 12 τ − q (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) − γ/p ( k ) H (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) γ/p ( k ) V p ( k ) ≤ − β/p ( k )2 ˆ C (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) − β/p ( k ) H ! / (1 − β/p ( k )) + β p ( k ) τ q ˆ C (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) β/p ( k ) V p ( k ) ! p ( k ) /β + 1 − γ/p ( k )2 ˜ C (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) − γ/p ( k ) H ! / (1 − γ/p ( k )) + γ p ( k ) τ − q ˜ C (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) γ/p ( k ) V p ( k ) ! p ( k ) /γ = C (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) H + C τ qp ( k ) β (cid:13)(cid:13)(cid:13)(cid:13) dudt (cid:13)(cid:13)(cid:13)(cid:13) V p ( k ) + C (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) H + C τ − q ) p ( k ) γ (cid:13)(cid:13)(cid:13)(cid:13) dvdt (cid:13)(cid:13)(cid:13)(cid:13) V p ( k ) . (2.9)This completes the proof of Lemma 2.1.Now we come to the energy stability. For simplicity, we denote k·k L ( t i ,t j ; V α ) by k · k L ( I i,j ; V α ) hereafter. 6 emma 2.2. For system (2.1) , the following energy estimation establishes: E ( u n +1 ) − E ( u n ) + (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + Aτ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) ≤ k − X j =0 C L τ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) − j − X i = − ℓ i ( t − t n ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ) (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V γ ) (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V β ) , (2.10) where the convention ℓ − ( t − t n ) = 0 has been used.Proof. To establish the desired energy estimates, we take inner product of (2.1)with du n +1 ( t ) dt , which gives (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) H + Aτ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) V p ( k ) + ddt E ( u n +1 ( t ))= k − X i =0 ℓ i ( t − t n ) N ( u n − i ) − N ( u n +1 ( t )) , du n +1 ( t ) dt ! H . (2.11)Integrating from t n to t n +1 gives (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + Aτ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + E ( u n +1 ) − E ( u n )= Z t n +1 t n k − X i =0 ℓ i ( t − t n ) N ( u n − i ) − N ( u n +1 ( t )) , du n +1 ( t ) dt ! H dt. (2.12)Note that the sum of the Lagrange basis functions equals one, i.e., P k − i =0 ℓ i ( t − t n ) = 1, thus terms within the integral in RHS of (2.11) (denoted by NLT) can7e rewritten asNLT = k − X i =0 ℓ i ( t − t n ) (cid:0) N ( u n − i ) − N ( u n +1 ( t )) (cid:1) , du n +1 ( t ) dt ! H = k − X i =0 ℓ i ( t − t n ) (cid:0) N ( u n − i ) − N ( u n − i +1 ) + · · · + N ( u n ) − N ( u n +1 ( t )) (cid:1) , du n +1 ( t ) dt ! H = k − X i =0 ℓ i ( t − t n ) (cid:0) N ( u n ) − N ( u n +1 ( t )) (cid:1) , du n +1 ( t ) dt ! H + k − X j =1 k − X i = j ℓ i ( t − t n ) (cid:0) N ( u n − j ) − N ( u n − j +1 ) (cid:1) , du n +1 ( t ) dt H = (cid:18) N ( u n ) − N ( u n +1 ( t )) , du n +1 ( t ) dt (cid:19) H + k − X j =1 − j − X i =0 ℓ i ( t − t n ) ! (cid:0) N ( u n − j ) − N ( u n − j +1 ) (cid:1) , du n +1 ( t ) dt ! H . (2.13)Using the Cauchy-Schwartz inequality and the Lipschitz continuity (1.4) makes Z t n +1 t n (cid:18) N ( u n ) − N ( u n +1 ( t )) , du n +1 ( t ) dt (cid:19) H dt ≤ C L Z t n +1 t n (cid:13)(cid:13) u n − u n +1 ( t ) (cid:13)(cid:13) V γ (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) V β dt ≤ C L Z t n +1 t n τ (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V γ ) (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) V β dt ≤ C L τ (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V γ ) (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V β ) :=NLT , (2.14)where the last inequality follows from the H¨older inequality.Similarly, the remaining terms in RHS of (2.12) can be estimated as Z t n +1 t n − j − X i =0 ℓ i ( t − t n ) ! (cid:0) N ( u n − j ) − N ( u n − j +1 ) (cid:1) , du n +1 ( t ) dt ! H dt ≤ C L τ (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V γ ) (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V β ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) − j − X i =0 ℓ i ( t − t n ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ) :=NLT j , ≤ j ≤ k − . (2.15)This completes the proof. 8ext we give an upper bound for the L -integral in time (cid:13)(cid:13)(cid:13) − P j − i =0 ℓ i ( t − t n ) (cid:13)(cid:13)(cid:13) L ( I n,n +1 ) and further provide the energy stability for scheme (2.1). Recall that in (2.2),the Lagrange basis ℓ i ( s ) is expressed as the polynomial of s with coefficients ξ i,j . According to the properties of ξ i,j , it’s easy to see (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) − j − X i =0 ℓ i ( t − t n ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) − j − X i =0 k − X r =0 ξ i,r ( t − t n ) r (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ) = C ∗ j τ / , ≤ j ≤ k − , (2.16)where the constants C ∗ j are independent of time step-size τ or the current time t . For convenience, we follow the convention of C ∗ = 1 hereafter.We now introduce the following notation for the sake of brevity in presenta-tion. C j := k − − j X r =0 C ∗ k − − r , j = 0 , · · · , k − . (2.17)It follows that C j = C j +1 + C ∗ j , C k − = C ∗ k − . (2.18)Next, we define the following modified energy˜ E ( u n ) = E ( u n ) + C L C k − X j =1 C j (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; H ) + C L C k − X j =1 C j τ k (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V p ( k ) ) , (2.19)where both C , C depend on ˆ C , ˜ C as specified in (2.7).Thanks to (2.7), C , C can be made as small as we need so long as weset ˆ C, ˜ C small enough. Therefore, for small enough constants ˆ C , ˜ C and largeenough constant A the following inequalities hold (cid:0) − C L (cid:0) C + C C (cid:1)(cid:1) ≥ C L C C , (2.20) (cid:0) A − C L (cid:0) C + C C (cid:1)(cid:1) ≥ C L C C . (2.21)Note that C ∗ = 1, and hence C = C + 1 according to (2.18). Therefore(2.20)–(2.21) are simplified as1 ≥ C L ( C + C ) C ,A ≥ C L ( C + C ) C . Pick ˆ C and ˜ C so that C L ( C + C ) C ≤
1, i.e., (1 − β/p ( k )) ˆ C / (1 − β/p ( k )) +(1 − γ/p ( k )) ˜ C / (1 − γ/p ( k )) ≤ / ( C L C ), and then let A = C L (cid:18) β p ( k ) ˆ C − p ( k ) /β + γ p ( k ) ˜ C − p ( k ) /γ (cid:19) C ,
9e have (2.20)–(2.21).We are now ready to prove the main result of the energy stability.
Theorem 2.3.
The numerical scheme (2.1) is energy stable in the sense that ˜ E ( u n +1 ) ≤ ˜ E ( u n ) , ∀ n ≥ k, (2.22) provided that (2.20) – (2.21) are satisfied, and p ( k ) = ( β + γ ) k .Proof. By (2.16), the estimation (2.15) for NLT j , ≤ j ≤ k − j = C L C ∗ j τ (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V γ ) (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V β ) . (2.23)Applying Lemma 2.1 to (2.14) and (2.23), these nonlinear terms can be boundedfurther:NLT ≤ C L " C (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + C τ qp ( k ) β (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + C (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + C τ − q ) p ( k ) γ (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) , (2.24)NLT j ≤ C L C ∗ j " C (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + C τ qp ( k ) β (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + C (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; H ) + C τ − q ) p ( k ) γ (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V p ( k ) ) . (2.25)Choosing indexes q and p ( k ) to satisfy2 qp ( k ) β = k, − q ) p ( k ) γ = k, (2.26)and simple calculation shows q = 11 + γ/β , p ( k ) = ( β + γ ) k . (2.27)10hen estimates (2.24)–(2.25) giveNLT ≤ C L " C (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + C τ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + C (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + C τ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) , (2.28)NLT j ≤ C L C ∗ j " C (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + C τ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + C (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; H ) + C τ k (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V p ( k ) ) . (2.29)Unify the expression in (2.24)–(2.25) with the convention of C ∗ = 1 and combine(2.28)–(2.29) with (2.12), it yields (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + Aτ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + E ( u n +1 ) − E ( u n ) ≤ C L C + C k − X j =0 C ∗ j (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + C L C + C k − X j =0 C ∗ j τ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + C L C k − X j =1 C ∗ j (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; H ) + C L C τ k k − X j =1 C ∗ j (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V p ( k ) ) . (2.30)Adding C L C P k − j =1 C j +1 (cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; H ) , C L C P k − j =1 C j +1 τ k (cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V p ( k ) )
11o both sides of (2.30) and utilizing (2.18), we deduce E ( u n +1 ) + − C L C + C k − X j =0 C ∗ j (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; H ) + A − C L C + C k − X j =0 C ∗ j τ k (cid:13)(cid:13)(cid:13)(cid:13) du n +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n,n +1 ; V p ( k ) ) + C L C k − X j =1 C j +1 (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; H ) + C L C k − X j =1 C j +1 τ k (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V p ( k ) ) ≤ E ( u n ) + C L C C (cid:13)(cid:13)(cid:13)(cid:13) du n ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − ,n ; H ) + C L C C τ k (cid:13)(cid:13)(cid:13)(cid:13) du n ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − ,n ; V p ( k ) ) + C L C k − X j =2 C j (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; H ) + C L C k − X j =2 C j τ k (cid:13)(cid:13)(cid:13)(cid:13) du n − j +1 ( t ) dt (cid:13)(cid:13)(cid:13)(cid:13) L ( I n − j,n − j +1 ; V p ( k ) ) . (2.31)Then the modified energy decaying property (2.22) follows from (2.20)–(2.21)and (2.31). Remark 2.4.
The requirement postulated in (2.20) – (2.21) is sufficient but notnecessary. It would be interesting to find the optimal choice of A that guaranteesenergy stability. Our numerical results presented in the next section suggest thatthe stability is insensitive to the choice of A for a certain range well below thetheoretical requirement that we have derived here. The theoretically optimal p ( k ) comes from two consideration: (1) it should be as small as possible toavoid large artificial error; (2) it should be large enough to control the nonlinearterm when combined with the original dissipation term. The critical p ( k ) isderived under the condition that we wish to take A to be independent of ǫ .Larger p ( k ) may lead to larger error, especially for high-frequency solutions.We can use a relatively low-order regularization to reduce the artificial errors,while this treatment may give us an A that depends on ǫ or τ . In [11], a Aτ (cid:0) ∆ u n +1 − ∆ u n (cid:1) regularization term is added in their third-order ETDscheme, the exponent p (3) = 1 is smaller than our theoretical optimal value, thusthe artificial error is supposed to be smaller, while their stabilized coefficient isof order O ( ǫ − ) . This is due to the fact that the regularization combined withoriginal dissipation term is insufficient to control the explicit nonlinear termsince temporal accuracy of the artificial regularization term has to be kept, thusthe surface diffusion term is involved in the energy stability analysis which is oforder O ( ǫ ) . Remark 2.5.
Note that the regularization term Aτ k ∂∂t L p ( k ) u can be replaced bya Dupont-Douglas type Aτ k − L p ( k ) (cid:0) u n +1 − u n (cid:1) as long as the exact solution issmooth enough. The continuous form adopted here is consistent with the spirit f ET D . As we mentioned before, the introduction of high-order regularizationmay lead to large truncation error for the high-frequency solution, while theconvergence in finite time is preserved as we can see in the following numericaltests. For long-time coarsening process with random initial data, two alternativetreatment could be considered in the initial evolution until a relatively smoothsolution is obtained: 1. use some high-order methods without stabilization; 2.use variable time step-size in the simulation, which is supposed to be very smallin the beginning in order to handle the effect of high-frequency components ofsolution.
In this section, we applied the abstract framework to the NSS equation with k =4. The two dimensional domain Ω = [0 , L ] with periodic boundary conditionis considered. In this case, L = ∆ , N ( u ) = −∇ · (cid:16) ∇ u |∇ u | (cid:17) in our generalframework. Thus the abstract functional spaces are specified to H = { f : f ∈ L with zero mean } , V / = { f : f ∈ H per with zero mean } , V α = { f : f ∈ H αper with zero mean } and the Lipschitz continuity takes the following form kN ( u ) − N ( v ) k V − ≤ k u − v k V , (3.1)that is β = γ = , then equation (2.27) indicates p ( k ) = k/ du n +1 ( t ) dt + ǫ ∆ u n +1 ( t ) + Aτ ddt ∆ u n +1 ( t ) = X i =0 ℓ i ( t − t n ) N ( u n − i ) , t ∈ [ t n , t n +1 ] . (3.2)The Lagrange basis functions { ℓ i ( s ) } i =0 , ≤ s ≤ τ are ℓ ( s ) = s τ + s τ + 11 s τ + 1 , ℓ ( s ) = − s τ − s τ − sτ ,ℓ ( s ) = s τ + 2 s τ + 3 s τ , ℓ ( s ) = − s τ − s τ − s τ . The corresponding constants in (2.16) are C ∗ = p / , C ∗ = p / , C ∗ = p /
945 and C = C ∗ + C ∗ + C ∗ = (cid:16) √ √ √ (cid:17) / √ . (3.3)The constants in (2.6) are C = ˆ C / , C = ˆ C − , C = ˜ C / , C = ˜ C − .Substitute these into (2.20)–(2.21), then we can take ˆ C = ˜ C = (cid:16) C ) (cid:17) / to13atisfy 1 − C / (cid:0) C (cid:1) ≥ C / C . Thus the regularization coefficient A = ˆ C − (cid:0) C (cid:1) = ( C ) with C given in (3.3) is required. However, we will show the numerical results do notdepend much on this relatively stringent condition in the following.Note that initial steps are needed to start the high-order multi-step scheme(3.2) and here we use the ETD-RK method to compute the first three numer-ical solutions. The general expression of the fourth-order ETD-RK method isillustrated in [33] as a n = e − ǫ L τ/ u n + ( − ǫ L ) − (cid:16) e − ǫ L τ/ − I (cid:17) N ( u n , t n ) ,b n = e − ǫ L τ/ u n + ( − ǫ L ) − (cid:16) e − ǫ L τ/ − I (cid:17) N ( u n , t n + τ / ,c n = e − ǫ L τ/ a n + ( − ǫ L ) − (cid:16) e − ǫ L τ/ − I (cid:17) (2 N ( b n , t n + τ / − N ( u n , t n )) ,u n +1 = e − ǫ L τ u n + τ − ( − ǫ L ) − nh − ǫ L τ + e − ǫ L τ (cid:0) ǫ L τ + ( − ǫ L τ ) (cid:1)i N ( u n , t n )+ 2 h − ǫ L τ − e − ǫ L τ (2 + ǫ L τ ) i ( N ( a n , t n + τ /
2) + N ( b n , t n + τ / h − ǫ L τ − ( − ǫ L τ ) + e − ǫ L τ (4 + ǫ L τ ) i N ( c n , t n + τ ) o . (3.4)To solve system (3.2), the spatial discretization is performed by the Fourierpseudo-spectral method with a resolution N = 128 and this can be efficientlyimplemented via the fast Fourier transform. The convergence and energy sta-bility tests are provided. In this subsection, the fourth order temporal convergence of scheme (3.2) isverified. The parameters are chosen as L = 2 π , ǫ = 0 .
01 and terminal time T = 1. To test the convergence, an artificial forcing term g is added in the righthand side of (3.2) to make the exact solution u ( t ) = e − t cos(2 x ) cos(2 y ): g = ( − ε ) u − u e − t [1 − cos(4 x ) cos(4 y )]+ 16 e − t u [1 + 2 e − t (1 − cos(4 x ) cos(4 y ))] [cos(4 x ) + cos(4 y ) − x ) cos(4 y )] , and the discrete L error is calculated. Firstly, the effect of regularized coef-ficient A on the convergence is examined in Table 1. From which we can see,the errors grow linearly with the increasing of A and the convergence order ispreserved for all of the choice of A . Then let A be fixed, results for p ( k ) = 1 . p ( k ) = 1 . p ( k ) = 2 . p ( k ) = 2 . p ( k ).Table 1: Temporal convergence of (3.2) with p ( k ) = 2. τ A = 1 A = 5 A = 10 A = ( C ) error order error order error order error order2.50E-03 6.53e-07 3.26e-06 6.52e-06 1.12e-041.25E-03 4.10e-08 3.993 2.05e-07 3.992 4.10e-07 3.991 7.18e-06 3.9596.25E-04 2.57e-09 3.997 1.28e-08 3.996 2.57e-08 3.996 4.50e-07 3.9943.13E-04 1.59e-10 4.009 8.04e-10 3.999 1.61e-09 4.000 2.82e-08 3.9981.56E-04 8.94e-12 4.157 4.96e-11 4.017 9.99e-11 4.007 1.76e-09 4.000Table 2: Temporal convergence of (3.2) with A = ( C ) . τ p ( k ) = 1 . p ( k ) = 1 . p ( k ) = 2 . p ( k ) = 2 . In this subsection, the physically interesting coarsening process is simulated.The parameters are now set as L = 12 . ǫ = 0 . T = 30000 and τ = 10 − .Two choices for the regularization coefficient A = ( C ) and A = 10 aretested. The scaling laws for the energy E , average surface roughness h and theaverage slope m will be shown, that is E ∼ O ( − ln( t )), h ∼ O ( t ) and m ∼ O ( t ) as t → ∞ . (See [21, 36, 37] and references therein). The correspondingdefinitions for these physical quantities are E ( u ) = (cid:18) −
12 ln(1 + |∇ u | ) , (cid:19) + ε k ∆ u k ,h ( u, t ) = s h | Ω | X x | u ( x, t ) − ¯ u ( t ) | , with ¯ u ( t ) := h | Ω | X x u ( x, t ) ,m ( u, t ) = s h | Ω | X x |∇ u ( x i,j , t ) | . The snapshots of the numerical solution (3.2) at time t = 1, 5000, 10000,15000, 20000, 30000 with A = ( C ) and A = 10 are displayed in Figures 115nd 2 respectively. The evolution of E , h and m with two choices of A is estab-lished in Figures 3–5, and the linear fitting results for the solution of (3.2) intime interval [1 , A with a uniform step-size τ = 10 − , in Figures 1–5. Since random initial datais used, the artificial error arising from high-order regularization term may belarge in the initial stage, resulting in the sensitivity to stabilized coefficient A .Therefore, an additional coarsening simulation is performed in the followingwith a variable time step-size, which is set as: τ = 10 − for t ≤ τ = 10 − for 1 < t ≤ τ = 10 − for 10 < t ≤
100 and τ = 10 − for t >
100 andother parameters are kept. The related results are displayed in Figures 6–7and Figures 9–11. From which we can see, the sensitivity to A is significantlyreduced. Also, minor differences of solutions between two choices of A areobserved in Figure 8. Remark 3.1.
Thanks to the anonymous reviewer for pointing out the differ-ences between snapshots of numerical solution for two choices of A at the sametime level, in Figures 1–2. In fact, the introduction of high-order regularizationterm may result in large truncation error in the beginning since random ini-tial data is chosen, and further lead to different steady phase states. The samesimulation of coarsening process is performed with variable time step-sizes andother parameters kept. The time step-size is taken to be small to handle theeffect of high-order regularization during the initial evolution of solution anda relatively large step-size is adopted after a smooth solution obtained. Withthis treatment being taken, the insensitivity to the stabilized coefficient A of ournumerical scheme can be observed in Figures 6–11. This exactly shows the neces-sity of variable step-size in the long-time simulation of coarsening process. Thechoice of variable step-size is empirical here, and similar idea has been appliedin [38, 51]. There are also various adaptive methods regarding to the choice oftime step-size, see [10, 20, 42, 54]. In fact, one of our authors, Xiaoming Wang,has already pointed out the feasibility and necessity of using adaptive strategiesor hybrid approach (utilize some alternative high-order methods without regular-ization for the initial stage) for long-time computation in our previous paper [9],and it would be our future work to conduct the adaptive time-stepping strategywith a posterior estimate. A highly efficient k th order accurate and linear numerical scheme is proposed forgradient flows with mild nonlinearity by utilizing ETD and MS methods togetherwith Lagrange interpolation and stabilization. The nonlinear term is treatedexplicitly and a k th order Dupont-Douglas type regularization Aτ k ∂∂t L p ( k ) u isadded for stability. The constant A is independent of the small parameter ǫ or the discretization in space or in time. Lipschitz continuity for the nonlinear16 = 5000 t = 10000 t = 15000 t = 20000 t = 30000 Figure 1: Snapshots of the numerical solutions of scheme (3.2) with A = ( C ) , τ = 10 − .term is assumed. The instability caused by the explicit treatment can be to-tally overcome with the aid of regularization. The exponent p ( k ) = ( β + γ ) k isdetermined explicitly by Lipschitz condition and the order of the scheme. Asan example, a fourth order ETD-MS method is applied to a thin film epitaxymodel without slope selection, some numerical experiments have been presentedto validate the fourth order convergence in time and the unconditional long-timeenergy stability. The method can be applied to a wide range of gradient flowsafter proper “preparation” of the original equation.It is easy to see that our method can be generalized to the case of variable-step without any difficulty. All we need to do is to use a more general Lagrangeinterpolation polynomial. This opens a door to the development of higher or-der in time efficient time-adaptive strategies [39]. The Lagrange interpolationpolynomials could be replaced by other appropriate interpolations so long asthe near decomposition of the identity and the bounded properties are satis-fied. The details together with the convergence analysis will be reported in asubsequent work. Acknowledgements
This work is supported in part by the following grants: NSFC12071090, Shang-hai Science and Technology Research Program 19JC1420101 and a 111 projectB08018 (W. Chen), NSFC11871159, Guangdong Provincial Key Laboratory forComputational Science and Material Design 2019B030301001 (X. Wang). Theauthors thank the anonymous reviewers for their comments and suggestions.17 = 5000 t = 10000 t = 15000 t = 20000 t = 30000
Figure 2: Snapshots of the numerical solutions of scheme (3.2) with A = 10, τ = 10 − . References [1] E. O. Asante-Asamani, A. Kleefeld, and B. A. Wade. A second-order expo-nential time differencing scheme for non-linear reaction-diffusion systemswith dimensional splitting.
Journal of Computational Physics , 109490,2020.[2] S. M. Allen and J. W. Cahn. A microscopic theory for antiphase boundarymotion and its application to antiphase domain coarsening.
Acta Metallur-gica , 27(6):1085–1095, 1979.[3] L. Ambrosio and N. Gigli and G. Svare. Gradient Flows in metric spacesand in the space of probability measures, 2nd ed. Birkhauser, 2008.[4] D. M. Anderson, G. B. McFadden, and A. A. Wheeler. Diffuse-interfacemethods in fluid mechanics.
Annual Review of Fluid Mechanics , 30(1):139–165, 1998.[5] G. Beylkin, J. M. Keiser, and L. Vozovoi. A new class of time discretizationschemes for the solution of nonlinear PDEs.
Journal of ComputationalPhysics , 147(2):362–387, 1998.[6] L. A. Caffarelli and N. E. Muler. An L ∞ bound for solutions of the Cahn-Hilliard equation. ArRMA , 133(2):129–144, 1995.[7] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. i.interfacial free energy.
The Journal of Chemical Physics , 28(2):258–267,1958.[8] W. Chen, W. Li, Z. Luo, C. Wang, and X.Wang. A stabilized second orderexponential time differencing multistep method for thin film growth model18 -1 Time -500-450-400-350-300-250-200-150-100-500 E ne r g y Evolution of energy
Figure 3: Semi-log plot of the energy E of scheme scheme (3.2) with τ = 10 − .Fitted line has the form a ln( t ) + b , with coefficients a = − . b = − . A = 10 and a = − . b = − .
35 for A = ( C ) .without slope selection. ESAIM: Mathematical Modelling and NumericalAnalysis , 54(3):727–750, 2020.[9] W. Chen, W. Li, C. Wang, S. Wang, and X. Wang. Energy stable higher-order linear ETD multi-step methods for gradient flows: application to thinfilm epitaxy. Special issue dedicated to Professor Andrew Majda on theoccasion of his seventieth birthday.
Research in the Mathematical Sciences ,7(3):1–27, 2020.[10] W. Chen, X. Wang, Y. Yan, and Z. Zhang. A Second Order BDF NumericalScheme with Variable Steps for the Cahn-Hilliard Equation.
SIAM Journalon Numerical Analysis , 57(1):495–525, 2019.[11] K. Cheng, Z. Qiao, and C. Wang. A third order exponential time differenc-ing numerical scheme for no-slope-selection epitaxial thin film model withenergy stability.
Journal of Scientific Computing , 81(1):154–185, 2019.[12] S. M. Cox and P. C. Matthews. Exponential time differencing for stiffsystems.
Journal of Computational Physics , 176(2):430–455, 2002.[13] M. Doi and S. F. Edwards.
The theory of polymer dynamics , volume 73.Oxford university press, 1988. 19 -1 Time -1 A v e r age r oughne ss Evolution of average roughness
Figure 4: The log-log plot of the average surface roughness h of (3.2) with τ = 10 − . Fitted lines have the form at b , with coefficients a = 0 . b = 0 . A = 10 and a = 0 . b = 0 . A = ( C ) .[14] Q. Du, L. Ju, X. Li, and Z. Qiao. Maximum principle preserving exponentialtime differencing schemes for the nonlocal Allen–Cahn equation. SIAMJournal on Numerical Analysis , 57(2):875–898, 2019.[15] Q. Du and W. Zhu. Stability analysis and application of the exponentialtime differencing schemes.
Journal of Computational Mathematics , 200–209, 2004.[16] Q. Du and W. Zhu. Analysis and applications of the exponential time differ-encing schemes and their contour integration modifications.
BIT NumericalMathematics , 45(2):307–328, 2005.[17] K. Elder, M. Katakowski, M. Haataja, and M. Grant. Modeling elasticityin crystal growth.
Physical Review Letters , 88(24):245701, 2002.[18] C. M. Elliott and A. Stuart. The global dynamics of discrete semilinearparabolic equations.
SIAM Journal on Numerical Analysis , 30(6):1622–1663, 1993.[19] D. J. Eyre. Unconditionally gradient stable time marching the Cahn-Hilliard equation. In
Materials Research Society Symposium Proceedings ,529:39–46, 1998. 20 -1 Time M ound s l ope Evolution of average mound slope
Figure 5: The log-log plot of the average slope m of (3.2) with τ = 10 − . Fittedlines have the form at b , with coefficients a = 2, b = 0 . A = 10 and a = 2 . b = 0 . A = ( C ) .[20] X. Feng, T. Tang, and J. Yang. Long time numerical simulations for phase-field problems using p-adaptive spectral deferred correction methods. SIAMJournal on Scientific Computing , 37(1):A271–A294, 2015.[21] L. Golubovic. Interfacial coarsening in epitaxial growth models withoutslope selection.
Physical Review Letters , 78(1):90–93, 1997.[22] Y. Gong, J. Zhao, and Q. Wang. Arbitrarily high-order unconditionallyenergy stable sav schemes for gradient flow models.
Computer PhysicsCommunications , 249:107033, 2020.[23] M. E. Gurtin, D. Polignone, and J. Vinals. Two-phase binary fluids andimmiscible fluids described by an order parameter.
Mathematical Modelsand Methods in Applied Sciences , 6(06):815–831, 1996.[24] E. Hairer, S. P. Noersett, and G. Wanner. Solving ordinary differentialequations i. nonstiff problems.
Springer Ser. Comput. Math , 8, 1993.[25] M. Hochbruck and A. Ostermann. Explicit exponential Runge-Kutta meth-ods for semilinear parabolic problems.
SIAM Journal on Numerical Anal-ysis , 43(3):1069–1090, 2005.[26] M. Hochbruck and A. Ostermann. Exponential integrators.
Acta Numer. ,19(May):209–286, 2010. 21 = 5000 t = 10000 t = 15000 t = 20000 t = 30000
Figure 6: Snapshots of the numerical solutions of scheme (3.2) with A = ( C ) and variable time step-size: τ = 10 − for t ≤ τ = 10 − for1 < t ≤ τ = 10 − for 10 < t ≤
100 and τ = 10 − for t > BIT Numerical Mathematics , 51(4):889–908, 2011.[28] J. Huang, L. Ju, and B. Wu. A fast compact exponential time differenc-ing method for semilinear parabolic equations with Neumann boundaryconditions.
Applied Mathematics Letters , 94, 257–265, 2019.[29] L. Ju, X. Li, Z. Qiao, and H. Zhang. Energy stability and error estimatesof exponential time differencing schemes for the epitaxial growth modelwithout slope selection.
Mathematics of Computation , 87(312):1859–1885,2018.[30] L. Ju, X. Liu, and W. Leng. Compact implicit integration factor methodsfor a family of semilinear fourth-order parabolic equations.
Discrete &Continuous Dynamical Systems-B , 19(6):1667–1687, 2014.[31] L. Ju, J. Zhang, and Q. Du. Fast and accurate algorithms for simulatingcoarsening dynamics of Cahn-Hilliard equations.
Computational MaterialsScience , 108:272–282, 2015.[32] L. Ju, J. Zhang, L. Zhu, and Q. Du. Fast explicit integration factor methodsfor semilinear parabolic equations.
Journal of Scientific Computing , 62(2):431–455, 2015.[33] A. Kassam and L. N. Trefethen. Fourth-order time-stepping for stiff PDEs.
SIAM Journal on Scientific Computing , 26(4):1214–1233, 2005.[34] R. V. Kohn and X. Yan. Upper bound on the coarsening rate for anepitaxial growth model.
Communications on Pure and Applied Mathemat- = 5000 t = 10000 t = 15000 t = 20000 t = 30000 Figure 7: Snapshots of the numerical solutions of scheme (3.2) with A = 10 andvariable time step-size: τ = 10 − for t ≤ τ = 10 − for 1 < t ≤ τ = 10 − for 10 < t ≤
100 and τ = 10 − for t > ics: A Journal Issued by the Courant Institute of Mathematical Sciences ,56(11):1549–1564, 2003.[35] F. M. Leslie. Theory of flow phenomena in liquid crystals. Advances inliquid crystals , 4:1–81, 1979.[36] B. Li and J. Liu. Thin film epitaxy with or without slope selection.
Euro-pean Journal of Applied Mathematics , 14(6):713–743, 2003.[37] B. Li and J. Liu. Epitaxial growth without slope selection: Energetics,coarsening, and dynamic scaling.
Journal of Nonlinear Science , 14(5):429–451, 2004.[38] W. Li, W. Chen, C. Wang, Y. Yan, and R. He. A second order energystable linear scheme for a thin film model without slope selection.
Journalof Scientific Computing , 76(3): 1905–1937, 2018.[39] F. Lou, T. Tang and H. Xie. Parameter-free time adaptivity based onenergy evolution for the Cahn-Hilliard equation.
Communications in Com-putational Physics , 19(5):1542–1563, 2016.[40] D. Moldovan and L. Golubovic. Interfacial coarsening dynamics in epitaxialgrowth with slope selection.
Physical Review E , 61(6):6190, 2000.[41] Y. Morita and K. Tachibana. An entire solution to the Lotka-Volterracompetition-diffusion equations.
SIAM Journal on Mathematical Analysis ,40(6):2217–2240, 2009.[42] Z. Qiao, Z. Zhang, and T. Tang. An adaptive time-stepping strategy for themolecular beam epitaxy models.
SIAM Journal on Scientific Computing ,33(3):1395–1414, 2011. 23 = 1 t = 5000 t = 10000 t = 15000 t = 20000 t = 30000
Figure 8: Differences of the numerical solutions of scheme (3.2) between A = 10and A = ( C ) at the same time level for variable time step-size: τ = 10 − for t ≤ τ = 10 − for 1 < t ≤ τ = 10 − for 10 < t ≤
100 and τ = 10 − for t > L -normin these six snapshots are 10 − , − , − , − , − , − respectively.The magnitude of numerical solution at time t = 1 is 1 and 100 for the rest.[43] J. Shen, C. Wang, X. Wang, and S. M. Wise. Second-order convex splittingschemes for gradient flows with Ehrlich-Schwoebel type energy: applicationto thin film epitaxy. SIAM Journal on Numerical Analysis , 50(1):105-125,2012.[44] J. Shen, J. Xu, and J. Yang. The scalar auxiliary variable (SAV) approachfor gradient flows.
Journal of Computational Physics , 353:407–416, 2018.[45] J. Shen, J. Xu, and J. Yang. A new class of efficient and robust energystable schemes for gradient flows.
SIAM Review , 61(3):474–506, 2019.[46] J. Shen and X. Yang. Numerical approximations of Allen-Cahn andCahn-Hilliard equations.
Discrete & Continuous Dynamical Systems-A ,28(4):1669, 2010.[47] J. Shin, H. G. Lee, and J. Y. Lee. Unconditionally stable methods for gra-dient flow using convex splitting Runge-Kutta scheme.
Journal of Compu-tational Physics , 347:367–381, 2017.[48] Y. Takeuchi.
Global dynamical properties of Lotka-Volterra systems . WorldScientific, 1996.[49] R. Temam.
Infinite-dimensional dynamical systems in mechanics andphysics , volume 68. Springer Science & Business Media, 1997.24 -1 Time -500-450-400-350-300-250-200-150-100-500 E ne r g y Evolution of energy
Figure 9: Semi-log plot of the energy E of scheme scheme (3.2) with τ = 10 − for t ≤ τ = 10 − for 1 < t ≤ τ = 10 − for 10 < t ≤
100 and τ = 10 − for t > a ln( t ) + b , with coefficients a = − . b = − .
48 for both A = 10 and A = ( C ) .[50] C. Wang, X. Wang, and S. M. Wise. Unconditionally stable schemes forequations of thin film epitaxy. Discrete & Continuous Dynamical Systems-A , 28(1):405, 2010.[51] S. Wang, W. Chen, H. Pan, and C. Wang. Optimal rate convergenceanalysis of a second order scheme for a thin film model with slope selection.
Journal of Computational and Applied Mathematics , 377(15):112855, 2020.[52] X. Wang, L. Ju, and Q. Du. Efficient and stable exponential time differ-encing Runge-Kutta methods for phase field elastic bending energy models.
Journal of Computational Physics , 316:21–38, 2016.[53] P. Yue, J. J. Feng, C. Liu, and J. Shen. A diffuse-interface method forsimulating two-phase flows of complex fluids.
Journal of Fluid Mechanics ,515:293, 2004.[54] Z. Zhang and Z. Qiao. An adaptive time-stepping strategy for the Cahn-Hilliard equation.
Communications in Computational Physics , 11(4):1261–1278, 2012. 25 -1 Time -1 A v e r age r oughne ss Evolution of average roughness
Figure 10: The log-log plot of the average surface roughness h of (3.2) with τ = 10 − for t ≤ τ = 10 − for 1 < t ≤ τ = 10 − for 10 < t ≤ τ = 10 − for t > at b , with coefficients a = 0 . b = 0 . A = 10 and A = ( C ) .26 -1 Time M ound s l ope Evolution of average mound slope
Figure 11: The log-log plot of the average slope m of (3.2) with τ = 10 − for t ≤ τ = 10 − for 1 < t ≤ τ = 10 − for 10 < t ≤
100 and τ = 10 − for t > at b , with coefficients a = 2 . b = 0 . A = 10 and A = ( C )512