Arbitrarily High-order Linear Schemes for Gradient Flow Models
AArbitrarily High-order Linear Schemes for Gradient Flow Models
Yuezheng Gong ∗ , Jia Zhao † and Qi Wang ‡ Abstract
We present a paradigm for developing arbitrarily high order, linear, unconditionally energy stablenumerical algorithms for gradient flow models. We apply the energy quadratization (EQ) technique toreformulate the general gradient flow model into an equivalent gradient flow model with a quadraticfree energy and a modified mobility. Given solutions up to t n = n ∆ t with ∆ t the time step size, welinearize the EQ-reformulated gradient flow model in ( t n , t n +1 ] by extrapolation. Then we employ analgebraically stable Runge-Kutta method to discretize the linearized model in ( t n , t n +1 ]. Then we usethe Fourier pseudo-spectral method for the spatial discretization to match the order of accuracy in time.The resulting fully discrete scheme is linear, unconditionally energy stable, uniquely solvable, and mayreach arbitrarily high order. Furthermore, we present a family of linear schemes based on prediction-correction methods to complement the new linear schemes. Some benchmark numerical examples aregiven to demonstrate the accuracy and efficiency of the schemes. Keywords:
Energy stable schemes, gradient flow models, Runge-Kutta methods, linear high-orderschemes, pseudo-spectral methods.
For many import phenomena in physics, life science, and engineering, the processes are driven byminimizing the free energy or maximizing entropy, i.e., dissipative dynamics. Gradient flow models areusually used to model these phenomena. The generic form of a gradient flow model is given by ∂∂t
Φ = G δFδ Φ , (1.1)with proper boundary conditions. Here Φ is the thermodynamical variable, F is the free energy for theisothermal system (or entropy for the nonisothermal system), and G is the mobility operator/matrix. Thegradient flow model is thermodynamically consistent if it yields a positive entropy production or negativeenergy dissipation rate. The classical Allen-Cahn equation [3] and Cahn-Hilliard equation [4] are twoexamples of gradient flow models (1.1). Other gradient flow models include the molecular beam epitaxymodel [8], the phase-field crystal model [10], the thermodynamically consistent dendritic growth model [35],the surfactant model [33], the diblock copolymer model [7] etc.Most gradient flow models are nonlinear, so that their analytical solutions are intractable. Hence,designing accurate, efficient, and stable algorithms to solve them becomes essential [11,13,32,34,37,43,44].A numerical scheme that preserves the energy dissipation property is known as an energy stable scheme[11]. It has been shown that schemes that are not energy stable could lead to instability or oscillatory ∗ College of Science, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China; Email:[email protected] † Department of Mathematics & Statistics, Utah State University, Logan, UT 84322, USA; Email: [email protected] ‡ Corresponding author. Beijing Computational Science Research Center, Beijing 100193, China; Department of Mathe-matics, University of South Carolina, Columbia, SC 29208, USA Email: [email protected] a r X i v : . [ m a t h . NA ] O c t olutions [13]. This is because non-energy-stable schemes may introduce truncation errors that destroy thephysical law numerically. Thus, developing energy stable numerical algorithms is necessary for accuratelyresolving the dynamics of gradient flow models [20, 21, 23, 30, 44].Over the years, the development of numerical algorithms has been done primarily on a specific gradientflow model, exploiting its specific mathematical properties and structures. The noticeable ones include theAllen-Cahn and Cahn-Hilliard equation [1, 9, 11, 20, 24, 32, 36, 42] as well as the molecular beam epitaxymodel [6,12,26–28,41]. As a result, most numerical algorithms developed for a specific gradient flow modelcan hardly be applied to another gradient flow model with a different free energy and mobility. The status-quo did not change until the energy quaratization (EQ) approach was introduced a few years ago [39, 44],which turns out to be rather general so that it can be readily applied to general gradient flow modelswith no restrictions on the specific form of the mobility and free energy [16, 40] (so long as it is boundedbelow, which is usually the case). Based on the idea of EQ, the scalar auxiliary variable (SAV) approachwas introduced later [31], where each time step only linear systems with constant coefficients need to besolved. Several other extensions of EQ and SAV approaches have been explored further. For instance,regularization terms [5] and stabilization terms [38] are added, and the modified energy quadratizationtechnique [22, 25] are introduced to improve the EQ methodology. We notice that most existing schemesrelated to the EQ methodology are up to 2nd order accurate in time. Some higher-order ETD schemeshave been introduced to solve the MBE model and Cahn-Hilliard models recently, but their theoreticalproofs of energy stability are still missing. Shen et al. [30] remarked on higher-order BDF schemes usingthe SAV approach, but no rigorous theoretical proofs are available.Recently, Gong et al. [15, 18, 19] introduced the arbitrarily high-order schemes for solving gradient flowmodels by combining the methodology of energy quadratization with quadratic invariant Runge-Kutta(RK) methods. This seminal idea sheds light on solving gradient flow models with an arbitrarily highorder of accuracy. We note that the proposed high order schemes are fully nonlinear so that the solutionexistence and uniqueness are not guaranteed when the time step is large. Moreover, the implementationcan be complicated compared to linear schemes. Often, it requires iterations at each time step, which addsup to the computational cost.In this paper, we will address these issues by proposing a new paradigm for developing arbitrarily highorder schemes that are unconditionally energy stable and linear. As a significant advance over our previouswork, the newly proposed paradigm would result in linear schemes, while preserving unconditional energystability. These newly proposed schemes will bring in significant improvement in numerical implementationand reduction in computational cost. First of all, as the schemes are all linear, only linear systems needto be solved at each time step. Therefore, they are easy to implement and computationally efficient. Inaddition, the existence and uniqueness of the numerical solutions can be guaranteed, which is actuallyindependent of the time steps. In other words, larger time steps can be readily applied so long as theaccuracy requirement is met. This is warranted by the existence of a unique solution and unconditionallyenergy stable. Equipped with the benefits of EQ method and RK method, the newly proposed schemesapply to general gradient flow models. For illustration purposes, we solve the Cahn-Hilliard model and themolecular beam epitaxy model to demonstrate the effectiveness of some selected new schemes. To showtheir accuracy and efficiency, we also compare two proposed schemes with the 2nd order convex splittingschemes.The rest of this paper is organized as follows. In Section 2, we briefly introduce the general gradientflow model and use the EQ method to reformulate it into an equivalent form. In Section 3, the arbitrarilyhigh-order linear energy stable schemes are introduced, and their energy stability and uniquely solvabilityare discussed. In Section 4, several numerical examples are shown to illustrate the power of our proposedarbitrarily high order schemes. In the end, some concluding remarks are given.2 Gradient Flow Models and Their EQ Reformulation
In this section, we present the general gradient flow model firstly and then apply the energy quadrati-zation technique to reformulate the model into an equivalent gradient flow form with a quadratic energyfunctional, a modified mobility matrix and the corresponding energy dissipation law, which is called theEQ reformulated model. The EQ reformulation for this class of gradient flow models provides an elegantplatform for developing arbitrarily high-order unconditionally energy stable schemes [15, 19].
Mathematically, the form of a general gradient flow model is given by [31, 44] ∂∂t
Φ = G δFδ Φ , (2.1)where Φ = ( φ , · · · , φ d ) T is the state variable vector, G is the d × d mobility matrix operator which candepend on Φ, F is the free energy, and δFδ Φ is the variational derivative of the free energy functional withrespect to the state variable, known as the chemical potential. The triple (Φ , G , F ) uniquely defines thegradient flow model. For (2.1) to be thermodynamically consistent, the time rate of change of the freeenergy must be non-increasing: dFdt = (cid:18) δFδ Φ , ∂ Φ ∂t (cid:19) = (cid:18) δFδ Φ , G δFδ Φ (cid:19) ≤ , (2.2)where the inner product is defined by ( f , g ) = d (cid:80) i =1 (cid:82) Ω f i g i d x , ∀ f , g ∈ (cid:0) L (Ω) (cid:1) d , which requires G to benegative semi-definite. The L norm is defined as (cid:107) f (cid:107) = (cid:112) ( f , g ). Note that the energy dissipation law(2.2) holds only for suitable boundary conditions. Such boundary conditions include periodic boundaryconditions and the other boundary conditions that make the boundary integrals resulted from the integra-tion by parts vanish in the calculation of variational derivatives. In this paper, we limit our study to theseboundary conditions. We reformulate the gradient flow model (2.1) by transforming the free energy into a quadratic formusing nonlinear transformations. For the purpose of illustration, we assume the free energy is given by thefollowing F = 12 (Φ , L Φ) + (cid:0) f (Φ , ∇ Φ) , (cid:1) , (2.3)where L is a linear, self-adjoint, positive semi-definite operator (independent of Φ), and f is the bulk partof the free energy density, which has a lower bound. Then the free energy F can be rewritten into aquadratic form F = 12 (Φ , L Φ) + (cid:107) q (cid:107) − C | Ω | , (2.4)by introducing an auxiliary variable q = (cid:112) f (Φ , ∇ Φ) + C , where C is a positive constant large enough tomake q real-valued for all Φ.We denote g [Φ] = (cid:112) f (Φ , ∇ Φ) + C . Then model (2.1) can be reformulated into the following equivalent3ystem Φ t = G (cid:18) L Φ + 2 q ∂g∂ Φ − ∇ · (cid:16) q ∂g∂ ∇ Φ (cid:17)(cid:19) ,q t = ∂g∂ Φ · Φ t + ∂g∂ ∇ Φ · ∇ Φ t , (2.5)with initial conditions Φ( x ,
0) = Φ ( x ) , q ( x ,
0) = (cid:113) f (cid:0) Φ ( x ) , ∇ Φ ( x ) (cid:1) + C. (2.6)It is readily to prove that the reformulated system (2.5) preserves the following energy dissipation law d F dt = (cid:32) L Φ + 2 q ∂g∂ Φ − ∇ · (cid:16) q ∂g∂ ∇ Φ (cid:17) , G (cid:18) L Φ + 2 q ∂g∂ Φ − ∇ · (cid:16) q ∂g∂ ∇ Φ (cid:17)(cid:19)(cid:33) ≤ . (2.7)We introduce u = (Φ , q ) T (2.8)and recast system (2.5) into a compact gradient flow form u t = M δ F δ u , (2.9)with a modified mobility operator M = ∂g∂ Φ + ∂g∂ ∇ Φ · ∇ G (cid:18) , ∂g∂ Φ − ∇ · ∂g∂ ∇ Φ (cid:19) . (2.10)The energy dissipation law given in (2.7) is recast into d F dt = (cid:16) δ F δ u , M δ F δ u (cid:17) ≤ . (2.11)Since the EQ-reformulated form in (2.5) has a quadratic free energy, we next discuss how to devise linearhigh-order energy stable schemes for it. In this section, we first derive a high-precision linear gradient-flow system to approximate EQ-reformulatedmodel (2.5) up to t n = n ∆ t , where ∆ t is the time step. In particular, the corresponding energy dissipationlaw is inherited. Then the algebraically stable RK method [2] is applied to the resulting linear gradient-flow system to develop a class of linear semi-discrete schemes in time. We name the schemes linear energyquadratizatized Runge-Kutta (LEQRK) methods. In order to improve accuracy and stability, a prediction-correction technique is proposed for the LEQRK schemes, leading to the LEQRK-PC methods. These newalgorithms are linear, unconditionally energy stable, and can be devised at any desired order in time. Assuming numerical solutions of Φ up to t ≤ t n have been obtained, we then solve system (2.5) in t ∈ ( t n , t n +1 ] approximately. We utilize the numerical solutions of Φ at t ≤ t n to obtain its interpolatingpolynomial approximation denoted by Φ N ( t ). Then we approximate model (2.5) in ( t n , t n +1 ] using the4ollowing linear, variable coefficient gradient flow system Φ t = G (cid:18) L Φ + 2 q ∂g N ∂ Φ − ∇ · (cid:16) q ∂g N ∂ ∇ Φ (cid:17)(cid:19) ,q t = ∂g N ∂ Φ · Φ t + ∂g N ∂ ∇ Φ · ∇ Φ t , (3.1)where ∂g N ∂ Φ = ∂g∂ Φ [Φ N ( t )] and ∂g N ∂ ∇ Φ = ∂g∂ ∇ Φ [Φ N ( t )] are independent of Φ. The linear gradient flow system(3.1) satisfies the following energy dissipation law d F dt = (cid:32) L Φ + 2 q ∂g N ∂ Φ − ∇ · (cid:16) q ∂g N ∂ ∇ Φ (cid:17) , G (cid:18) L Φ + 2 q ∂g N ∂ Φ − ∇ · (cid:16) q ∂g N ∂ ∇ Φ (cid:17)(cid:19)(cid:33) ≤ . (3.2)Applying a s -stage RK method for the linear system (3.1), we obtain the following LEQRK scheme. Scheme 3.1 ( s -stage LEQRK Scheme) . Let b i , a ij ( i, j = 1 , · · · , s ) be real numbers and let c i = s (cid:80) j =1 a ij .For given (Φ n , q n ) and Φ N ( t n + c i ∆ t ) , ∀ i , the following intermediate values are first calculated by Φ ni = Φ n + ∆ t s (cid:80) j =1 a ij k nj ,Q ni = q n + ∆ t s (cid:80) j =1 a ij l nj ,k ni = G (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19) ,l ni = (cid:16) ∂g∂ Φ (cid:17) n, ∗ i · k ni + (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i · ∇ k ni , i = 1 , · · · , s, (3.3) where (cid:16) ∂g∂ Φ (cid:17) n, ∗ i = ∂g∂ Φ [Φ N ( t n + c i ∆ t )] and (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i = ∂g∂ ∇ Φ [Φ N ( t n + c i ∆ t )] . Then (Φ n +1 , q n +1 ) is updatedvia Φ n +1 = Φ n + ∆ t s (cid:88) i =1 b i k ni , (3.4) q n +1 = q n + ∆ t s (cid:88) i =1 b i l ni . (3.5) Definition . Denote a symmetric matrix M with the elements M ij = b i a ij + b j a ji − b i b j . A RK method is said to be algebraically stable if its RK coefficients satisfystability condition b i ≥ , ∀ i = 1 , , · · · , s, and M is positive semi-definite . (3.6)Next, we show that the algebraically stable LEQRK scheme is unconditionally energy stable. Theorem 3.1.
The LEQRK scheme with their RK coefficients satisfying stability condition (3.6) is un-conditionally energy stable, i.e., it satisfies F n +1 ≤ F n , (3.7) where F n = (Φ n , L Φ n ) + (cid:107) q n (cid:107) − C | Ω | . roof. Denoting Φ n +1 = Φ n + ∆ t s (cid:80) i =1 b i k ni and noticing that operator L is linear and self-adjoint, we have12 (Φ n +1 , L Φ n +1 ) −
12 (Φ n , L Φ n ) = ∆ t s (cid:88) i =1 b i ( k ni , L Φ n ) + ∆ t s (cid:88) i,j =1 b i b j ( k ni , L k nj ) . (3.8)Applying Φ n = Φ ni − ∆ t s (cid:80) j =1 a ij k nj to the right of (3.8), we deduce12 (Φ n +1 , L Φ n +1 ) −
12 (Φ n , L Φ n ) = ∆ t s (cid:88) i =1 b i ( k ni , L Φ ni ) − ∆ t s (cid:88) i,j =1 M ij ( k ni , L k nj ) , (3.9)where s (cid:80) i,j =1 b i a ij ( k ni , L k nj ) = s (cid:80) i,j =1 b j a ji ( k ni , L k nj ) and M ij = b i a ij + b j a ji − b i b j are used. Note that L canbe denoted as L = A ∗ A , where A is a linear operator and A ∗ is the adjoint operator of A . Since M ispositive semi-definite, we have s (cid:88) i,j =1 M ij ( k ni , L k nj ) = s (cid:88) i,j =1 M ij ( A k ni , A k nj ) ≥ . (3.10)Combining eqs. (3.9) and (3.10) leads to12 (Φ n +1 , L Φ n +1 ) −
12 (Φ n , L Φ n ) ≤ ∆ t s (cid:88) i =1 b i ( k ni , L Φ ni ) . (3.11)Similarly, we have (cid:107) q n +1 (cid:107) − (cid:107) q n (cid:107) ≤ t s (cid:88) i =1 b i ( l ni , Q ni ) . (3.12)Adding (3.11) and (3.12) and noticing that l ni = (cid:16) ∂g∂ Φ (cid:17) n, ∗ i · k ni + (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i · ∇ k ni , we obtain F n +1 − F n ≤ ∆ t s (cid:88) i =1 b i (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17) , k ni (cid:19) . (3.13)Substituting k ni = G (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19) in (3.13) and noticing the negativesemi-definite property of G and b i ≥ , ∀ i , we arrive at F n +1 − F n ≤ . This completes the proof.
Remark . Note that the Gauss method is a special kind of algebraically stable RK method, whose RKcoefficients satisfy M = . Thus the Gauss method preserves the discrete energy dissipation law F n +1 − F n = ∆ t s (cid:88) i =1 b i (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17) , k ni (cid:19) ≤ . (3.14) Remark . After appropriate spatial discretization that satisfies the discrete integration-by-parts formula(see [16, 17] for details), the algebraically stable LEQRK scheme naturally leads to a fully discrete energystable scheme. In this paper, we employ the Fourier pseudo-spectral method for spatial discretization. We6mit the details here due to space limitation. Interested readers are referred to our ealier work [16, 19] fordetails.Next, we discuss the solvability of the resulting fully discrete scheme.
Theorem 3.2.
If RK coefficient matrix A = ( a ij ) is positive semi-definite and mobility operator G satisfies G = −B ∗ B , the fully discrete scheme derived by applying the Fourier pseudo-spectral method to Scheme3.1 is uniquely solvable.Proof.
Without loss of generality, we still use the notations G , L and ∇ to denote the corresponding discreteoperators in the fully discrete scheme. We consider the homogeneous linear equation system of (3.3) Φ ni = ∆ t s (cid:80) j =1 a ij k nj ,Q ni = ∆ t s (cid:80) j =1 a ij l nj ,k ni = G (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19) ,l ni = (cid:16) ∂g∂ Φ (cid:17) n, ∗ i · k ni + (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i · ∇ k ni , i = 1 , · · · , s, (3.15)where Φ ni , Q ni , k ni , l ni are unknown. To prove unique solvability of the fully discrete scheme, we need toprove that the homogeneous linear equation system (3.15) admits only a zero solution.Computing the discrete inner product of the third equation in (3.15) with L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17) and sum over i , we deduce from (3.15)∆ t s (cid:88) i,j =1 a ij ( A k ni , A k nj ) + 2∆ t s (cid:88) i,j =1 a ij ( l ni , l nj ) + s (cid:88) i =1 (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) B (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) = 0 , where G = −B ∗ B and L = A ∗ A are used. Since A = ( a ij ) is positive semi-definite, which implies that thefirst two terms of the above equation are nonnegative, thus we have B (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19) = 0 , ∀ i, (3.16)which leads to G (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19) = 0 , ∀ i. (3.17)Therefore, according to (3.15), we arrive at k ni = 0 , l ni = 0 , Φ ni = 0 , Q ni = 0 , ∀ i. (3.18)This completes the proof. Theorem 3.3.
If diagnally implicit RK coefficients satisfy a ii > , then the fully discrete scheme derivedby applying the Fourier pseudo-spectral method to Scheme 3.1 is uniquely solvable.Proof.
For the diagnally implicit RK (DIRK) scheme, we solve Φ ni , Q ni , k ni , l ni in turn from i = 1 to s .7herefore, we here consider the following homogeneous linear equation system Φ ni = ∆ ta ii k ni ,Q ni = ∆ ta ii l ni ,k ni = G (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19) ,l ni = (cid:16) ∂g∂ Φ (cid:17) n, ∗ i · k ni + (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i · ∇ k ni , (3.19)where Φ ni , Q ni , k ni , l ni are unknown. To prove unique solvability of the fully discrete scheme, we need toprove that homogeneous linear equation system (3.19) admits only a zero solution.Similar to the proof of Theorem 3.2 , we have∆ ta ii (cid:107)A k ni (cid:107) + 2∆ ta ii (cid:107) l ni (cid:107) = (cid:32) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17) , G (cid:18) L Φ ni + 2 Q ni (cid:16) ∂g∂ Φ (cid:17) n, ∗ i − ∇ · (cid:16) Q ni (cid:16) ∂g∂ ∇ Φ (cid:17) n, ∗ i (cid:17)(cid:19)(cid:33) ≤ , which leads to A k ni = 0 , l ni = 0 . (3.20)Combining (3.19) and (3.20), we deduce in turn l ni = 0 , Q ni = 0 , L k ni = 0 , L Φ ni = 0 , k ni = 0 , Φ ni = 0 . (3.21)This completes the proof. Remark . In this paper, we give examples in two 4th-order algebraically stable RK methods, i.e.,Gauss4th and DIRK4th given by the Butcher tables, respectively, − √
36 14 14 − √ + √
36 14 + √
36 1412 12 , σ σ
12 12 − σ σ − σ σ − σ σµ − µ µ , with σ = cos( π/ / √ / µ = 1 / (cid:0) σ − (cid:1) . We note that Gauss4th and DIRK4th satisfy theconditions in Theorem 3.2 and
Theorem 3.3 , respectively. Therefore, after an appropriate spatialdiscretization, the LEQRK schemes equipped with Gauss4th or DIRK4th are uniquely solvable.
Remark . Noticing that Φ mi approximates Φ( t m + c i ∆ t ), we can choose the time nodes t m , t m + c i ∆ t ( m < n ) and t n as the interpolation points to obtain the interpolation polynomial Φ N ( t ). However, toomany interpolation points will cause the interpolation polynomial to be highly oscillating, which may makeΦ N ( t n + c i ∆ t ) an inaccurate extrapolation for Φ( t n + c i ∆ t ). Therefore, we only take t n − , t n − + c i ∆ t, ∀ i and t n as the interpolation points in this paper. For example, for the Gauss4th method, we choosethe interpolation points ( t n − , Φ n − ), ( t n − + c ∆ t, Φ n − ), ( t n − + c ∆ t, Φ n − ), ( t n , Φ n ) and obtain the8orresponding interpolation polynomialΦ N ( t n − + s ∆ t ) = ( s − c )( s − c )( s − − c c Φ n − + s ( s − c )( s − c ( c − c )( c −
1) Φ n − + s ( s − c )( s − c ( c − c )( c −
1) Φ n − + s ( s − c )( s − c )(1 − c )(1 − c ) Φ n , where c = 1 / − √ / c = 1 / √ /
6. Thus we haveΦ N ( t n + c ∆ t ) = (2 √ − n − + (7 √ − n − + (6 − √ n − + (10 − √ n , (3.22)Φ N ( t n + c ∆ t ) = − (2 √ n − + (6 + 5 √ n − − (7 √ n − + (10 + 4 √ n . (3.23)Replacing n with n − n = Φ n − − √ n − + √ n − . (3.24)According to (3.22)-(3.24), we obtainΦ N ( t n + c ∆ t ) = (6 − √ n − + (1 − √ n − + (5 √ − n − , (3.25)Φ N ( t n + c ∆ t ) = (6 + 2 √ n − − (5 √ n − + (1 + 3 √ n − . (3.26)Note that if we take t n − , t n − + c i ∆ t ( i = 1 ,
2) as the interpolation points, we can also derive (3.25)-(3.26), which implies that the LEQRK scheme induced by the Gauss4th method and the interpolations(3.22)-(3.23) or (3.25)-(3.26) may achieve third order accuracy.
To improve the accuracy as well as stability of
Scheme 3.1 , we propose a prediction-correction schememotivated by the works in [14, 19, 29]. Employing the prediction-correction strategy to
Scheme 3.1 , weobtain the following prediction-correction method:
Scheme 3.2 ( s -stage LEQRK-PC Scheme) . Let b i , a ij ( i, j = 1 , · · · , s ) be real numbers and let c i = s (cid:80) j =1 a ij .For given (Φ n , q n ) and Φ N ( t n + c i ∆ t ) , Q N ( t n + c i ∆ t ) , ∀ i , the following intermediate values are first calculatedby the following prediction-correction strategy1. Prediction: we set Φ n, i = Φ N ( t n + c i ∆ t ) , Q n, i = Q N ( t n + c i ∆ t ) . Let M > be a given integer. For m = 0 to M − , we compute Φ n,m +1 i , k n,m +1 i , l n,m +1 i , Q n,m +1 i using Φ n,m +1 i = Φ n + ∆ t s (cid:80) j =1 a ij k n,m +1 j ,k n,m +1 i = G (cid:18) L Φ n,m +1 i + 2 Q n,mi ∂g∂ Φ [Φ n,mi ] − ∇ · (cid:16) Q n,mi ∂g∂ ∇ Φ [Φ n,mi ] (cid:17)(cid:19) ,l n,m +1 i = ∂g∂ Φ [Φ n,m +1 i ] · k n,m +1 i + ∂g∂ ∇ Φ [Φ n,m +1 i ] · ∇ k n,m +1 i ,Q n,m +1 i = q n + ∆ t s (cid:80) j =1 a ij l n,m +1 j , i = 1 , · · · , s. (3.27) If max i (cid:107) Φ n,m +1 i − Φ n,mi (cid:107) ∞ < T OL , we stop the iteration and set Φ n, ∗ i = Φ n,m +1 i ; otherwise, we set Φ n, ∗ i = Φ n,Mi . . Correction: for the predicted Φ n, ∗ i , we compute the intermediate values Φ ni , Q ni , k ni , l ni via Φ ni = Φ n + ∆ t s (cid:80) j =1 a ij k nj ,Q ni = q n + ∆ t s (cid:80) j =1 a ij l nj ,k ni = G (cid:18) L Φ ni + 2 Q ni ∂g∂ Φ [Φ n, ∗ i ] − ∇ · (cid:16) Q ni ∂g∂ ∇ Φ [Φ n, ∗ i ] (cid:17)(cid:19) ,l ni = ∂g∂ Φ [Φ n, ∗ i ] · k ni + ∂g∂ ∇ Φ [Φ n, ∗ i ] · ∇ k ni , i = 1 , · · · , s. (3.28) Then (Φ n +1 , q n +1 ) is updated via Φ n +1 = Φ n + ∆ t s (cid:88) i =1 b i k ni , (3.29) q n +1 = q n + ∆ t s (cid:88) i =1 b i l ni . (3.30) Remark . Note that Q N ( t ) of Scheme 3.2 denotes the interpolation polynomial of q . If we takeΦ N ( t ) = Φ n , Q N ( t ) = q n , then Scheme 3.1 reduces to first order while
Scheme 3.2 with appropriatepredictions can achieve the desired high order. In numerical computations, we will apply
Scheme 3.2 tofigure out the necessary initial information. In addition, linear system (3.27) is constant coefficient andthus can be readily solved by using the fast Fourier transform (FFT).
Remark . If we choose M = 0 , then Scheme 3.2 reduces to
Scheme 3.1 . If M is large enough, theLEQRK-PC scheme approximates the IEQ-RK scheme proposed in [15]. There is no theoretical resulton the choice of iteration step M . From our numerical experience, several iteration steps M ≤ Remark . Similar to
Scheme 3.1 , we can also establish energy stability and solvability for the LEQRK-PC scheme, which is omitted here to save space.
In the previous sections, we present some high-order linear energy stable schemes for general gradientflow models. In this section, we apply the proposed schemes to two benchmark gradient flow models:the Cahn-Hilliard model for binary fluids and the molecular beam epitaxial (MBE) growth model. Forconvenience, the LEQRK schemes equipped with Gauss4th and DIRK4th are abbreviated respectively asLEQGRK and LEQDIRK, while their corresponding LEQRK-PC schemes with the prediction iteration M are denoted by LEQGRK-PC- M and LEQDIRK-PC- M . We consider the Cahn-Hilliard model for immiscible binary fluids given as follows φ t = λ ∆( − ε ∆ φ + φ − φ ) , (4.1)with the double-well bulk energy F = ε (cid:107)∇ φ (cid:107) + 14 (cid:107) φ − (cid:107) , (4.2)10here λ is the mobility parameter and ε controls the interfacial thickness. If we introduce the auxiliaryvariable q = ( φ − − γ ), where γ ≥ F = 12 (cid:16) φ, − ε ∆ φ + γφ (cid:17) + (cid:107) q (cid:107) − γ + 2 γ | Ω | . (4.3)Then the Cahn-Hilliard equation (4.1) is equivalently transform into the following system φ t = λ ∆( − ε ∆ φ + γφ + 2 qφ ) ,q t = φφ t , (4.4)which satisfies the following energy dissipation law d F dt = − λ (cid:13)(cid:13) ∇ ( − ε ∆ φ + γφ + 2 qφ ) (cid:13)(cid:13) ≤ . (4.5)Applying the LEQRK-PC scheme to system (4.4), we obtain Scheme 4.1.
Let b i , a ij ( i, j = 1 , · · · , s ) be real numbers and c i = s (cid:80) j =1 a ij . For given ( φ n , q n ) and Φ N ( t n + c i ∆ t ) , Q N ( t n + c i ∆ t ) , ∀ i , the following intermediate values are first calculated by the following prediction-correction strategy.1. Prediction: we set Φ n, i = Φ N ( t n + c i ∆ t ) , Q n, i = Q N ( t n + c i ∆ t ) and M > as a given positiveinteger. For m = 0 to M − , we compute Φ n,m +1 i , k n,m +1 i , l n,m +1 i , Q n,m +1 i using Φ n,m +1 i = φ n + ∆ t s (cid:80) j =1 a ij k n,m +1 j ,k n,m +1 i = λ ∆ (cid:0) − ε ∆Φ n,m +1 i + γ Φ n,m +1 i + 2 Q n,mi Φ n,mi (cid:1) ,l n,m +1 i = Φ n,m +1 i k n,m +1 i ,Q n,m +1 i = q n + ∆ t s (cid:80) j =1 a ij l n,m +1 j , i = 1 , · · · , s. (4.6) Given an error tolerance
T OL > , if max i (cid:107) Φ n,m +1 i − Φ n,mi (cid:107) ∞ < T OL , we stop the iteration and set Φ n, ∗ i = Φ n,m +1 i ; otherwise, we set Φ n, ∗ i = Φ n,Mi .2. Correction: for the predicted Φ n, ∗ i , we compute the intermediate values Φ ni , Q ni , k ni , l ni via Φ ni = φ n + ∆ t s (cid:80) j =1 a ij k nj ,Q ni = q n + ∆ t s (cid:80) j =1 a ij l nj ,k ni = λ ∆ (cid:0) − ε ∆Φ ni + γ Φ ni + 2 Q ni Φ n, ∗ i (cid:1) ,l ni = Φ n, ∗ i k ni , i = 1 , · · · , s. (4.7)11 hen ( φ n +1 , q n +1 ) is updated via φ n +1 = φ n + ∆ t s (cid:88) i =1 b i k ni , (4.8) q n +1 = q n + ∆ t s (cid:88) i =1 b i l ni . (4.9)First of all, we present the time mesh refinement tests to show the order of accuracy of the proposedschemes. We consider the domain as [0 2 π ] and choose model parameter values λ = 0 . ε = 1 and γ = 1.Note that the analytical solution for the Cahn-Hilliard equation is usually unknown. To better calculatethe errors in time mesh refinement tests, we create an exact solution φ ( x, y, t ) = sin( x ) sin( y ) cos( t ), byadding a corresponding forcing term on the right-hand side of the Cahn-Hilliard equation. Then, we solveit in a 2D spatial domain with periodic boundary conditions. The equation is discretized spatially usingthe Fourier pseudo-spectral method with 128 spatial meshes.The numerical solution of φ at t = 1 is calculated using a set of different numerical schemes with varioustime steps. Both the L and L ∞ errors in the solution are calculated, and the results are summarized inFigure 4.1. We observe that, due to the low-order extrapolation, LEQDIRK only reaches 2nd orderaccuracy, but it can reach its 4th order accuracy with only two prediction iterations. Similarly, due tothe low-order extrapolation, LEQGRK only has 3rd order accuracy, and it can easily reach its 4th orderaccuracy with one prediction iteration. From Figure 4.1, we also see that the LEQDIRK-PC schemewith only three prediction iterations can reach similar accuracy as IEQDIRK proposed in [15], while theLEQGRK-PC scheme only requires two prediction iterations.To further compare the DIRK4th and Gauss4th schemes, we summarize their L and L ∞ errors inthe same plot, as shown in Figure 4.2(a)-(b), respectively. We observe that the Gauss4th scheme reachesits order of accuracy even with a larger time step size. After a few iterations, the DIRK4th scheme alsoreaches its order of accuracy quickly. Also, with the same time step size, the Gauss4th scheme is moreaccurate than the DIRK4th scheme.To further benchmark these two schemes, we conduct several numerical tests. For comparison, we alsoimplement the widely used 2nd order convex splitting scheme (which we refer to as the 2nd-CS scheme inthis paper), φ n +1 − φ n ∆ t = λ ∆ (cid:104) − ε ∆ φ n + + 12 (( φ n ) + ( φ n +1 ) ) φ n + − ( 32 φ n − φ n − ) (cid:105) . (4.10)We emphasis that there is no theoretical proofs for energy dissipation for the 2nd-CS scheme above, thoughit is more accurate than the first-order convex splitting scheme.For the first example, we choose the domain as [0 1] , and parameters λ = 1, (cid:15) = 0 .
01, and γ = 1.Then, we use the initial condition [34] φ ( x, y, t = 0) = 0 . (cid:16) cos(3 x ) cos(4 y ) + (cos(4 x ) cos(3 y )) + cos( x − y ) cos(2 x − y ) (cid:17) . (4.11)This initial profile would drive a fast coarsening dynamics, such that the algorithm would predict ’wrong’dynamics if it is not accurate or robust enough. In this example, we intend to find the maximum possibletime step that one can capture the correct dynamics numerically. Various numerical schemes with differenttime steps are implemented and compared. The numerical results are summarized in Figure 4.3, where thepredicted profile of φ ( x, y ) at t = 0 . t = 6 . × − . For the DIRK4th scheme with 5 prediction iterations, the maximum time step size is12 a) L and L ∞ errors using the DIRK4th scheme(b) L and L ∞ errors using the Gauss4th scheme Figure 4.1: Time step refinement tests with the proposed numerical schemes for the Cahn-Hilliard equation. (a) L error (b) L ∞ error Figure 4.2: L and L ∞ errors using DIRK4th and Gauss4th schemes for the Cahn-Hilliard equation.13pproximately ∆ t = 1 . × − . For the Gauss4th method with five prediction iterations, it is 2 . × − .Notice, the prediction steps could be easily solved with FFT, so the computational cost is negligiblecompared to the correction step.These results indicate that the DIRK4th and Gauss4th schemes are superior over the 2nd order convexsplitting scheme in this simulation. In addition, one should notice that there is no theoretical guarantee ofmonotonic energy decay with the 2nd order convex splitting scheme, and the implementation of the convexsplitting scheme is relatively complicated, as nonlinear equations have to be solved at each time step. Incontrast, the proposed high-order schemes here are linear and easy to implement. Also, they are rathergeneral so that they can be applied to a broad class of gradient flow models.To further confirm these findings, we conduct an additional numerical experiment with random initialconditions. Specifically, we use φ ( x, y, t = 0) = 0 . x, y ) , (4.12)where rand( x, y ) generates random numbers between − In this subsection, we focus on the molecular beam epitaxial growth model with slope selection givenas follows φ t = − λ (cid:16) ε ∆ φ − ∇ · (cid:0) ( |∇ φ | − ∇ φ (cid:1)(cid:17) , (4.13)where the free energy functional is given by F = ε (cid:107) ∆ φ (cid:107) + 14 (cid:13)(cid:13) |∇ φ | − (cid:13)(cid:13) . (4.14)We let q = (cid:0) |∇ φ | − − γ (cid:1) and rewrite the energy functional as F = 12 (cid:16) φ, ε ∆ φ − γ ∆ φ (cid:17) + (cid:107) q (cid:107) − γ + 2 γ | Ω | . (4.15)Using the EQ reformulation, we have the following equivalent system (cid:40) φ t = − λ (cid:16) ε ∆ φ − γ ∆ φ − ∇ · (cid:0) q ∇ φ (cid:1)(cid:17) ,q t = ∇ φ · ∇ φ t , (4.16)with the consistent initial condition (cid:40) φ ( t = 0) = φ ,q ( t = 0) = (cid:0) |∇ φ | − − γ (cid:1) . (4.17)It is readily to show that new system (4.16) obeys the following energy dissipation law d F dt = − λ (cid:13)(cid:13)(cid:13) ε ∆ φ − γ ∆ φ − ∇ · (cid:0) q ∇ φ (cid:1)(cid:13)(cid:13)(cid:13) ≤ . (4.18)Applying the LEQRK-PC scheme for system (4.16), we have the following scheme.14 a) The profile of φ at t = 0 . t = 2 . × − , . × − , . × − with the 2nd-order convex splitting scheme.(b) The profile of φ at t = 0 . t = 5 × − , . × − , . × − with the DIRK4th scheme.(c) The profile of φ at t = 0 . t = 5 × − , . × − , . × − with the Gauss4th scheme. Figure 4.3: A comparison of the three schemes on predicting accurate Cahn-Hilliard dynamics at varioustime step sizes. The figures show the numerical results of φ at time t = 0 . a) The profile of φ at t = 0 . t = 2 . × − , . × − , . × − with the 2nd order convex splitting method(b) The profile of φ at t = 0 . t = 5 × − , . × − , . × − with the DIRK4th method.(c) The profile of φ at t = 0 . t = 5 × − , . × − , . × − with the Gauss4th method. Figure 4.4: A comparison of the three schemes on predicting accurate Cahn-Hilliard dynamics with randominitial conditions. These figures show the numerical results of φ at time t = 0 . cheme 4.2. Let b i , a ij ( i, j = 1 , · · · , s ) be real numbers and let c i = s (cid:80) j =1 a ij . For given ( φ n , q n ) and Φ N ( t n + c i ∆ t ) , Q N ( t n + c i ∆ t ) , ∀ i , the following intermediate values are first calculated by the followingprediction-correction strategy1. Prediction: we set Φ n, i = Φ N ( t n + c i ∆ t ) , Q n, i = Q N ( t n + c i ∆ t ) and M > as a given integer. For m = 0 to M − , we compute Φ n,m +1 i , k n,m +1 i , l n,m +1 i , Q n,m +1 i using Φ n,m +1 i = φ n + ∆ t s (cid:80) j =1 a ij k n,m +1 j ,k n,m +1 i = − λ (cid:16) ε ∆ Φ n,m +1 i − γ ∆Φ n,m +1 i − ∇ · (cid:0) Q n,mi ∇ Φ n,mi (cid:1)(cid:17) ,l n,m +1 i = ∇ Φ n,m +1 i · ∇ k n,m +1 i ,Q n,m +1 i = q n + ∆ t s (cid:80) j =1 a ij l n,m +1 j , i = 1 , · · · , s. (4.19) Given the error tolerance
T OL > , if max i (cid:107) Φ n,m +1 i − Φ n,mi (cid:107) ∞ < T OL , we stop the iteration and set Φ n, ∗ i = Φ n,m +1 i ; otherwise, we set Φ n, ∗ i = Φ n,Mi .2. Correction: for the predicted Φ n, ∗ i , we compute the intermediate values Φ ni , Q ni , k ni , l ni via Φ ni = φ n + ∆ t s (cid:80) j =1 a ij k nj ,Q ni = q n + ∆ t s (cid:80) j =1 a ij l nj ,k ni = − λ (cid:16) ε ∆ Φ ni − γ ∆Φ ni − ∇ · (cid:0) Q ni ∇ Φ n, ∗ i (cid:1)(cid:17) ,l ni = ∇ Φ n, ∗ i · ∇ k ni , i = 1 , · · · , s. (4.20) Then ( φ n +1 , q n +1 ) is updated via φ n +1 = φ n + ∆ t s (cid:88) i =1 b i k ni , (4.21) q n +1 = q n + ∆ t s (cid:88) i =1 b i l ni . (4.22)We apply the proposed arbitrarily high order schemes to solve MBE model (4.13). We repeat the timestep refinement test first. Here we use domain [0 2 π ] and choose parameters λ = 0 . γ = 1 and ε = 1.By adding the proper force term on the right-hand side of the equation, we create the real solution φ ( x, y, t ) = sin( x ) sin( y ) cos( t ) . (4.23)for the MBE model. Then we solve the modified model in the domain with a periodic boundary usingthe pseudo-spectral method for spatial discretization on 128 meshes. The L errors and L ∞ errors usingdifferent schemes and various time steps are summarized in Figure 4.5. Here we observe similar results,i.e., the DIRK4th scheme reaches 2nd order accuracy without prediction, but obtain 4th order accuracywith only two iteration steps for both the L and L ∞ norms. This is due to the low-order approximationfor extrapolating the explicit terms. Analogously, the Gauss4th scheme is 3rd order accurate without anyprediction steps and reaches 4th order accuracy with one iteration step.17 a) L and L ∞ errors using the DIRK4th scheme(b) L and L ∞ errors using the Gauss4th scheme Figure 4.5: Time step refinement tests with the proposed numerical schemes for the MBE model.18 a) L errors (b) L ∞ errors Figure 4.6: L and L ∞ errors using the DIRK4th method and the Gauss4th method for the MBE model.To compare the accuracy of the DIRK4th scheme with that of the Gauss4th scheme in solving themolecular beam epitaxy (MBE) model, we summarize their L and L ∞ errors in the same figure as shownin 4.6. We observe that the Gauss4th method has smaller errors than the DIRK4th scheme if using thesame time step sizes.Next, we use the proposed DIRK4th and Gauss4th schemes to solve two benchmark problems associatedto the MBE model (4.13). As before, we introduce the 2nd-order convex splitting scheme φ n +1 − φ n δt = − λ (cid:16) ε ∆ φ n + − ∇ · (cid:0) ( |∇ φ n +1 | + |∇ φ n | ) ∇ φ n + (cid:1) + ∆( 32 φ n − φ n − ) (cid:17) , (4.24)which will be used for comparison with the proposed linear high-order schemes.Following [34], we choose the domain as [0 2 π ] , parameters λ = 1, ε = 0 .
1, and γ = 1. We solve theMBE model in a periodic domain using the pseudo-spectral method with 128 meshes. All the numericalschemes (i.e., the DIRK4th, Gauss4th, and 2nd CS scheme) are implemented. Five prediction iterationsare used for both the DIRK4th and Gauss4th scheme. The energy from t = 0 to t = 15 are calculatedwith different time steps and the results are summarized in Figure 4.7. We observe the maximum timesteps to obtain accurate solutions are ∆ t = 0 . t = 0 . t = 0 . . and ε = 0 .
03. The rest parameters are the same as inprevious examples. We use 256 meshes. It is known that the MBE coarsening dynamics follows a powerlaw, where the energy decreases as O ( t − ), and the roughness increases as O ( t ) [34]. The numericalresults are summarized in Figure 4.8, showing a strong agreement with the expected power law.The profile of φ and ∆ φ at different times are summarized in Figure 4.9 and 4.10, respectively. Theseprofiles look qualitatively similar to the reported results. These results strongly support our claim that19 a) Convex splitting scheme (b) LEQDIRK-PC-5 (c) LEQGRK-PC-5 Figure 4.7: Numerical results of energy evolution for the MBE model using different schemes with varioustime steps. (a) Energy (b) Rougness
Figure 4.8: The numerical results show proper power law dynamics for the decreasing energy as O ( t − / )and increasing roughness as O ( t / ). 20 a) φ at t = 0 (b) φ at t = 5 (c) φ at t = 10(d) φ at t = 20 (e) φ at t = 30 (f) φ at t = 45 Figure 4.9: The profile of φ at different times.the general arbitrarily high order linear schemes can be applied to predict accurate dynamics for the MBEmodel. In this paper, we present a new paradigm for developing arbitrarily high order, fully discrete numericalalgorithms. These newly proposed algorithms have several advantageous properties: (1) the schemesare all linear such that they are easy to implement and computationally efficient; (2) the schemes areunconditionally energy stable and uniquely solvable such that large time steps can be used in some long-time dynamical simulations; (3) the schemes can reach arbitrarily high-order of accuracy spatial-temporallysuch that relatively large meshes can guarantee the desired accuracy of numerical solutions; (4) the schemesdo not depend on the specific expression of the free energy explicitly such that it can be readily appliedto a large class of general gradient flow models. The proofs for energy stability and uniquely solvabilityare given, and numerical tests with benchmark problems are shown to illustrate the effectiveness of theproposed schemes.
Acknowledgments
Yuezheng Gong’s work is partially supported by the Natural Science Foundation of Jiangsu Province(Grant No. BK20180413) and the National Natural Science Foundation of China (Grant No. 11801269).Jia Zhao’s work is partially supported by National Science Foundation under grant number NSF DMS-1816783; and National Institutes of Health under grant number R15-GM132877. Jia Zhao would also like21 a) ∆ φ at t = 0 (b) ∆ φ at t = 5 (c) ∆ φ at t = 10(d) ∆ φ at t = 20 (e) ∆ φ at t = 30 (f) ∆ φ at t = 45 Figure 4.10: The profile of ∆ φ at different times.to acknowledge NVIDIA Corporation for their donation of a Quadro P6000 GPU for conducting some ofthe numerical simulations in this paper. Qi Wang’s work is partially supported by DMS-1815921, OIA-1655740 and a GEAR award from SC EPSCoR/IDeA Program. Supports by NSFC awards References [1] J. W. Barrett, J. F. Blowey, and H. Garcke. On fully practical finite element approximations ofdegenerate cahn-hilliard systems.
ESAIM: M2AN , 35:713–748, 2001.[2] K. Burrage and J. C. Butcher. Stability criteria for implicit Runge-Kutta methods.
SIAM Journal onNumerical Analysis , 16(1):46–57, 1979.[3] J. W. Cahn and S. M. Allen. A microscopic theory for domain wall motion and its experimentalvarification in fe-al alloy domain growth kinetics.
J. Phys. Colloque , C7:C7–51, 1977.[4] J. W. Cahn and J. E. Hilliard. Free energy of a nonuniform system. I. interfacial free energy.
Journalof Chemical Physics , 28:258–267, 1958.[5] L. Chen, J. Zhao, and X. Yang. Regularized linear schemes for the molecular beam epitaxy modelwith slope selection.
Applied Numerical Mathematics , 128:138–156, 2018.[6] W. Chen, S. Conde, C. Wang, X. Wang, and S. Wise. A linear energy stable scheme for a thin filmmodel without slope selection.
Journal of Scientific Computing , 52:546–562, 2012.227] Q. Cheng, X. Yang, and J. Shen. Efficient and accurate numerical schemes for a hydrodynamicallycoupled phase field diblock copolymer model.
Journal of Computational Physics , 341:44–60, 2017.[8] S. Clarke and D. D. Vvedensky. Origin of reflection high-energy electron-diffraction intensity os-cillations during molecular-beam epitaxy: A computational modeling approach.
Phys. Rev. Lett. ,58:2235–2238, 1987.[9] S. Dong, Z. Yang, and L. Lin. A family of second-order energy-stable schemes for cahn-hilliard typeequations.
Journal of Computational Physics , 383:24–54, 2019.[10] M. Elsey and B. Wirth. A simple and efficient scheme for phase field crystal simulation.
ESAIM:M2AN , 47:1413–1432, 2013.[11] D. Eyre. Unconditionally gradient stable time marching the Cahn-Hilliard equation.
Computationaland mathematical models of microstructural evolution (San Francisco, CA, 1998) , 529:39–46, 1998.[12] W. Feng, C. Wang, S. Wise, and Z. Zhang. A second-order energy stable backward differentiationformula method for the epitaxial thin film equation with slope selection.
Numerical Methods for PartialDifferential Equations , 34(6):1975–2007, 2018.[13] D. Furihata and T. Matsuo.
Discrete variational derivative method. A structure-preserving numericalmethod for partial differential equations . CRC Press, 2011.[14] K. Glasner and S. Orizaga. Improving the accuracy of convexity splitting methods for gradient flowequations.
Journal of Computational Physics , 315:52–64, 2016.[15] Y. Gong and J. Zhao. Energy-stable Runge-Kutta schemes for gradient flow models using the energyquadratization approach.
Applied Mathematics Letters , 94:224–231, 2019.[16] Y. Gong, J. Zhao, and Q. Wang. Linear second order in time energy stable schemes for hydrody-namic models of binary mixtures based on a spatially pseudospectral approximation.
Advances inComputational Mathematics , 44:1573–1600, 2018.[17] Y. Gong, J. Zhao, and Q. Wang. Second order fully discrete energy stable methods on staggered gridsfor hydrodynamic phase field models of binary viscous fluids.
SIAM Journal on Scientific Computing ,40(2):B528–B553, 2018.[18] Y. Gong, J. Zhao, and Q. Wang. Arbitrarily high-order unconditionally energy stable sav schemes forgradient flow models.
ArXiv , pages 1–22, 2019.[19] Y. Gong, J. Zhao, and Q. Wang. Arbitrary high order energy stable schemes for gradient flow modelsbased on energy quadratization.
SIAM Journal on Scientific Computing , page accepted, 2019.[20] F. Guillen-Gonzalez and G. Tierra. Second order schemes and time-step adaptivity for Allen-Cahnand Cahn-Hilliard models.
Computers Mathematics with Applications , 68(8):821–846, 2014.[21] Z. Guo, P. Lin, J. Lowengrub, and S. Wise. Mass conservative and energy stable finite difference meth-ods for the quasi-incompressible navier-stokes-cahn-hilliard system: primitive variable and projection-type schemes.
Computer Methods in Applied Mechanics and Engineering , 326:144–174, 2017.[22] D. Hou, M. Azaiez, and C. Xu. A variant of scalar auxiliary variable approaches for gradient flows.
Journal of Computational Physics , 395:307–332, 2019.2323] W. Hu and M. lai. Unconditionally energy stable immersed boundary method with application tovesicle dyanmics.
East Asian Journal on Applied Mathematics , 3(3):247–262, 2013.[24] J. Kim, K. Kang, and J. Lowengrub. Conservative multigrid methods for ternary cahn-hilliard systems.
Commun. Math. Sci. , 2:53–77, 2004.[25] Z. Liu and X. Li. Efficient modified techniques of invariant energy quadratization approach for gradientflows.
Applied Math. Letters , 98:206–214, 2019.[26] F. Luo, H. Xie, M. Xie, and F. Xu. Adaptive time-stepping algorithms for molecular beam epitaxy:based on energy or roughness.
Applied Mathematics Letters , 99:105991, 2020.[27] Z. Qiao, Z. Sun, and Z. Zhang. Stability and convergence of second-order schemes for the nonlinearepitaxial growth model without slope selection.
Mathematics of Computation , 84(292):653–674, 2015.[28] Z. Qiao, Z. Zhang, and Tao Tang. An adaptive time-stepping strategy for the molecular beam epitaxymodels.
SIAM Journal on Scientific Computing , 33(3):1395–1414, 2011.[29] J. Shen and J. Xu. Stabilized predictor-corrector schemes for gradient flows with strong anisotropicfree energy.
Communications in Computational Physics , 24(3):635–654, 2018.[30] J. Shen, J. Xu, and J. Yang. The scalar auxiliary variable (sav) approach for gradient flows.
Journalof Computational Physics , 353:407–416, 2018.[31] J. Shen, J. Xu, and J. Yang. A new class of efficient and robust energy stable schemes for gradientflows.
SIAM Review , 61(3):474–506, 2019.[32] J. Shen and X. Yang. Numerical approximations of Allen-Cahn and Cahn-Hilliard equations.
Disc.Conti. Dyn. Sys.-A , 28:1669–1691, 2010.[33] C. Theng, I. Chern, and M. Lai. Simulating binary fluid-surfactant dynamics by a phase field model.
Discrete and Continuous Dynamical Systems B , 17(4):1289–1307, 2012.[34] C. Wang, X. Wang, and S. Wise. Unconditionally stable schemes for equations of thin film epitaxy.
Discrete and Continuous Dynamic Systems , 28(1):405–423, 2010.[35] S.-L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun, and G.B. McFadden.Thermodynamically-consistent phase-field models for solidification.
Physica D , 69:189–200, 1993.[36] S. Wise, J. Kim, and J. Lowengrub. Solving the regularized strongly anisotropic cahn-hilliard equationby an adaptive nonlinear multigrid method.
Journal of Computational Physics , 226(1):414–446, 2007.[37] S. Wise, C. Wang, and J. S. Lowengrub. An energy-stable and convergent finite-difference scheme forthe phase field crystal equation.
SIAM Journal of Numerical Analysis , 47(3):2269–2288, 2009.[38] Z. Xu, X. Yang, H. Zhang, and Z. Xie. Efficient and linear schemes for anisotropic Cahn-Hilliard modelusing the stabilized-invariant energy quadratization (S-IEQ) approach.
Computer Physics Communi-cations , 238:36–49, 2019.[39] X. Yang. Linear, first and second-order, unconditionally energy stable numerical schemes for the phasefield model of homopolymer blends.
Journal of Computational Physics , 327:294–316, 2016.[40] X. Yang and L. Ju. Efficient linear schemes with unconditional energy stability for the phase fieldelastic bending energy model.
Journal of Computational Physics , 315:691–712, 2017.2441] X. Yang, J. Zhao, and Q. Wang. Numerical approximations for the molecular beam epitaxial growthmodel based on the invariant energy quadratization method.
J. Comput. Phys. , 333:104–127, 2017.[42] X. Ye. The Legendre collocation method for the cahn-hilliard equation.
Journal of Computationaland Applied Mathematics , 150:87–108, 2003.[43] Z. Zhang, Y. Ma, and Z. Qiao. An adaptive time-stepping strategy for solving the phase field crystalmodel.
Journal of Computational Physics , 249:204–215, 2013.[44] J. Zhao, X. Yang, Y. Gong, X. Zhao, J. Li, X. Yang, and Q. Wang. A general strategy for numericalapproximations of thermodynamically consistent nonequilibrium models-part I: Thermodynamicalsystems.