A linearly second-order energy stable scheme for the phase field crystal model

Applied Numerical Mathematics 140 (2019) 134–164

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A linearly second-order energy stable scheme for the phase field crystal model Shuaichao Pei, Yanren Hou ∗ , Bo You School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China

a r t i c l e

i n f o

Article history: Received 27 October 2018 Received in revised form 29 January 2019 Accepted 30 January 2019 Available online 14 February 2019 Keywords: Phase field crystal model Energy stability Linear scheme Second-order accuracy

a b s t r a c t In this paper, we propose a linear, unconditionally energy stable and second-order (in time) numerical scheme based on a convex splitting scheme and the semi-implicit spectral deferred correction (SISDC) method for the phase field crystal equation. The convex splitting scheme we use is linear, uniquely solvable and unconditionally energy stable but is of first-order, so we take the SISDC method to improve the rate of convergence. The resulted scheme inherits the advantages of the convex splitting scheme and thus leads to linear equations at each time step, which is easy to implement. We also prove that the scheme is unconditionally weak energy stable and of second-order accuracy in time. Numerical experiments are presented to validate the accuracy, efficiency and energy stability of the proposed numerical strategy. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The phase field crystal (PFC) model proposed by Elder et al. [8,9] has shown great versatility to describe the phenomena of crystal growth on the atomic length and diffusive time scales. In the model, a phase field variable, which represents the concentration field of a coarse-grained temporal average of the density of atoms, is introduced to describe the phase transition from the liquid phase to the crystal phase. The model can account for elastic and plastic deformations of the lattice, dislocations, grain boundaries, and many other observable phenomena. The reference [22] describes the variety of micro-structures that can be modeled using the PFC approach. The PFC equation is a sixth-order nonlinear partial differential equation, which satisfies an energy law (see from section 2). It is desirable to design efficient numerical schemes that preserve the energy law at the discrete level. There have been lots of numerical analysis works on the energy stable schemes for the PFC model recently. In [6,18,31], the authors proposed first-order and second-order numerical schemes based on the convex splitting method. Both schemes are unconditionally stable, uniquely solvable, and an efficient nonlinear multi-grid method is designed to solve the resulted nonlinear systems. In [28], the authors developed numerical strategies based on a new convex splitting for the PFC model. The strategies are unconditionally stable, uniquely solvable, and a nonlinear iterative method is applied. Glasner et al. [13] introduced a linear, unconditionally energy stable algorithm based on a linear convex splitting. Yang and Han [33] presented a series of linear, unconditionally energy stable numerical schemes based on the invariant energy quadratization approach (see [34]). In [25,26], Shen et al. proposed the scalar auxiliary variable approach for a large class of gradient flows, which leads to linear, energy stable schemes. Vignal et al. [30] developed a nonlinear, second-order time accurate, and energy

*

Corresponding author. E-mail addresses: [email protected] (S. Pei), [email protected] (Y. Hou), [email protected] (B. You).

https://doi.org/10.1016/j.apnum.2019.01.017 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

135

stable strategy by using the Crank–Nicolson method and the additional stabilization approach. Gomez et al. [14] introduced a nonlinear, second-order time accurate, and unconditionally energy stable algorithm with the modified Crank–Nicolson method. In [35], the authors gave nonlinear, unconditionally energy stable schemes with an adaptive time step strategy. Guo et al. [16] presented some unconditionally energy stable schemes using the local discontinuous Galerkin method. The goal of this paper is to develop an efficient and energy stable numerical scheme. Inspired by the work in [12], we first consider a linear implicit-explicit scheme by using the convex splitting method in [13], and we prove rigorously that the scheme is unconditionally energy stable and is of first-order accuracy in time. Then we apply the semi-implicit spectral deferred correction (SISDC) method to achieve a second-order scheme. Rigorously proof is given to demonstrate the unconditionally weak energy stability and the second-order convergence in time. The reasons for us to employ the SISDC method are as follows: iteration loops can improve the order of accuracy in a flexible and simple way, for instance, we only need to iterate once to get the desirable second-order accurate approximate solution; since the resulted scheme is a one-step method, and all we need is the information from the first-order scheme, it is easy to analyze and implement. Based on the convex splitting and the SISDC method, we get first-order and second-order fully discrete numerical schemes, both linear. So we just need to solve two linear constant coefficient equations at each time step, which is easy to implement and efficient. Some numerical experiments are presented to validate the accuracy and energy stability of our scheme. The rest of this paper is organized as follows. In section 2, we introduce the PFC equation, the closely-related Swift– Hohenberg (SH) equation and their corresponding energy laws. In section 3, the first-order convex splitting scheme is presented. We prove the unconditional energy stability and give the error estimate in the time discrete case. In section 4, we give the second-order scheme, which follows a brief review on the SISDC method. Rigorous proof is given to demonstrate the scheme is unconditionally weak energy stable and of second-order accuracy in time. The fully discrete schemes of the first-order and the second-order are established in section 5. We present some numerical experiments to demonstrate the accuracy and energy stability of our scheme in section 6. Some concluding remarks are presented in section 7. 2. Phase field crystal model Herein we consider a dimensionless energy of the form [9,29]

 

E (φ) =

1 4

4

φ +

1−ε 2

2

2

1



2

φ − |∇φ| + (φ) 2

dx,

(2.1)



where  is a bounded domain in Rd (d = 2, 3) with smooth boundary ∂ , φ :  → R is the density field, constant. Taking the variational derivative of (2.1) in L 2 and H −1 , we get two common types of gradient dynamics:

ε ∈ (0, 1) is a

• non-conserved dynamics

φt = −G (φ)μ,

(2.2)

• conserved dynamics

φt = ∇ · (G (φ)∇ μ) , where G (φ) > 0 is a mobility,

μ :=

(2.3)

μ is the chemical potential defined as

δ E (φ) = φ 3 + (1 − ε )φ + 2φ + 2 φ, δφ

δ E (φ)

and δφ denotes the variational derivative with respect to φ . We impose either the periodic boundary conditions or the homogeneous Neumann boundary conditions. Equation (2.2) is the Swift–Hohenberg (SH) equation (see [29]) and is of fourth order in space. Equation (2.3) is the phase field crystal(PFC) equation and is a sixth-order equation. d 2 The PFC system is mass-conservative in the sense that dt  φ dx = 0. Moreover, by taking the L inner product of (2.3) with μ, one finds that the system satisfies an energy law

d dt



E (φ) = − G (φ)∇ μ2L 2 ≤ 0.

(2.4)

One could find that the SH equation (2.2) also satisfies an energy law

d dt



E (φ) = − G (φ)μ2L 2 ≤ 0,

but (2.2) does not conserve mass. In this paper, we primarily concern with the PFC equation (2.3), and most of the theoretical results and numerical algorithms can be applied to the SH equation (2.2) as well.

136

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

In the PFC system, the density is relatively homogeneous in the liquid phase and spatially periodic in the solid phase (cf. [8,9]), which implies that the density field φ is bounded. Besides, the numerical solutions always stay well bounded in practical numerical simulations. So, just like a common practice (cf. [2,5,19,27]) of considering the Cahn–Hilliard equation with a truncated potential, we replace the function 14 φ 4 in E (φ) by

⎧ 2 3M 3 2 3 4 ⎪ ⎨ 2 φ − 2M φ + 4 M , φ > M , 1 4 F (φ) = 4 φ , φ ∈ [− M , M ], ⎪ ⎩ 3M 2 2 3 3 4 φ + 2M φ + M , φ < −M , 2 4

where M is a positive number to be chosen. Its derivative is

⎧ 2 3 ⎪ ⎨3M φ − 2M , φ > M , 3 f (φ) = φ , φ ∈ [− M , M ], ⎪ ⎩ 2 3M φ + 2M 3 , φ < − M .

Remark 2.1. With the replacement of

1 4 φ 4

by F (φ), we get

max | f  (φ)| ≤ 3M 2 .

(2.5)

φ∈R

Then, we can obtain the discrete energy law and derive a H 2 estimate for the numerical solution (see section 3), which implies an L ∞ bound for the numerical solution. For convenience, we assume that G (φ) = 1. The PFC model can be rewritten as follows:

φt = μ,

in  × (0, T ],

(2.6a)

μ = f (φ) + (1 − ε)φ + 2φ +  φ,

in  × (0, T ],

(2.6b)

φ(x, 0) = φ0 (x),

in ,

(2.6c)

2

2

∂n φ|∂  = 0, ∂n (φ)|∂  = 0, ∂n ( φ)|∂  = 0,

(2.6d)

where n is the unit outward normal vector on the boundary ∂ . In the following, we use (·, ·) and  ·  to denote the standard L 2 inner product in L 2 () and the corresponding induced L 2 norm. 3. The first-order convex splitting scheme 3.1. A linear convex splitting scheme The convex splitting method proposed by Eyre [11] has been applied to various gradient flows (see, e.g., [1,15,24,31]). The idea of the energy convex splitting is to split the free energy (2.1) into the difference of two convex functionals E c (φ) and E e (φ). The ways of splitting are not unique, they can be either linear or nonlinear. Here we want to design linear energy stable schemes. So we take the energy splitting formulation (see [13]) as follows:

E (φ) = E c (φ) − E e (φ),

 

E c (φ) =

1−ε+β 2

1



 

φ 2 + (φ)2 dx, E e (φ) = 2



β 2 φ + |∇φ|2 − F (φ) dx, 2

(3.1a) (3.1b)



where β is a positive number to be chosen. Lemma 3.1. If β ≥ 3M 2 , then the splitting (3.1) satisfies the standard convex splitting conditions, i.e., E c (φ) and E e (φ) are convex functionals. Proof. For E c (φ),

E c (φ + sψ) = E c (φ) + s





(1 + + β)φ +  φ ψ dx +



It shows

2

s2

 (1 − + β)ψ 2 + (ψ)2 dx + O (s3 ).

2 

(3.2)

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

dE c ds

(φ + sψ)|s=0 =

137



 (1 + + β)φ + 2 φ ψ dx,



so we get

δ E c (φ) = (1 + + β)φ + 2 φ. δφ

(3.3)

The equality (3.2) also shows

d2 E c ds2

 (1 − + β)ψ 2 + (ψ)2 dx ≥ 0,

(φ + sψ)|s=0 = 

which proves that E c (φ) is convex. For E e (φ),

 (βφ − 2φ − f (φ)) ψ dx +

E e (φ + sψ) = E e (φ) + s

s2 2



It shows

dE e ds

  βψ 2 − f  (φ)ψ 2 + 2|∇ψ|2 dx + O (s3 ).

(3.4)



 (βφ − 2φ − f (φ)) ψ dx,

(φ + sψ)|s=0 = 

so we get

δ E e (φ) = βφ − 2φ − f (φ). δφ

(3.5)

Thanks to (2.5) and β ≥ 3M 2 , we have

d2 E e ds2

(φ + sψ)|s=0 =



 βψ 2 − f  (φ)ψ 2 + 2|∇ψ|2 dx ≥ 0.



So E e (φ) is convex.

2

Since the goal of this work is to improve the rate of convergence in time, we first present related schemes at the discrete-time, continuous-space level. Let N be a positive integer and 0 = t 0 < t 1 < · · · < t N = T be a uniform partition of [0, T ], with the time step length τ = NT . We denote the approximate solution as φn , μn , w n at tn , and take the initial data φ0 = φ0 (x). Based on the convex splitting (3.1), we give a first-order numerical scheme. For 0 ≤ n ≤ N − 1, given φn , we compute φn+1 , μn+1 , w n+1 as follows:

φn+1 − φn

τ

= μn+1 ,

(3.6a)

μn+1 = (1 − ε + β)φn+1 +  w n+1 − (βφn − 2φn − f (φn )) ,

(3.6b)

w n+1 = φn+1 ,

(3.6c)

∂n φn+1 |∂  = 0, ∂n (φn+1 )|∂  = 0, ∂n (2 φn+1 )|∂  = 0.

(3.6d)

Remark 3.1. The scheme (3.6) can also be viewed as a first-order stabilized scheme (cf. [27,32]). We view it a convex splitting scheme here to easily obtain its uniquely solvability. The scheme (3.6) is mass-conservative, uniquely solvable and unconditionally energy stable. Indeed, we have the following theorems. Theorem 3.1. The scheme (3.6) is mass-conservative and uniquely solvable for any τ > 0. Proof. The mass conservation follows from (φn+1 − φn , 1) = τ (μ, 1) = −(∇ μ, ∇ 1) = 0. The proof of the uniquely solvability is based on convexity arguments, one can find a detail proof in [28]. 2

138

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

Theorem 3.2. If β ≥ 3M 2 , the scheme (3.6) is unconditionally energy stable, i.e.,

E (φn+1 ) ≤ E (φn ), 0 ≤ n ≤ N − 1, ∀τ > 0.

(3.7)

Proof. Taking the inner product of (3.6a), (3.6b) with

μn+1 , φn+1 − φn , respectively, we obtain

(φn+1 − φn , μn+1 ) = −τ (∇ μn+1 , ∇ μn+1 ), (μn+1 , φn+1 − φn ) =(1 − ε )(φn+1 , φn+1 − φn ) + β(φn+1 − φn , φn+1 − φn ) + (φn+1 , φn+1 − φn ) + 2(φn , φn+1 − φn ) + ( f (φn ), φn+1 − φn ). Furthermore,

 φn+1 2 − φn 2 + φn+1 − φn 2 , 2  1 φn+1 2 − φn 2 + φn+1 − φn 2 , (φn+1 , φn+1 − φn ) = (1 − ε )(φn+1 , φn+1 − φn ) =

1−ε

2 2 (φn , φn+1 − φn ) = −2 (∇φn , ∇φn+1 − ∇φn )

= ∇φn 2 − ∇φn+1 2 + ∇φn+1 − ∇φn 2 . For the term ( f (φn ), φn+1 − φn ), we use the Taylor expansion for F (φ):

F (φn+1 ) = F (φn ) + f (φn )(φn+1 − φn ) +

f  (ξn ) 2

(φn+1 − φn )2 ,

Then,

 ( f (φn ), φn+1 − φn ) = ( F (φn+1 ), 1) − ( F (φn ), 1) −

f  (ξn ) 2

ξn is between φn and φn+1 .

(φn+1 − φn )2 , 1 .

Combining the above relations together, we arrive at

−τ ∇ μn+1 2 = E (φn+1 ) − E (φn ) +

1 − ε + 2β − f  (ξn ) 2

φn+1 − φn 2

(3.8)

1

+ φn+1 − φn 2 + ∇φn+1 − ∇φn 2 . 2

Thanks to (2.5) and β ≥ 3M 2 , we have

1 − ε + 2β − f  (ξn ) 2

φn+1 − φn 2 ≥ 0,

the equality (3.8) implies E (φn+1 ) ≤ E (φn ), ∀τ > 0. This completes the proof.

2

Thanks to the discrete energy law (3.7), we can get an uniform-in-time H 2 estimate of the numerical solution as follows. Lemma 3.2. Suppose that the initial data is sufficiently regular that E (φ0 ) ≤ C 1 , for some constant C 1 , then

φn  H 2 ≤ C 2 , 0 ≤ n ≤ N ,

(3.9)

where C 2 is a constant depending on C 1 and . Proof. From the discrete energy law (3.7), we obtain E (φn ) ≤ E (φ0 ) ≤ C 1 , 0 ≤ n ≤ N. If φn ∈ [− M , M ], we have

1 1−ε 1 E (φn ) = φn 4L 4 + φn 2 − ∇φn 2 + φn 2 4 2 2  1 1 5 2 2 2 2

=

4

Using the inequality ∇φn  = −(φn , φn ) ≤ 2

E (φn ) ≥ Since φn  ≤ 2

C0 4

+ φn  − ∇φn 2 +

φn  + ∇φn  + φn 

φn 2H 2 −

φn 2L 4

1 2

21 + 8ε

· || ≤

16 1

4

4

5 φn 2 4

1 φn 2 , 5

+

1

φn 2 + φn 4L 4 .

φn 4L 4 21+8ε

4

+

21+8ε ||, 4

we obtain

we get

1 − 2ε 4

1

φn 2 + φn 4L 4 . 4

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

E (φn ) ≥

C0 4

φn 2H 2 +

which implies φn 2H 2 ≤ If |φn | > M, we have

E (φn ) =

=

3M 2 2

4C 1 C0

3 16

φn 4L 4 −

(21+8ε )2

+

16C o

(21 + 8ε )2 64

||,

||.

 φn 2 ± 2M 3

139

φn dx +

3 4

1−ε

M 4 || +

2

1

φn 2 − ∇φn 2 + φn 2 2



1

 15 φn 2 + ∇φn 2 + φn 2 + φn 2 − 32 32  48M 2 + 16(1 − ε ) − 1 + φn 2 ± 2M 3 φ0 dx + 32

33

∇φn 2

32 3 4 M ||. 4



Using the inequality ∇φn 2 = −(φn , φn ) ≤

E (φn ) ≥

C0 32

φn 2H 2 +

11 φn 2 20

5 φn 2 , 11

+

960M 2 + 320(1 − ε ) − 383 640

we get

 φn 2 ± 2M 3

φ0 dx +

3 4

M 4 ||.



Letting M ≥ 1, we have 960M + 320(1 − ε ) − 383 ≥ 0. Then (3.9) follows. 2 2

φn 2H 2

≤

32C 1 C0

∓

64M 3 C0



24

4  φ0 dx − C 0 M ||. The desired result

3.2. Error estimate for the first-order scheme In this subsection, we derive the error estimate in time for the first-order scheme (3.6). We suppose that the exact solution satisfies the following regularity:

  

φ ∈ W 2,2 0, T ; L 2 () ∩ W 1,2 0, T ; H 2 () ∩ L 2 0, T ; H 3 () .

(3.10)

We need the following conclusion [10] for our analysis.





Lemma 3.3. Assume that  is open, bounded, and ∂  is smooth. Take s to be a nonnegative integer. Suppose u ∈ L 2 0, T ; H s+2 ()     with ut ∈ L 2 0, T ; H s () . Then u ∈ C [0, T ]; H s+1 () and



max u (t ) H s+1 () ≤ C 3 u  L 2 0, T : H s+2 () + ut  L 2 (0, T ; H s ()) ,

0≤t ≤ T

where C 3 is a constant depending only on T , s and .









With the regularity assumption (3.10), it is clear that φ ∈ L 2 0, T ; H 3 () , φt ∈ L 2 0, T ; H 1 () . With the help of the above Lemma, we have

max φ(t ) H 2 ≤ C 4 .

(3.11)

0≤t ≤ T

Before further investigation, we introduce a discrete Gronwall inequality [17]. Lemma 3.4. Fix T ≥ 0. Let N be a positive integer and τ ≤ such that τ

 N −1

an + τ

m=0 cm

n 

T . N

N N N Suppose {am }m =0 , {bm }m=0 and {cm }m=0 are non-negative sequences

≤ C 5 , where C 5 is independent of τ and N. Further suppose that

bm ≤ C 6 + τ

m =0

n −1 

am cm ,

∀1 ≤ n ≤ N ,

m =0

where C 6 > 0 is a constant independent of τ and N. Then, for all τ > 0,

an + τ

n  m =0

 bm ≤ C 6 exp

τ

n −1 

 cm

≤ C 6 exp (C 5 ) ,

∀1 ≤ n ≤ N .

m =0

We introduce a useful inequality which will be used in the later analysis.

140

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

Lemma 3.5. Denoting e φ,n = φ(tn ) − φn , n ≥ 1, we have

∇ e φ,n 2 ≤ for 4γ − α > 0,

γ 4αγ 2 ∇ e φ,n 2 , e φ,n 2 + α (4γ − α ) 4γ − α

α > 0, γ > 0.

Proof. Since φ(t ) and φn satisfy the boundary conditions (2.6d), we have ∂n e φ,n |∂  = 0. Using Cauchy’s inequality, for α > 0, γ > 0, we have

∇ e φ,n 2 = − (e φ,n , e φ,n ) ≤

1 4α

e φ,n 2 = − (∇ e φ,n , ∇ e φ,n ) ≤

e φ,n 2 + α e φ,n 2 , 1 4γ

∇ e φ,n 2 + γ ∇ e φ,n 2 .

Combining the two inequalities above, we get

∇ e φ,n 2 ≤

1 4α

e φ,n 2 +

then the desired result follows.

α 4γ

∇ e φ,n 2 + αγ ∇ e φ,n 2 ,

2

Now let us give the error estimate of the first-order scheme (3.6). Theorem 3.3. Suppose that the exact solution satisfies (3.10). Then





φ(tn ) − φn  + 4τ 1 − ε −

 n

1 10(1 − ε )2

1

2

2

∇ (φ(tk ) − φk ) 

≤ C7τ , 1 ≤ n ≤ N ,

k =1

where C 7 is a constant independent of τ . Proof. We denote e φ,n = φ(tn ) − φn , e μ,n = μ(tn ) − μn , e w ,n = w (tn ) − w n . At t = tn+1 , the exact solution satisfies

φ(tn+1 ) − φ(tn )

τ

= μ(tn+1 ) − R 1,n+1 ,

(3.12)

μ(tn+1 ) =(1 − ε + β)φ(tn+1 ) +  w (tn+1 ) − (βφ(tn ) − 2φ(tn ) − f (φ(tn ))) − R 2,n+1 + R 3,n+1 + R 4,n+1 , w (tn+1 ) = φ(tn+1 ),

(3.13) (3.14)

where

R 1,n+1 = φt (tn+1 ) −

φ(tn+1 ) − φ(tn )

τ

(3.15)

,

R 2,n+1 = β (φ(tn+1 ) − φ(tn )) ,

(3.16)

R 3,n+1 = 2 (φ(tn+1 ) − φ(tn )) ,

(3.17)

R 4,n+1 = f (φ(tn+1 )) − f (φ(tn )) .

(3.18)

By the Taylor expansion and the Cauchy–Schwarz inequality, we have the following estimates 2

 R 1,n+1  ≤

τ 3

t n +1  φtt (s)2 ds, tn

t n +1   R 2,n+1  ≤ β τ φt (s)2 ds, 2

2

tn

t n +1  2  R 3,n+1  ≤ 4τ φt (s)2 ds, tn

⎧  ⎨9τ ttn+1 φ 2 (s)φt (s)2 ds, φ ∈ [− M , M ], n  R 4,n+1 2 ≤ ⎩9τ M 4  tn+1 φ (s)2 ds, |φ| > M . t tn

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

141

Using the embedding of H 2 () to L ∞ () and the inequality (3.11), we have t n +1 t n +1   2 2 9τ φ (s)φt (s) ds ≤ 9τ φ(s)4L ∞ φt (s)2 ds tn

tn t n +1 t n +1   4 2 4 ≤ 9τ φ(s) H 2 φt (s) ds ≤ 9τ C 4 φt (s)2 ds. tn

tn

So, t n +1   R 4,n+1  ≤ C 8 τ φt (s)2 ds, 2

tn

where C 8 = max{9C 44 , 9M 4 }. Then combining (3.6a)–(3.6c) and (3.12)–(3.14) reaches

e φ,n+1 − e φ,n

τ

= e μ,n+1 − R 1,n+1 ,

(3.19)

e μ,n+1 =(1 − ε )e φ,n+1 + β(e φ,n+1 − e φ,n ) + e w ,n+1 + 2e φ,n

(3.20)

+ f (φ(tn )) − f (φn ) − R 2,n+1 + R 3,n+1 + R 4,n+1 , e w ,n+1 = e φ,n+1 .

(3.21)

Taking the inner product of (3.19) with e φ,n+1 , the inner product of (3.20) with e φ,n+1 , and using the equality (3.21), we have



e φ,n+1 − e φ,n

τ

, e φ,n+1 = (e μ,n+1 , e φ,n+1 ) − ( R 1,n+1 , e φ,n+1 ),

(e μ,n+1 , e φ,n+1 ) =(1 − ε )(e φ,n+1 , e φ,n+1 ) + β(e φ,n+1 − e φ,n , e φ,n+1 )   − (∇ e φ,n+1 , ∇ e φ,n+1 ) + 2(e φ,n , e φ,n+1 ) + f (φ(tn )) − f (φn ), e φ,n+1 − ( R 2,n+1 , e φ,n+1 ) + ( R 3,n+1 , e φ,n+1 ) + ( R 4,n+1 , e φ,n+1 ). After integration by parts and adding the two relations above, we arrive at

1 2τ

 e φ,n+1 2 − e φ,n 2 + e φ,n+1 − e φ,n 2 + (1 − ε )∇ e φ,n+1 2

β

(3.22)



∇ e φ,n+1 2 − ∇ e φ,n 2 + ∇ e φ,n+1 − ∇ e φ,n 2 + ∇ e φ,n+1 2   =2(e φ,n , e φ,n+1 ) + f (φ(tn )) − f (φn ), e φ,n+1 − ( R 2,n+1 , e φ,n+1 ) +

2

+ ( R 3,n+1 , e φ,n+1 ) + ( R 4,n+1 , e φ,n+1 ) − ( R 1,n+1 , e φ,n+1 )  J 1 + J 2 + J 3 + J 4 + J 5 + J 6. Now, we investigate the six terms J i , 1 ≤ i ≤ 6, one by one. For J 1 , using Cauchy’s inequality, we have

J 1 = 2(e φ,n , e φ,n+1 ) = −2(∇ e φ,n , ∇ e φ,n+1 ) ≤ 10∇ e φ,n 2 + Applying Lemma 3.5 with

J1 ≤

121 4

α=

1 , 11

1 10

∇ e φ,n+1 2 .

γ = 14 , we get

1

1

4

10

e φ,n 2 + ∇ e φ,n 2 +

∇ e φ,n+1 2 .

1 For J 2 , defining φ 0 = | | (φ0 , 1) and noticing (φ(tn ), 1) = (φ0 , 1) = (φn , 1), we have (φ(tn ) − φ 0 , 1) = (φn − φ 0 , 1) = 0. Using Hölder’s inequality and the embedding of H1 () to L 6 () yields

142

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

(φ 3 (tn ) − φn3 , e φ,n+1 )

 ≤  φ(tn )2 + φ(tn )φn + φn2 (φ(tn ) − φn )  · e φ,n+1 

 ≤ 2 φ(tn )2 + φn2  L 3 e φ,n L 6 e φ,n+1 

 ≤ 2 φ(tn )2L 6 + φn 2L 6 e φ,n L 6 e φ,n+1 

 1 ≤ 4 φ(tn ) − φ 0 2L 6 + φn − φ 0 2L 6 + 2|| 3 |φ 0 |2 e φ,n L 6 e φ,n+1 

 ≤ C 9 ∇φ(tn )2 + ∇φn 2 + C 10 ∇ e φ,n  · e φ,n+1  ≤ 2C 11 ∇ e φ,n  · e φ,n+1  2 ≤ 10(1 − ε )C 11 ∇ e φ,n 2 + 2 ≤ 10(1 − ε )C 11 ∇ e φ,n 2 + 2 ∇ e φ,n 2 + ≤ 10(1 − ε )C 11

Applying Lemma 3.5 with

α=

1

e φ,n+1 2

10(1 − ε )  1 10(1 − ε ) 1 40(1 − ε )2

1 , 2 10(1−ε )C 11 +1

 (φ 3 (tn )−φn3 , e φ,n+1 ) ≤

1 4(1 − ε )



∇ e φ,n+1 2 + (1 − ε )∇ e φ,n+1 2

∇ e φ,n+1 2 +

1 10

∇ e φ,n+1 2 .

γ = 14 , we have

2 1 + 10(1 − ε )C 11

2

4

1

1

4

40(1 − ε )2

e φ,n 2 + ∇ e φ,n 2 +

∇ e φ,n+1 2 +

On the other hand, we have



3M 2 e φ,n , e φ,n+1 ≤

≤

45(1 − ε ) M 4 2 45(1 − ε ) M 4 2

e φ,n 2 + e φ,n 2 +

1 10(1 − ε )

e φ,n+1 2

1 40(1 − ε )2

∇ e φ,n+1 2 +

1 10

∇ e φ,n+1 2 .

So,

J2 ≤  C 11 e φ,n 2 + 

1 4

∇ e φ,n 2 +

2 1+10(1−ε )C 11 4

where  C 11 = max{ For J 3 , there holds

2

1 40(1 − ε )2

∇ e φ,n+1 2 +

1 10

∇ e φ,n+1 2 ,

, 45(1−2ε)M }. 4

5(1 − ε )

1

 R 2,n+1 2 + e φ,n+1 2 10(1 − ε )  5(1 − ε ) 1 2 2 2 ≤  R 2,n+1  + ∇ e φ,n+1  + (1 − ε )∇ e φ,n+1  2 10(1 − ε ) 4(1 − ε ) 5(1 − ε ) 1 1  R 2,n+1 2 + ∇ e φ,n+1 2 + ∇ e φ,n+1 2 . ≤ 2 2 10 40(1 − ε )

J 3 = −( R 2,n+1 , e φ,n+1 ) ≤

2 1

By the similar analysis of J 3 , we get

J4 ≤ J5 ≤

5(1 − ε ) 2 5(1 − ε ) 2

 R 3,n+1 2 +  R 4,n+1 2 +

1 40(1 − ε

)2

1 40(1 − ε )2

∇ e φ,n+1 2 + ∇ e φ,n+1 2 +

1 10 1 10

For J 6 , we have

J 6 ≤  R 1,n+1 2 +

1 4

e φ,n+1 2 .

Combining the above inequalities with (3.22), we arrive at

∇ e φ,n+1 2 , ∇ e φ,n+1 2 .

1 10

∇ e φ,n+1 2 .

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164



1 2τ

+

 e φ,n+1 2 − e φ,n 2 + e φ,n+1 − e φ,n 2 + 1 − ε −

β

143



1 10(1 − ε )2

∇ e φ,n+1 2

(3.23)

 1  ∇ e φ,n+1 2 − ∇ e φ,n 2 + ∇ e φ,n+1 − ∇ e φ,n 2 + ∇ e φ,n+1 2 − ∇ e φ,n 2

2

2

t n +1   1 121 τ 2 2 e φ,n  + C 11 + φtt 2 dt ≤ e φ,n+1  + 

4

4

3

tn

+

5τ β 2 (1 − ε )

t n +1 

t n +1 t n +1   5C 8 τ (1 − ε ) 2 φt  dt + 10τ (1 − ε ) φt  dt + φt 2 dt . 2

2

2

tn

tn

tn

Summing (3.23) from n = 0 to n = m and noticing e φ,0 = 0, we arrive at

 2

2

e φ,m+1  + β τ ∇ e φ,m+1  + τ ∇ e φ,m+1  + 2τ 1 − ε − ≤

m τ

2

1

2

 m

10(1 − ε )2

 m 121  e φ,n+1 2 + 2τ  C 11 + e φ,n 2 + R , 4

n =0

∇ e φ,n+1 2

n =0

n =0

where

R=

2τ 2

T

T 2

2

2

φtt  dt + 5τ β (1 − ε )

3 0

T 2

2

0

T + 5C 8 τ 2 (1 − ε )

φt 2 dt ≤

φt 2 dt

φt  dt + 20τ (1 − ε ) 0

C 12 + β 2 (1 − ε )C 13 + (1 − ε )C 14 2

τ 2.

0

If

τ satisfies τ < 1, the inequality (3.25) implies that 1 2

 e φ,m+1 2 + 2τ 1 − ε − 

≤ τ 2 C 11 + 61

m 

m +1

1 10(1 − ε )2

e φ,n 2 +

∇ e φ,n 2

n =1

C 12 + β (1 − ε )C 13 + (1 − ε )C 14 2

2

n =1

τ 2.

Then by using the discrete Gronwall inequality, we get

 e φ,m+1 2 + 4τ 1 − ε −

1

m +1

10(1 − ε )2

∇ e φ,n 2

n =1



  ≤ C 12 + β (1 − ε )C 13 + (1 − ε )C 14 τ 2 exp 4 C 11 T + 122T , 2

which implies the desired result.

2

Furthermore, we can get the error estimate for ∇ e φ,n , which is useful in the following analysis. Theorem 3.4. Suppose that the exact solution satisfies (3.10). Then

 ∇ e φ,n  + 4τ (

19 20

− ε)

n 

1

2

e φ,n 2

≤ C 15 τ ,

k =1

where C 15 is a constant independent of τ . Proof. The proof is similar to that in Theorem 3.3, we omit it for brevity.

2

(3.24)

144

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

4. The second-order scheme based on the SISDC method The spectral deferred correction (SDC) method was first introduced by Dutt et al. [7] to solve initial value ordinary differential equations. The basic idea of the SDC method is to convert the original ODEs into the corresponding Picard equations and then apply a deferred correction procedure in the integral formulation, aiming to achieve higher-order accuracy in an iterative way. Then, a semi-implicit spectral deferred correction (SISDC) method was introduced in [21]. In this section, we first give a review about the SISDC method, and then propose the second-order scheme based on the SISDC method for the PFC equation. 4.1. The SISDC method Consider the ODE equation:

 φt = F (t , φ), t ∈ (0, T ], φ(0) = φ0 .

(4.1)

Recall that 0 = t 0 < t 1 < · · · < t N = T is a uniform partition of [0, T ], τ = NT . The SISDC method described in this subsection will be used in every interval [tn , tn+1 ], 0 ≤ n ≤ N − 1. We use the first-order scheme (3.6) to obtain an approximate solution φ 0 , let φn0 represent φ 0 (tn ), φn1 represent a corrected approximation at tn . The corresponding Picard integral equation of (4.1) is

t φ(t ) = φ0 +

F (s, φ(s)) ds.

(4.2)

0

Then a residual function is defined by

t 0



F s, φ 0 (s) ds − φ 0 (t ).

r (t , φ ) = φ0 +

(4.3)

0

By denoting the error function by δ(t ) = φ(t ) − φ 0 (t ), we can rewrite (4.2) in the form

t



F s, φ 0 (s) + δ(s) ds − φ 0 (t ).

δ(t ) = φ0 +

(4.4)

0

Combining (4.3) with (4.4), we get the correction equation

t







F s, φ 0 (s) + δ(s) − F s, φ 0 (s) ds + r (t , φ 0 ).

δ(t ) =

(4.5)

0

Then, we use some numerical method to discretize the correction equation (4.5), and get a higher-order approximate solution φ 1 (t ) by φ 1 (t ) = φ 0 (t ) + δ(t ). Suppose F (t , φ) = F I (t , φ) + F E (t , φ). Then we treat F I implicitly and F E explicitly. Applying the semi-implicit Euler method to the correction equation (4.5), we have



δn+1 = δn + τ F I (tn+1 , φn0+1 + δn+1 ) + F E (tn , φn0 + δn ) 

− τ F I (tn+1 , φn0+1 ) + F E (tn , φn0 ) + rn+1 (φ 0 ) − rn (φ 0 ),

(4.6)

where δi = δ(τi ), r i (φ 0 ) = r (τi , φ 0 ). We take the notation

I nn+1 (φ 0 )

t n +1 

F (s, φ 0 (s))ds.

=

(4.7)

tn

Then it follows from (4.3) that

I nn+1 (φ 0 ) = rn+1 (φ 0 ) − rn (φ 0 ) + φn0+1 − φn0 .

(4.8)

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

145

The formulas (4.6), (4.8) together with φn0+1 + δn+1 = φn1+1 yield



  φn1+1 = φn1 + τ F I (tn+1 , φn1+1 ) + F E (tn , φn1 ) − τ F I (tn+1 , φn0+1 ) + F E (tn , φn0 ) + I nn+1 (φ 0 ).

(4.9)

Next, we introduce one way for the computation of I nn+1 (φ 0 ) in (4.7). Denoting the p + 1 Gauss–Lobatto nodes (cf. [23]) on [tn , tn+1 ] by tn = ξ0 < ξ1 < · · · < ξ p = tn+1 . The Lagrange interpolap tion polynomials associated with {ξi }i =0 are given by

L j (s) =

p 

s − ξi

i =0,i = j

ξ j − ξi

, 0 ≤ j ≤ p.

Then,

I nn+1 (φ 0 ) ≈

 tn+1

p 

F (ξ j , φ 0 (ξ j ))ω j ,

j =0

where ω j = t L j (s)ds. n In the rest of this section, we apply the SISDC method to the PFC equation. 4.2. The second-order scheme With the first-order convex splitting numerical method (3.6), we obtain an approximate solution φ 0 . As the first-order method is indeed a semi-implicit scheme, it is natural to apply the same explicit-implicit strategy to the SISDC method. In the PFC equation,



F (t , φ) :=  f (φ) + (1 − ε )φ + 2φ + 2 φ . We take



F I (t , φ) =  (1 − ε + β)φ + 2 φ , F E (t , φ) = − (βφ − 2φ − f (φ)) . Then, given φ01 = φ00 , the SISDC scheme for the PFC equation reads:



  φn1+1 = φn1 + τ F I (tn+1 , φn1+1 ) + F E (tn , φn1 ) − τ F I (tn+1 , φn0+1 ) + F E (tn , φn0 ) + I nn+1 (φ 0 ).

We take the following formula to compute I nn+1 (φ 0 ),

I nn+1 (φ 0 ) =

τ 2



F (tn , φn0 ) + F (tn+1 , φn0+1 ) .

Now we give an SISDC correction scheme for the first-order scheme (3.6):



φn1+1 = φn1 + τ  (1 − ε + β)φn1+1 + 2 φn1+1 − (βφn1 − 2φn1 − f (φn1 )) 

− τ  (1 − ε + β)φn0+1 + 2 φn0+1 − (βφn0 − 2φn0 − f (φn0 ))  τ + F (tn , φn0 ) + F (tn+1 , φn0+1 ) ,

(4.10)

2

with the boundary conditions ∂n φn1+1 |∂  = 0, ∂n (φn1+1 )|∂  = 0, ∂n (2 φn1+1 )|∂  = 0, which will be shown to be a secondorder correction scheme. We only need to solve the resulted linear system above, it’s easy-to-implement and efficient. We present the unconditional weak energy stability of the scheme (4.10). Theorem 4.1. If β ≥ 3M 2 , the SISDC correction scheme (4.10) is unconditionally energy stable, i.e.,

E (φn1 ) ≤ 3E (φ01 ), 0 ≤ n ≤ N , ∀τ > 0. Proof. The scheme (4.10) can be written as

  φn1+1 = φn1 + τ  (1 − ε )φn1+1 + β(φn1+1 − φn1 ) + 2 φn1+1 + 2φn1 + f φn1   

 f φ 0  − f φ 0  φ 0 − φn0+1 1 − ε + 2β 0 n n +1 0 2 n 0 0 − τ −  φn+1 − φn − φn+1 − φn −  . 2

2

2

(4.11)

146

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164





Taking the inner product of (4.12) with (−)−1 φn1+1 − φn1 and by the similar analysis as that of Theorem 3.2, we have

E (φn1+1 ) − E (φn1 ) +

=−

1

1 − ε + 2β − f  (αn )



φn1+1 − φn1 2−1 +

τ

2 2

≤ 

2

2

ηn is between φn0 and φn0+1 .

 φn0+1 − φn0 , φn1+1 − φn1

1 − ε + 2β − f  (ηn ) 2

≤

2



 φ 0 − φn0+1 φn0+1 − φn0 1 n) 0 0 2 n 1 − 2 , φn+1 − φn , φn+1 − φn − 

f  (η

1 − ε + 2β −

where αn is between φn1 and φn1+1 , Thanks to the inequalities



1

φn1+1 − φn1 2 + φn1+1 − φn1 2 + ∇φn1+1 − ∇φn1 2

1 − ε + 2β − f  (ηn ) 2 1−ε+β 2

2

φn0+1 − φn0 2 +

φn0+1 − φn0 2 + 

φn0 − φn0+1 2

1−ε+β 8

1 − ε + 2β − f  (ηn ) 8

φn1+1 − φn1 2

φn1+1 − φn1 2 ,

1

1

2

8

, φn1+1 − φn1 ≤ φn0+1 − φn0 2 + φn1+1 − φn1 2 ,

  1  φn0+1 − φn0 , φn1+1 − φn1 ≤∇φn0+1 − ∇φn0 2 + ∇φn1+1 − ∇φn1 2 , 4

we get

E (φn1+1 ) − E (φn1 ) +

3 (1 − ε + β) 8 3

3

φn1+1 − φn1 2 +

β − f  (αn ) 2

φn1+1 − φn1 2

(4.12)

+ φn1+1 − φn1 2 + ∇φn1+1 − ∇φn1 2 ≤−

8 1

φn1+1

τ

− φn1 2−1

+

4 1−ε+β 2

1

φn0+1 − φn0 2 + φn0+1 − φn0 2 + ∇φn0+1 − ∇φn0 2 . 2

Summing (4.12) from n = 0 to n = m, (m ≤ N − 1), we arrive at 1 1 E (φm +1 ) − E (φ0 ) + m   3

+

n =0

≤−

m 1

τ

n =0

8

m   3 (1 − ε + β)

8

n =0

φn1+1

− φn1 2

φn1+1 − φn1 2−1 +

+

3 4

φn1+1 − φn1 2 +

β − f  (αn ) 2

φn1+1 − φn1 2

(4.13)

∇φn1+1

− ∇φn1 2

m   1−ε+β

2

n =0

1 φn0+1 − φn0 2 + φn0+1 − φn0 2 + ∇φn0+1 − ∇φn0 2 . 2

On the other hand, summing (3.8) from n = 0 to n = m, we get

−τ

m  n =0

+

0 0 ∇ μn0+1 2 = E (φm +1 ) − E (φ0 ) +

m   1−ε+β

2

n =0

φn0+1

− φn0 2

+

1 2

m  β − f  (ξn ) n =0

φn0+1

2

φn0+1 − φn0 2

(4.14)

− φn0 2

+ ∇φn0+1

− ∇φn0 2

.

With the help of (3.7) and (4.15), we have m   1−ε+β n =0

2

φn0+1

− φn0 2

+

1 2

φn0+1

− φn0 2

+ ∇φn0+1

1 1 Thanks to (4.14) and (4.15), we arrive at E (φm +1 ) ≤ 3E (φ0 ), ∀τ > 0.

− ∇φn0 2

2

≤ 2E (φ00 ).

(4.15)

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

147

4.3. Error estimate for the second-order scheme In this subsection, we derive the error estimate for the correction scheme (4.10). We assume that the exact solution satisfies the following regularity:

 

φ ∈ W 3,2 0, T ; L 2 () ∩ W 2,2 0, T ; H 6 () .

(4.16)    0, T ; H 5 () , φt ∈ L 2 0, T ; H 3 () ; and φt ∈

 2

 With the regularity  assumption (4.16), it is clear that φ ∈ L L 2 0, T ; H 3 () , φtt ∈ L 2 0, T ; H 1 () . With the help of Lemma 3.3, we obtain max φ H 4 ≤ C 16 ;

0≤t ≤ T

max φt  H 2 ≤ C 17 .

(4.17)

0≤t ≤ T

Let S n+1 = e φ,n+1 − e φ,n , 0 ≤ n ≤ N − 1. First, we try to get an estimate for S n+1 , which is useful to the error estimate for the correction scheme (4.10). Lemma 4.1. Suppose that the exact solution satisfies (4.16). Then





 S n +1  + 4τ 1 − ε −

n +1

9 80(1 − ε )2

 12 2

∇ S k 

≤ C 18 τ 2 ,

0 ≤ n ≤ N − 1,

k =2

where C 18 is a constant independent of τ . Proof. The equations (3.19), (3.20) and (3.21) imply

e φ,n+1 − e φ,n

τ

 = (1 − ε )e φ,n+1 + β(e φ,n+1 − e φ,n ) + 2 e φ,n+1 + 2e φ,n

 +  f (φ(tn )) − f (φn0 ) − R 2,n+1 + R 3,n+1 + R 4,n+1 − R 1,n+1 ,

(4.18)

and at t = tn , it reads

e φ,n − e φ,n−1

τ

 = (1 − ε )e φ,n + β(e φ,n − e φ,n−1 ) + 2 e φ,n + 2e φ,n−1

  +  f (φ(tn−1 )) − f φn0−1 − R 2,n + R 3,n + R 4,n − R 1,n .

(4.19)

Subtracting (4.19) from (4.18), we get

S n +1 − S n

τ

 = (1 − ε ) S n+1 + β( S n+1 − S n ) + 2 S n+1 + 2 S n



  +  f (φ(tn )) − f φn0 − f (φ(tn−1 )) + f φn0−1 − R 2 + R 3 + R 4 − R 1 ,

where

R 1 = R 1,n+1 − R 1,n = φt (tn+1 ) − φt (tn ) −

φ(tn+1 ) − 2φ(tn ) + φ(tn−1 )

τ

,

R 2 = R 2,n+1 − R 2,n = β (φ(tn+1 ) − 2φ(tn ) + φ(tn−1 )) , R 3 = R 3,n+1 − R 3,n = 2 (φ(tn+1 ) − 2φ(tn ) + φ(tn−1 )) , R 4 = R 4,n+1 − R 4,n = f (φ(tn+1 )) − 2 f (φ(tn )) + f (φ(tn−1 )) . By the Taylor expansion and the Cauchy–Schwarz inequality, we have the following estimates t n +1   R 1  ≤2τ φttt (s)2 ds, 2

3

t n −1

t n +1   R 2  ≤ 2β τ φtt (s)2 ds, 2

2 3

t n −1

t n +1 

 R 3 2 ≤8τ 3 t n −1

φtt (s)2 ds,

⎧  ⎨6τ 3 ttn+1 2φ(s)φt2 (s) + φ 2 (s)φtt (s)2 ds, φ ∈ [− M , M ], n −1 2 R4 ≤ ⎩18τ 3 M 4  tn+1 φ (s)2 ds, |φ| > M . tt t n −1

(4.20)

148

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

With the help of the embedding of H 2 () to L ∞ () and (4.17), we obtain t n +1 t n +1   2

2 2 2 3 6τ 2φ(s)φt (s) + φ (s)φtt (s) ds ≤ 6τ 2φ(s)φt2 (s) + φ 2 (s)φtt (s) ds 3

t n −1

t n −1

≤ 6τ

3

t n +1 

2

1

2φ(s) L ∞ φt (s)2L ∞ || 2 + φ(s)2L ∞ φtt (s)

ds

t n −1

≤ C 19 τ So,  R 4 2 ≤ C 20 τ

3

t n +1 

2 φ(s) H 2 φt (s)2H 2 + φ(s)2H 2 φtt (s) ds.

t n −1

 3 tn+1

φtt (s)2 ds + C 21 τ 4 .     Next we estimate f (φ(tn )) − f φn0 − f (φ(tn−1 )) + f φn0−1 . When |φ| > M, it’s easy to handle. Here we just concern about the case of φ ∈ [− M , M ] for brevity. Using a3 − b3 = (a − b)3 + 3ab(a − b), we get tn−1

φ 3 (tn ) − (φn0 )3 − φ 3 (tn−1 ) + (φn0−1 )3

3

 3

 = φ(tn ) − φn0 + 3φ(tn )φn0 φ(tn ) − φn0 − φ(tn−1 ) − φn0−1 − 3φ(tn−1 )φn0−1 φ(tn−1 ) − φn0−1 =e 3φ,n − e 3φ,n−1 + 3φ(tn )φn0 e φ,n − 3φ(tn−1 )φn0−1 e φ,n−1    = e φ,n − e φ,n−1 e 2φ,n + e φ,n e φ,n−1 + e 2φ,n−1 − 3φ(tn )e φ,n e φ,n + 3φ(tn )e φ,n−1 e φ,n + 3φ(tn ) (φ(tn ) − φ(tn−1 )) e φ,n + 3φn0−1 (φ(tn ) − φ(tn−1 )) e φ,n + 3φ(tn−1 )φn0−1 (e φ,n − e φ,n−1 ) I1 + I2 + I3 + I4 + I5 + I6. For I 1 , there holds

    I 1  = e φ,n − e φ,n−1 e 2φ,n + e φ,n e φ,n−1 + e 2φ,n−1  ≤e 2φ,n + e φ,n e φ,n−1 + e 2φ,n−1 L 3 e φ,n − e φ,n−1 L 6

 ≤2 e φ,n 2L 6 + e φ,n−1 2L 6 e φ,n − e φ,n−1 L 6

 ≤C 22 ∇ e φ,n 2 + ∇ e φ,n−1 2 ∇ e φ,n − ∇ e φ,n−1  ≤ C 23 τ 3 .

Furthermore, we have

 I 2  =3φ(tn )e φ,n e φ,n  ≤ 3φ(tn ) L ∞ e φ,n L 3 e φ,n L 6 ≤ 3C 24 φ(tn ) H 2 ∇ e φ,n 2 ≤ C 25 τ 2 ,  I 3  =3φ(tn )e φ,n−1 e φ,n  ≤ 3φ(tn ) L ∞ e φ,n−1 L 3 e φ,n L 6 ≤ C 26 τ 2 ,  I 4  = 3φ(tn ) (φ(tn ) − φ(tn−1 )) e φ,n  ≤ 3φ(tn ) L ∞ φ(tn ) − φ(tn−1 ) L 3 e φ,n L 6 ⎛ ⎞1 2 tn 3 ⎜ ⎟ 2 ≤ C 27 ∇φ(tn ) − ∇φ(tn−1 ) · ∇ e φ,n  ≤ C 28 τ 2 ⎝ |∇φt (s)| ds⎠ , t n −1

I5 = ≤

3φn0−1 (φ(tn ) − φ(tn−1 )) e φ,n 

≤

3φn0−1 L ∞ φ(tn ) − φ(tn−1 ) L 3 e φ,n L 6 ⎛

3C 29 φn0−1  H 2 φ(tn ) − φ(tn−1 ) L 3 e φ,n  L 6

⎜ ≤ C 30 τ ⎝ 3 2

tn

⎞1

2

⎟ |∇φt (s)| ds⎠ , 2

t n −1

I6 =

3φ(tn−1 )φn0−1 (e φ,n

− e φ,n−1 ) ≤

3φ(tn−1 ) L ∞ φn0−1  L ∞ e φ,n

− e φ,n−1  ≤ C 31  S n .

Taking the inner product of (4.20) with S n+1 , we get



S n +1 − S n

τ



 , S n+1 = (1 − ε ) S n+1 + β( S n+1 − S n ) + 2 S n+1 + 2 S n ,  S n+1 + ( I 1 + I 2 + I 3 + I 4 + I 5 + I 6 − R 2 + R 3 + R 4 ,  S n +1 ) − ( R 1 , S n +1 ) .

(4.21)

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

149

Furthermore, using Cauchy’s inequality, we have

1

( R 1 , S n+1 ) ≤  R 1 2 +  S n+1 2 , 4

( R i ,  S n+1 ) ≤ 5(1 − ε ) R i 2 + ( I i ,  S n+1 ) ≤ 5(1 − ε ) I i 2 +

1 80(1 − ε )2 1

80(1 − ε )2

2 (1 − ε ) S n 2 + ( I 6 ,  S n+1 ) ≤ 5C 31

121

2 ( S n ,  S n+1 ) ≤

2

1

∇ S n+1 2 +

∇ S n+1 2 +

1 80(1 − ε

)2

1 20

1

2

20

∇  S n+1 2 ,

∇  S n+1 2 ,

∇ S n+1 2 +

1

 S n 2 + ∇  S n 2 +

20

1 20

i = 2, 3, 4, i = 1, 2, 3, 4, 5,

∇  S n+1 2 ,

∇  S n+1 2 .

Combining the above inequalities with (4.21), we arrive at

   S n+1 2 −  S n 2 +  S n+1 − S n 2 + 1 − ε −

1 2τ

80(1 − ε )2

∇ S n+1 2

(4.22)

 1  ∇ S n+1 2 − ∇ S n 2 + ∇ S n+1 − ∇ S n 2 + ∇  S n+1 2 − ∇  S n 2 2 2  4   5   1 121 2 2 2 2 2 =  S n +1  + + 5C 31 (1 − ε )  S n  + 5(1 − ε ) Ri  +  I i  +  R 1 2 . +

β

9

4

2

i =2

i =1

Summing (4.22) from n = 1 to n = m, we get

 2

2

 S m +1  −  S 1  + 2τ 1 − ε −

m +1

9 80(1 − ε )2

∇ S n 2

(4.23)

n =2



 + β τ ∇ S m+1 2 − ∇ S 1 2 + τ ∇  S m+1 2 − ∇  S 1 2  m τ 1  2 (1 − ε ) +  S n 2 + Re , ≤  S m+1 2 + τ 121 + 10C 31 2

2

n =1

where

T 2

Re ≤20β (1 − ε )τ

4

T 2

φtt (s) ds + 80(1 − ε )τ 0

4

φtt (s)2 ds + 10T (1 − ε )C 23 τ 6 0

T + 10(1 − ε )C 20 τ 4

φtt (s)2 ds + 10(1 − ε )C 21 τ 5 + 10T (1 − ε ) (C 25 + C 26 ) τ 4 0

T + 10(1 − ε ) (C 28 + C 30 ) τ 4

T ∇φt (s)2 ds + 4τ 4

0

φttt (s)2 ds ≤ C 32 τ 4 . 0

Next, we try to bound the term  S 1 2 + β τ ∇ S 1 2 + τ ∇  S 1 2 . Since S 1 = e φ,1 − e φ,0 = e φ,1 , we have the following equation from (4.18):

S1

τ

 =  (1 − ε ) S 1 + β S 1 + 2 S 1 − R 2,1 + R 3,1 + R 4,1 − R 1,1 ,

(4.24)

where R 1,1 , R 2,1 , R 3,1 , R 4,1 is defined in (3.15)–(3.18).       regularity assumption (4.16), we have φtt ∈ L 2 0, T ; H 2 () , φttt ∈ L 2 0, T ; L 2 () ; φt ∈ L 2 0, T ; H 5 () , φtt ∈  With the 2 3 L 0, T ; H () , then

max φtt  H 1 () ≤ C 33 ;

0≤t ≤ T

max φt  H 4 ≤ C 34 .

0≤t ≤ T

150

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

It follows that

2

 R 1,1  ≤

t1

τ

2

φtt  dt ≤

3

t1

τ

φtt 2H 1 dt ≤ C 35 τ 2 ,

3

t0

t0

t1 2

2

t1 2

 R 2,1  ≤β τ

2

φt 2H 2 dt ≤ C 36 τ 2 ,

φt  dt ≤ β τ t0

t0

t1  R 3,1 2 ≤4τ

t1 2 φt 2 dt ≤ 4τ

t0

φt 2H 4 dt ≤ C 37 τ 2 , t0



 φ 3 (t 1 ) − φ 3 (t 0 ) 2 ≤τ

t1 6φt (∇φ)2 + 12φ (∇φ) (∇φt ) + 6φφt φ + 3φ 2 φt 2 dt t0

t1 4φt 2L ∞ ∇φ4L ∞ || + 16φ2L ∞ ∇φ2L ∞ ∇φt 2L ∞

≤36τ t0

+ 4φ2L ∞ φt 2L ∞ φ2L ∞ + φ4L ∞ φt 2L ∞ dt t1 φt 2H 2 ∇φ4H 2 + φ2H 2 ∇φ2H 2 ∇φt 2H 2

≤C 38 τ t0

+ φ2H 2 φt 2H 2 φ2H 2 + φ4H 2 φt 2H 2 dt , t1 3M 2  (φ(t 1 ) − φ(t 0 )) 2 ≤ 9M 4 τ

t1 φt 2 dt ≤ 9M 4 τ

t0 2

φt 2H 2 dt , t0

2

 R 4,1  ≤ C 39 τ . Taking the inner product of (4.24) with S 1 , we get

 S 1 2 + τ (1 − ε )∇ S 1 2 + β τ ∇ S 1 2 + τ ∇  S 1 2         = − τ  R 2,1 , S 1 + τ  R 3,1 , S 1 + τ  R 4,1 , S 1 − τ R 1,1 , S 1 1

1

2

2

≤2τ 2  R 2,1 2 + 2τ 2  R 3,1 2 + 2τ 2  R 4,1 2 + 2τ 2  R 1,1 2 +  S 1 2 ≤ C 40 τ 4 +  S 1 2 , which implies that

 S 1 2 + β τ ∇ S 1 2 + τ ∇  S 1 2 ≤ C 41 τ 4 . If

(4.25)

τ satisfies τ < 1, (4.23) and (4.25) yield that 

9

2

 S m +1  + 4τ 1 − ε −

80(1 − ε )2

m +1

m 

2 ∇ S n 2 ≤ τ 243 + 20C 31 (1 − ε )  S n 2 + C 42 τ 4 . n =1

n =2

Then we apply the discrete Gronwall inequality and get

  S m+1 2 + 4τ 1 − ε −

9 80(1 − ε )2

which implies the desired result.

m +1



2 ∇ S n 2 ≤ C 42 τ 4 exp T 243 + 20C 31 (1 − ε ) ,

n =2

2

The following theorem is the main result of this subsection.

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

151

Theorem 4.2. Suppose that the exact solution satisfies (4.16). Then

 φ(tn ) − φn1  +

 4τ

1−ε−

 n

1 16(1 − ε )2

  ∇ φ(tk ) − φk1 2

1

2

≤ C 43 τ 2 ,

1 ≤ n ≤ N,

k =1

where C 43 is a constant independent of τ . Proof. We denote e 1φ,n = φ(tn ) − φn1 , e 1μ,n = μ(tn ) − μn1 , e 1w ,n = w (tn ) − w n1 . At t = tn+1 , the exact solution satisfies

τ

φ(tn+1 ) = φ(tn ) + ( F (tn+1 , φ) + F (tn , φ)) + R T 2 

=τ  (1 − ε )φ(tn+1 ) + β (φ(tn+1 ) − φ(tn )) + 2 φ(tn+1 ) + 2φ(tn ) + f (φ(tn ))   φ(tn ) − φ(tn+1 ) φ(tn ) − φ(tn+1 ) + τ  (1 − ε ) − β (φ(tn+1 ) − φ(tn )) + 2 2 2  φ(tn+1 ) − φ(tn ) f (φ(tn+1 )) − f (φ(tn )) + φ(tn ) + R T , + τ  2 + 2

(4.26)

2

3

where R T =C 44 τ F tt (ηn ), ηn ∈ (tn , tn+1 ). The correction scheme (4.10) implies



φn1+1

=

φn1

+ τ  (1 − ε )

φn0 − φn0+1 2

−β

φn0+1

− φn0



+

2

φn0 − φn0+1

 (4.27)

2

  + τ  (1 − ε )φn1+1 + β(φn1+1 − φn1 ) + 2 φn1+1 + 2φn1 + f φn1      φn0+1 − φn0 f φn0+1 − f φn0 + τ  2 + . 2

2

Subtracting (4.27) from (4.26) admits



e 1φ,n+1 = e 1φ,n + τ  (1 − ε )

e φ,n − e φ,n+1 2

  e φ,n − e φ,n+1 − β e φ,n+1 − e φ,n + 2



2

(4.28)

  + τ  (1 − ε )e 1φ,n+1 + β(e 1φ,n+1 − e 1φ,n ) + 2 e 1φ,n+1 + 2e 1φ,n + f (φ(tn )) − f φn1      f (φ(tn+1 )) − f φn0+1 − f (φ(tn )) + f φn0 e φ,n+1 − e φ,n + τ  2 + + RT . 2

2

Taking the inner product of (4.28) with e 1φ,n+1 , we get

1

 e 1φ,n+1 2 − e 1φ,n 2 + e 1φ,n+1 − e 1φ,n 2 + τ (1 − ε )∇ e 1φ,n+1 2 2  τβ + τ ∇ e 1φ,n+1 2 + ∇ e 1φ,n+1 2 − ∇ e 1φ,n 2 + ∇ e 1φ,n+1 − ∇ e 1φ,n 2 2 



   S n +1 1 1 , e 1φ,n+1 =2τ e φ,n , e φ,n+1 + τ f (φ(tn )) − f φn1 , e 1φ,n+1 − τ (1 − ε ) 2  

 S S n +1 n +1 − τ β S n+1 , e 1φ,n+1 + τ ∇  , ∇ e 1φ,n+1 + 2τ  , e 1φ,n+1 2 2       

f (φ(tn+1 )) − f φn0+1 − f (φ(tn )) + f φn0 +τ , e 1φ,n+1 + R T , e 1φ,n+1 2

K1 + K2 + K3 + K4 + K5 + K6 + K7 + K8. For K 1 , using Cauchy’s inequality, we have



1 , 21

γ = 14 , we get



τ

K 1 =2τ e 1φ,n , e 1φ,n+1 = −2τ ∇ e 1φ,n , ∇ e 1φ,n+1 ≤ 20τ ∇ e 1φ,n 2 + ∇ e 1φ,n+1 2 . 20 Applying Lemma 3.5 with

α=

(4.29)

152

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

K1 ≤

441τ 4

τ

e 1φ,n 2 +

4

∇ e 1φ,n 2 +

τ 20

∇ e 1φ,n+1 2 .

For K 2 , with similar analysis of J 2 in Theorem 3.3, we have



τ φ 3 (tn ) − (φn1 )3 , e1φ,n+1 ≤ 2C 45 τ ∇ e1φ,n  · e1φ,n+1  τ 2 ≤20(1 − ε )C 45 τ ∇ e1φ,n 2 + e 1φ,n+1 2 20(1 − ε ) 

≤

2 τ 1 + 20(1 − ε)C 45

2

τ

e 1φ,n 2 +

4

4



∇ e 1φ,n 2 +

τ 80(1 − ε )2

τ 3M 2 φ(t 1 ) − 3M 2 φn1 , eφ,n+1 = 3τ M 2 e1φ,n , e1φ,n+1



≤ 45τ (1 − ε ) M 4 e 1φ,n 2 + ≤ 45τ (1 − ε ) M 4 e 1φ,n 2 + K2 ≤  C 45 τ e 1φ,n 2 + 

τ 4

∇ e 1φ,n 2 +

2 1+20(1−ε )C 45

where  C 45 = max{ 4 By Lemma 4.1, we have

 S n+1  ≤ C 18 τ 2 ,

2

τ 80(1 − ε )2

∇ e 1φ,n+1 2 +

τ 20(1 − ε )

∇ e 1φ,n+1 2 +

τ 20

∇ e 1φ,n+1 2 ,

e φ,n+1 2

τ 80(1 − ε

τ 20

)2

∇ e 1φ,n+1 2 +

∇ e 1φ,n+1 2 ,

, 45(1 − ε ) M 4 }.

0 ≤ n ≤ N − 1.

By the similar analysis as that of Lemma 4.1, we can derive

∇ S n+1  ≤ C 46 τ 2 ,

 S n+1  ≤ C 47 τ 2 ,

∇  S n+1  ≤ C 48 τ 2 .

For K 3 , K 4 , K 5 , K 6 , K 8 , using Cauchy’s inequality, we get

 K 3 = − τ (1 − ε )

≤ ≤

5τ (1 − ε )3 4 5τ (1 − ε )3 4

S n +1 2

5τ (1 − ε )3 , e 1φ,n+1 ≤  S n+1 2 + 4

 S n+1 2 + C 18 τ 4 +

τ 80(1 − ε

τ 80(1 − ε



)2

)2

∇ e 1φ,n+1 2 +

∇ e 1φ,n+1 2 +

τ 20

τ

τ

≤5τ (1 − ε )β 2 C 18 τ 4 +  ∇

K5 = τ

S n +1 2

80(1 − ε )2

τ 80(1 − ε )2

τ





S n +1 2

4

≤5(1 − ε )C 47 τ 5 +

τ 20

∇ e 1φ,n+1 2

∇ e 1φ,n+1 2 ,

20

5

τ

4

20

∇ e 1φ,n+1 2 ,

, e 1φ,n+1 ≤ 5τ (1 − ε ) S n+1 2 +

≤5τ (1 − ε ) S n+1 2 +

K8 =

τ 20

e 1φ,n+1 2

5τ τ , ∇ e 1φ,n+1 ≤ ∇  S n+1 2 + ∇ e 1φ,n+1 2 ≤ C 48 τ 5 +

K 6 =2τ

e 1φ,n+1 2

∇ e 1φ,n+1 2

20(1 − ε )

∇ e 1φ,n+1 2 +

∇ e 1φ,n+1 2 +

20(1 − ε )

∇ e 1φ,n+1 2 ,

20

K 4 = − τ β S n+1 , e 1φ,n+1 ≤ 5τ (1 − ε )β 2  S n+1 2 +

≤5τ (1 − ε )β 2  S n+1 2 +

τ



R T , e 1φ,n+1 ≤

τ 80(1 − ε )2

τ 80(1 − ε )2 1

τ

 R T 2 +

∇ e 1φ,n+1 2 +

∇ e 1φ,n+1 2 +

τ 4

τ 20

τ 20(1 − ε )

τ 20

e 1φ,n+1 2

∇ e 1φ,n+1 2

∇ e 1φ,n+1 2 ,

e 1φ,n+1 2 ≤ C 49 τ 5  F tt (ηn )2 +

τ 4

e 1φ,n+1 2 .

τ 20

∇ e 1φ,n+1 2 ,

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

153

For K 7 , by the similar analysis as that in Lemma 4.1, we have



τ

φ 3 (tn+1 ) − (φn0+1 )3 − φ 3 (tn ) + (φn0 )3

≤ ≤

2

5τ (1 − ε ) 4

, e 1φ,n+1

φ 3 (tn+1 ) − (φn0+1 )3 − φ 3 (tn ) + (φn0 )3 2 + ⎛

⎞ t n +1  ⎝C 50 τ 4 + C 51 τ 3 ∇φt (s)2 ds⎠ +

5τ (1 − ε ) 4

tn

3τ M 2

τ 20(1 − ε )

τ 80(1 − ε

)2

e 1φ,n+1 2

∇ e 1φ,n+1 2 +

τ

5

τ +

4

1

80(1 − ε )2

τ

2

∇ e φ,n+1  +

t n +1  ∇φt (s)2 ds +

4

tn

20

∇ e 1φ,n+1 2 ,

∇ e 1φ,n+1 2 ,

τ

∇ e 1φ,n+1 2 +

80(1 − ε )2

τ 20

∇ e 1φ,n+1 2 ,

45(1−ε ) M 4 C 2

5(1−ε )C

50 18 where  C 50 = max{ , }. 4 4 Combining the above inequalities with (4.29), we arrive at

  e 1φ,n+1 2 − e 1φ,n 2 + e 1φ,n+1 − e 1φ,n 2 + τ 1 − ε −

1

+

20

2

2 45(1 − ε ) M 4 C 18

5(1 − ε )C 51 τ 4 K7 ≤  C 50 τ 5 +

2

τ

 3τ M 2  S n+1 , e 1φ,n+1 φ(tn+1 ) − φn0+1 − φ(tn ) + φn0 , e 1φ,n+1 =

2

≤



τβ



1 16(1 − ε )2

∇ e 1φ,n+1 2

(4.30)

 τ  ∇ e 1φ,n+1 2 − ∇ e 1φ,n 2 + ∇ e 1φ,n+1 − ∇ e 1φ,n 2 + ∇ e 1φ,n+1 2 − ∇ e 1φ,n 2

2

2

t n +1   τ 441 5C 51 (1 − ε )τ 4 e 1φ,n 2 + 5C 18 (1 − ε )β 2 τ 5 + C 45 + ∇φt 2 dt ≤ e 1φ,n+1 2 + τ 

4

4

4

tn

+ 5(1 − ε )C 47 τ 5 +

5(1 − ε ) C 18 + 5C 48 3

4

τ 5 + C 49 τ 5  F tt (ηn )2 +  C 50 τ 5 .

Summing (4.30) from n = 0 to n = m, m ≤ N − 1, noticing that e 1φ,0 = 0 and mτ ≤ T , we get

1 2

 e 1φ,m+1 2 + τ 1 − ε − ≤

τ 4

m +1

1 16(1 − ε )2

∇ e 1φ,n 2 +

n =1

βτ 2

∇ e 1φ,m+1 2 +

τ 2

∇ e 1φ,m+1 2

 m 221  1 2  e φ,m+1  + τ C 45 + e φ,n  + Rc , 1

2

2

n =1

where

Rc ≤

5T (1 − ε )3 C 18 + 5T C 48 4

+ C 49 τ

4

m 

τ 4 + 5T C 18 (1 − ε)β 2 τ 4 + 5T (1 − ε)C 47 τ 4

5C 51 (1 − ε )τ 4 τ  F tt (ηn ) + T  C 50 τ 4 +

1

m +1

16(1 − ε )2

≤ τ 4 C 45 + 442

4

0

τ < 1, (4.32) can be rewritten as  e 1φ,m+1 2 + 4τ 1 − ε − 

1

∇φt 2 dt ≤ C 52 τ 4 .

4

n =0

If

T

2

m 

n =1

∇ e 1φ,n 2 + 2β τ ∇ e 1φ,m+1 2 + 2τ ∇ e 1φ,m+1 2

e 1φ,n 2 + 4Rc .

n =1

Then the application of the discrete Gronwall inequality leads to

(4.31)

154

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

 1

e φ,m+1  + 4τ 1 − ε − which implies the desired result.

m +1

1

2

16(1 − ε )2

  ∇ e 1φ,n 2 ≤ C 52 τ 4 exp 4 C 45 T + 442T ,

n =1

2

5. The fully discrete schemes We use the finite element method for the spatial discretization. Suppose h is a quasi-uniform “triangulation” of the domain  with mesh size h, and  = ∪ K ∈h K . For a positive integer r, let P r ( K ) denote the space of polynomials of degree less than or equal to r on K . We introduce the finite element space





M h = ψh ∈ C 0 (); ψh | K ∈ P r ( K ), r ≥ 1 . We take the notation (φh0 , μh0 , w h0 ) to represent numerical approximate solution solved by the first-order convex splitting

scheme, and φh0,n = φh0 (tn ), μh0,n = μh0 (tn ), w h0,n = w n0 (tn ). The notation (φh1 , μh1 , w h1 ) denotes the numerical approximate solution gotten from the second-order SISDC method, with φh1,n = φh1 (tn ), μh1,n = μh1 (tn ), w h1,n = w n1 (tn ). Given the initial condition φh0,0 = φ0 , the first-order convex splitting scheme reads: for any 0 ≤ n ≤ N − 1, given φh0,n ∈ M h , find φh0,n+1 , μh0,n+1 , w h0,n+1 ∈ M h such that

(φh0,n+1 − φh0,n , ξh ) = −τ (∇ μh0,n+1 , ∇ξh ), 0 h,n+1 ,

(μ

ηh ) =(1 − ε

+ β)(φh0,n+1 ,

∀ξh ∈ M h ,

0 h ) − (∇ w h,n+1 , ∇ h )

η



η

  − β(φh0,n , ηn ) − 2(∇φh0,n , ∇ ηh ) + f φh0,n , ηh ,

( w h0,n+1 , ψh ) = −(∇φh0,n+1 , ∇ψh ),

∀ηh ∈ M h , ∀ψh ∈ M h .

Take the initial conditions φh1,0 = φ0 , μ0 = φ03 + (1 − ε )φ0 + 2φ0 +  w 0 , w 0 = φ0 . Then, we apply the SISDC method to get a second-order scheme: for any 0 ≤ n ≤ N − 1, given (φh1,n , μh1,n , w h1,n ) ∈ M h × M h × M h , find (φh1,n+1 , μh1,n+1 , w h1,n+1 ) ∈ M h × M h × M h such that for all (ξh , ηh , ψh ) ∈ M h × M h × M h

(φh1,n+1 , ξh ) =(φh1,n , ξh ) + τ G I (φh1,n+1 , μh1,n+1 , w h1,n+1 , ξh , ηh , ψh ) + τ G E (φh1,n , μh1,n , w h1,n , ξh , ηh , ψh ) −τ

G I (φh0,n+1 ,

0 0 h,n+1 , w h,n+1 , ξh ,

G E (φh0,n ,

(5.1)

0 0 h,n , w h,n , ξh ,

ηh , ψh ) − τ μ ηh , ψh ) τ + G (φh0,n+1 , μh0,n+1 , w h0,n+1 , ξh , ηh , ψh ) + G (φh0,n , μh0,n , w h0,n , ξh , ηh , ψh ), μ

τ

2

2

where

G I (φ, μ, w , ξ, η, ψ) = − (∇ μ, ∇ξ ) + (μ, η) − (1 − ε + β)(φ, η) + (∇ w , ∇ η) + ( w , ψ) + (∇φ, ∇ψ), G E (φ, μ, w , ξ, η, ψ) =β(φ, η) + 2(∇φ, ∇ η) − ( f (φ), η),

G = GI + GE.

We only need to solve linear systems above to get an approximate solution with the second-order of accuracy. It is easy-to-implement and efficient. 6. Numerical experiments In this section, we present some numerical experiments to validate that the scheme (5.1) is of second-order accuracy, mass-conservative, unconditionally energy stable. We also carry out some experiments to simulate the growth of a crystal and the phase transition. We use the P 1 finite element space for the variables φ, μ, w. All numerical experiments are performed by using the MATLAB software package iFEM [3]. 6.1. Convergence test Example 6.1. We test the convergence rates of the scheme (5.1) in space and time. The computational domain is set to be  = [0, 1]2 , and the parameters are taken to be ε = 0.05, M = 1, β = 3. We choose the suitable forcing term such that the exact solution is given by

φ(x, y , t ) = sin(t ) cos2 (2π x) cos2 (2π y ),

(x, y ) ∈  × [0, T ].

(6.1)

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

155

Table 1 The L 2 numerical errors and the H 1 numerical errors of the phase variable φ at T = 0.1. The exact solution is given in (6.1) with the parameters ε = 0.05, M = 1, β = 3, and the time step τ = 10−3 . h

φ − φh, N  L 2

order

φ − φh, N  H 1

order

1/8 1/16 1/32 1/64 1/128

5.1981e-02 1.7074e-02 4.6404e-03 1.1869e-03 2.9863e-04

— 1.61 1.88 1.97 1.99

3.4019e-01 1.4828e-01 6.2773e-02 2.9032e-02 1.4174e-02

— 1.20 1.24 1.11 1.03

Table 2 The L 2 numerical errors of the phase variable φ at T under the condition (A)

τ = h, T = 10

(B)

τ = Ch.

τ = 0.1h, T = 1

h

φ − φh, N  L 2

order

h

φ − φh, N  L 2

order

1/8 1/16 1/32 1/64 1/128

2.8127e-01 9.2793e-02 2.5272e-02 6.4708e-03 1.6288e-03

— 1.60 1.88 1.97 1.99

1/8 1/16 1/32 1/64 1/128

4.3785e-01 1.4399e-01 3.9137e-02 1.0010e-02 2.5179e-03

— 1.60 1.88 1.97 1.99

(C)

τ = 10h, T = 100

(D)

τ = 5h, T = 50

h

φ − φh, N  L 2

order

h

φ − φh, N  L 2

order

1/8 1/16 1/32 1/64 1/128

2.8294e-01 9.0109e-02 2.4013e-02 6.0692e-03 1.5130e-03

— 1.65 1.91 1.98 2.00

1/8 1/16 1/32 1/64 1/128

1.4897e-01 4.6941e-02 1.2464e-02 3.1458e-03 7.8359e-04

— 1.67 1.91 1.99 2.01

Table 3 The L 2 norm of φτ − φτ /2 at T and the rate of convergence in time.

τ

φτ − φτ /2  L 2

order

1 1/2 1/4 1/8 1/16

1.5303e-02 3.7657e-03 9.3819e-04 2.3485e-04 5.9274e-05

— 2.02 2.00 2.00 1.99

Firstly, we test the convergence order with respect to the spatial grid size h. The time step length τ should be chosen small enough so that the time discretization error is negligible compared with the spatial discretization error. Therefore, we take τ = 10−3 , the final time T = 0.1, N = τT . The L 2 errors and the H 1 errors of the phase variable φ between the exact

solution and the numerical solution are listed in Table 1, which shows the optimal convergence rates of P 1 element in L 2 and H 1 norms. Then, we test the convergence order with respect to the time step length τ . The strategy used in testing the spatial convergence order is a straight way, but taking h too small leads to a massive increasing in computational efforts. Here, we test the rate of convergence in time in another way. Note that, the L 2 numerical error of the phase variable φ is expected to be O (h2 ) + O (τ 2 ). Since we have validated that the spatial convergence order of φ in the L 2 norm approximate to 2, we can take τ = Ch so that O (h2 ) + O (τ 2 ) = O (h2 ). So, by studying the spatial L 2 numerical errors of φ at T under the condition τ = Ch, we could confirm the convergence order with respect to τ . The results of the L 2 numerical errors of φ at T under the condition τ = Ch are listed in Table 2. Without loss of generality, we take τ = 10h, τ = 5h, τ = h, τ = 0.1h, respectively. It’s clear that the rates of convergence suggested in Table 2 are consistent with our theoretical analysis. Furthermore, we provide another way to test the convergence order with respect to the time step length τ . Thenotation φτ −φ



φτ denotes the numerical approximate solution with the time step length τ . We use the formula log φτ /2 −φττ/2/4 L22 log(2) L to get the rate of convergence in time (cf. [20]). We take T = 5, h = 1/256. The results are listed in Table 3, which demonstrates that the convergence order with respect to the time step length

τ is 2.

156

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

Fig. 6.1. Time evolution of the energy functional and the mass for time steps of small differences in the energy evolution.

τ = 10, 5, 2, 1, 0.1, 0.01, 0.001 at [0, 30]. The small inset figures show the

6.2. Stability test Example 6.2. We try to demonstrate that the second-order scheme (5.1) is unconditionally energy stable and mass conservative. We use a domain  = [0, 128]2 , the parameters are taken to be ε = 0.05, M = 1, β = 3, the final time T = 30, and the spatial grid size h = 1. We take the initial condition

φ(x, y , 0) = cos(

π 32

x) cos(

π 32

y ),

(x, y ) ∈ .

We choose the approximation solution of τ = 0.001 as the exact solution. In Fig. 6.1, we plot the time evolution of the energy functional and the mass with different time step size of τ = 10, 5, 2, 1, 0.1, 0.01, 0.001 using the second-order scheme (5.1). From Fig. 6.1, one could find that the energy decreases over time, and when the time step gets small, such as τ = 1, 0.1, 0.01, the numerical solution matches well with the exact solution of τ = 0.001. We can also get that the mass conserves over time. Example 6.3. In this example, we try to verify that the scheme (5.1) with a random initial condition is still unconditionally energy stable and mass conservative. The computational domain is  = [0, 256]2 , the parameters are taken to be ε = 0.05, M = 1, β = 3, a random initial value around 0.5 is taken. We set the spatial grid size h = 1, the final time T = 60. Let τ = 10, 5, 2, 1, 0.1, 0.01, respectively, and we take the same random initial value among these different time steps. We plot the time evolution of the energy functional and the mass in Fig. 6.2, from which one could find that the energy is non-increasing and the mass is conserved. It demonstrates that the scheme proposed is unconditionally energy stable. 6.3. Crystal growth simulation and comparison with other schemes √π ]2 . Example 6.4. We simulate the growth of a crystal in a supercooled liquid in 2D. The computational domain is  = [0, 40 3

Similar numerical examples may be found in [13,30]. The initial seed crystal is represented with the following initial condition

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

157

Fig. 6.2. Time evolution of the energy functional and the mass with random initial conditions around 0.5 at different time steps The small inset figures show the small differences during the energy evolution.

τ = 10, 5, 2, 1, 0.1, 0.01.

 q 1 2q φ(x, y , 0) = φ + A ω(x, y ) cos( √ y ) cos(qx) − cos( √ y ) , 2

3

3

where

ω(x, y ) =

⎧ ⎪ ⎨

1−

(x−x0 )2 +( y − y 0 )2

⎪ ⎩0 ,

d20

2 , (x − x0 )2 + ( y − y 0 )2 ≤ d20 , otherwise.

Here O (x0 , y 0 ) is the center of the domain, d0 = 20√π , φ =

√

3 3

0.325 , 2

A=

√

12+8 6 φ, 15

and q =

√

3 . 2

We take the spatial grid size

h = 0.4 and take relatively small the time step length τ = 0.01 for better accuracy. The other parameters take the values ε = 0.5, M = 1, β = 3, T = 200. Fig. 6.3 shows snapshots of the numerical solution at different time t = 0, 10, 20, 30, 100, 200. We plot the time evolution of the energy functional and the mass in Fig. 6.4. It can be observed that the energy decreases at all times and the mass conserves during the process, which confirms that our algorithm is energy stable. Next, we test the efficiency of the SISDC scheme (5.1). For comparison, we take the convex splitting scheme (CS) in [18], and the scalar auxiliary approach (SAV) (see [25,26]). The CS scheme is

φn+1 − φn

τ

=μn+1/2 , 1

(6.2a)



μn+1/2 = (φn+1 + φn ) (φn+1 )2 + (φn )2 + w n+1/2 =

4 1 2

1−ε 2

(φn+1 + φn ) +  w n+1/2 + 3φn − φn−1 ,

(φn+1 + φn ) ,

with φ−1 = φ0 . The scheme is nonlinear. We take Newton iteration method to deal with the nonlinear term.

(6.2b) (6.2c)

158

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

Fig. 6.3. The simulation of the two-dimensional growth of a crystal in a supercooled liquid. Snapshots of the numerical solution φ are taken at t = 0, 10, 20, 30, 100, 200.

Fig. 6.4. Time evolution of the energy functional and the mass during the process of the growth of a crystal.

The SAV scheme is

φn+1 − φn

τ

=μn+1/2 ,

μn+1/2 = w n+1/2 + φn+1 + φn + w n+1/2 =

1

(6.3a) 1 2

(φn+1 + φn ) +

r n +1 + r n



2F 1 φ n+1/2

(φn+1 + φn ) ,   3  1   r n +1 − r n = φn+1/2 − ε φn+1/2 (φn+1 − φn ) dx, 2F 1 φ n+1/2 2





φn+1/2

3

 − ε φn+1/2 ,

(6.3b) (6.3c) (6.3d)

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

159

Fig. 6.5. The simulations of the two-dimensional growth of a crystal in a supercooled liquid using CS (row 1), SAV (row 2), SISDC (row 3). Snapshots of the numerical solution φ are taken at t = 10 (column 1), 30 (column 2), 100 (column 3), 200 (column 4).

! where φ n+1/2 =

3 φ 2 n

−

1 φ , 2 n−1

F 1 (φ) =

 (φ 2 −ε)2 

4

 + 1 dx. At the first step, we use a first-order SAV scheme to get the

numerical solution at t = τ2 , which is taken as φ 1/2 in (6.3) (cf. [25]). For the CS scheme, it needs to solve a nonlinear system at each time step. Newton iteration method is applied. At each time step, it needs to solve several variable coefficient equations. The SAV scheme and the SISDC scheme are linear and it involves solving two linear constant coefficient problems at each time step. Among all the three schemes, we take P 1 finite element space for the spatial discretization. We take the same parameters and initial condition with Example 6.4. In Fig. 6.5, we plot the simulations using the CS scheme (row 1), the SAV scheme (row 2) and the SISDC scheme (row 3) at t = 10 (column 1), t = 30 (column 2), t = 100 (column 3), and t = 200 (column 4). The evolutions of the energy functional and the mass are presented in Fig. 6.6. For either simulations in Fig. 6.5 or the evolutions in Fig. 6.6, no visible difference is observed. In Table 4, we present the elapsed time from t = 0 to t = 200 by using the three schemes. From Table 4, we can see that the computational time of the linear schemes, the SAV scheme and the SISDC scheme, is less than that of the nonlinear CS scheme. The computational time of the SISDC scheme is quite closed to that of the SAV scheme. It corresponds with the analysis for their algorithms, since they both need to solve two linear coefficient equations at each time step. It demonstrates that the SISDC scheme is efficient.

Remark 6.1. It is potential that the computational time of the SAV scheme and the SISDC scheme in Table 4 could be further reduced. From t = 0 to t = T , the matrix of the linear problems got from both schemes doesn’t vary at all numerical time steps. If the feature could be taken use of when executing computer programs, the computational efficiency of both schemes will be largely improved. On the other hand, for the CS scheme, the matrices of the equations vary in each iteration step. Limited by our current knowledge, we didn’t find an efficient way in MATLAB to take the feature into the implementing of the both schemes. We will learn more implementations of highly effective algorithms for solving large scale linear equations in MATLAB. Once a suitable way is available, the computational time of the SAV scheme and the SISDC scheme will be largely reduced.

160

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

Fig. 6.6. Time evolution of the energy functional and the mass during the process of the growth of a crystal using the CS scheme, the SAV scheme and the SISDC scheme. Table 4 Comparison of elapsed time during the process of crystal growth using the CS scheme, the SAV scheme and the SISDC scheme. Scheme

CS

SAV

SISDC

Elapsed Time

1.4018e+05s

5.7429e+04s

5.7670e+04s

6.4. Phase transition simulation Example 6.5. In this example, we simulate the phase transition behavior of the PFC model. The similar numerical examples can be found in [4,18,33]. Let the computation domain be  = [0, 128]2 , the initial condition be φ(x, y , 0) = 0.07 + 0.07 · rand(x, y ), where rand(x, y ) is a random number between −1 and 1. We take the spatial grid size h = 0.5 and the time step length τ = 0.1. The other parameters are ε = 0.025, M = 1, β = 3. We show the phase transition behavior of the density field φ in Fig. 6.7, the results are qualitatively consistent with the simulations in [18,33]. Fig. 6.8 shows the time evolution of the energy functional and the mass. It can be seen that the energy decreases at all times and the mass conserves over time, which demonstrates numerically that our scheme is energy stable. Just like the practices in [4,18], we research on the effect of using different time steps in the dynamics for both the second-order scheme (5.1) and the first-order scheme (3.6) during the simulation of the phase separation. The random initial condition and all other parameters are the same with those in Example 6.5. The simulations using the second-order scheme and the first-order scheme are shown in Fig. 6.9 and Fig. 6.10 respectively. For both schemes, simulations are presented with 4 different time steps: τ = 0.1 (row 1, Fig. 6.9 and Fig. 6.10), τ = 1 (row 2, Fig. 6.9 and Fig. 6.10), τ = 2 (row 3, Fig. 6.9 and Fig. 6.10), τ = 5 (row 1, Fig. 6.9 and Fig. 6.10). The first row in Fig. 6.9 are obtained at t = 250, 500, 750, 1000. We compare the simulations using different time steps at the same discrete energy levels, which means the results in each column in both Fig. 6.9 and Fig. 6.10 share the same discrete energy. The associated times are listed in the figures. First, we focus on the effect of using different time steps for each scheme. By comparing every column in Fig. 6.9 or in Fig. 6.10, we could find that using a smaller time step length needs a smaller numerical time to reach the same energy level. To compare the second-order scheme (5.1) and the first-order scheme (3.6), we can compare each row in Fig. 6.9 with the corresponding one in Fig. 6.10. It is easy to find that it takes less time for the second-order scheme than the first-order scheme to reach the same energy level, and the gap in time increases rapidly when the time step length becomes large. In Table 5, we present the elapsed time when executing codes. Each row in Table 5 show the computational time spent between t = 0 and the corresponding state in Fig. 6.9. From Table 5, we can find that using a bigger time step length costs less computational time to reach a corresponding energy level, though it costs a larger numerical time from Fig. 6.9. The reason may be that a smaller step length brings more steps. Based on this observation, it is reasonable to take a large time step length to save actual time. However, the evolution of the energy in Fig. 6.8 shows that the energy functional has a rapid decrease in early time. A large time step may fail to capture this character correctly. So, an adaptive time step strategy is desirable for a long-time simulation. The relevant research is currently under study. 7. Concluding remarks In this paper, a linear, unconditionally energy stable and second-order (in time) numerical scheme for the PFC equation is proposed. The scheme is based on a linear convex splitting scheme and the SISDC method. We prove that the convex splitting scheme is mass conservative, uniquely solvable, unconditionally energy stable and of first-order accuracy in time. Then the SISDC method is applied to improve the order of accuracy. The resulted second-order scheme inherits the advantages of the first-order convex splitting scheme and thus leads to linear equations at each time step, so it is simple and

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

161

Fig. 6.7. The simulation of the phase separation. Snapshots of the numerical approximation of φ are taken at t =10, 200, 500, 800, 1000, 2000.

Fig. 6.8. Time evolution of the energy functional and the mass during the process of the phase separation. The small inset figure changes the value of ordinate axis to make the evolution more visible.

Table 5 The elapsed time corresponding to states in Fig. 6.9. Each row presents the elapsed time from t = 0 to the associate state in Fig. 6.9.

τ τ τ τ

= 0.1 =1 =2 =5

column 1 in Fig. 6.9

column 2 in Fig. 6.9

column 3 in Fig. 6.9

column 4 in Fig. 6.9

3.6354e+04s 1.0381e+04s 1.0686e+04s 9.1979e+03s

7.2860e+04s 2.0931e+04s 2.1049e+04s 1.8014e+04s

1.0940e+05s 3.3045e+04s 3.1914e+04s 2.8453e+04s

1.4598e+05s 4.4883e+04s 4.3001e+04s 3.9411e+04s

162

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

Fig. 6.9. The simulation of the phase separation with the second-order scheme (5.1) using four different time steps: τ = 0.1 (row 1), τ = 1 (row 2), τ = 2 (row 3), τ = 5 (row 4). The energy associate with each column is equal. The associate numerical time of each snapshot is shown on the bottom of the plot.

efficient. Besides, the second-order scheme is a one-step method, which makes it easier to analyze and implement than the usual two-step methods. And we prove rigorously that the resulted scheme is unconditionally weak energy stable and second-order accurate. In section 6, numerical experiments are presented to verify that the resulted scheme is of secondorder accuracy, energy stable, mass conservative. We also present a few numerical results from some standard numerical tests of the PFC model, which demonstrate that the scheme proposed performs well. Acknowledgements The authors would like to thank the reviewers for their valuable comments and suggestions. This work is partially supported by NSFC grant No. 11571274.

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164

163

Fig. 6.10. The simulation of the phase separation with the first-order scheme (3.6) using four different time steps: τ = 0.1 (row 1), τ = 1 (row 2), τ = 2 (row 3), τ = 5 (row 4). The energy associate with each column is equal to that in Fig. 6.9. The associate numerical time of each snapshot is shown on the bottom of the plot.

Appendix A. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.apnum.2019.01.017. References [1] A. Baskaran, J.S. Lowengrub, C. Wang, S.M. Wise, Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation, SIAM J. Numer. Anal. 51 (5) (2013) 2851–2873.

164

[2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]

S. Pei et al. / Applied Numerical Mathematics 140 (2019) 134–164 L.A. Caffarelli, N.E. Muler, An l∞ bound for solutions of the Cahn–Hilliard equation, Arch. Ration. Mech. Anal. 133 (2) (Dec 1995) 129–144. L. Chen, iFEM: An Integrated Finite Element Methods Package in MATLAB, Technical report, University of California at Irvine, 2009. M. Cheng, J.A. Warren, An efficient algorithm for solving the phase field crystal model, J. Comput. Phys. 227 (12) (2008) 6241–6248. N. Condette, C. Melcher, E. Süli, Spectral approximation of pattern-forming nonlinear evolution equations with double-well potentials of quadratic growth, Math. Comput. 80 (273) (2011) 205–223. L. Dong, W. Feng, C. Wang, S.M. Wise, Z. Zhang, Convergence analysis and numerical implementation of a second order numerical scheme for the three-dimensional phase field crystal equation, Comput. Math. Appl. 75 (6) (2018) 1912–1928. A. Dutt, L. Greengard, V. Rokhlin, Spectral deferred correction methods for ordinary differential equations, BIT 40 (2) (2000) 241–266. K.R. Elder, M. Grant, Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals, Phys. Rev. E 70 (2004) 051605. K.R. Elder, M. Katakowski, M. Haataja, M. Grant, Modeling elasticity in crystal growth, Phys. Rev. Lett. 88 (2002) 245701. L.C. Evans, Partial Differential Equations, second edition, American Mathematical Society, 2010. D.J. Eyre, Unconditionally gradient stable time marching the Cahn–Hilliard equation, MRS Proc. 529 (1998) 39. X. Feng, T. Tang, J. Yang, Long time numerical simulations for phase-field problems using p-adaptive spectral deferred correction methods, SIAM J. Sci. Comput. 37 (1) (2015) A271–A294. K. Glasner, S. Orizaga, Improving the accuracy of convexity splitting methods for gradient flow equations, J. Comput. Phys. 315 (2016) 52–64. H. Gomez, X. Nogueira, An unconditionally energy-stable method for the phase field crystal equation, Comput. Methods Appl. Mech. Eng. 249–252 (2012) 52–61. Z. Guan, C. Wang, S.M. Wise, A convergent convex splitting scheme for the periodic nonlocal Cahn–Hilliard equation, Numer. Math. 128 (2) (2014) 377–406. R. Guo, Y. Xu, Local discontinuous Galerkin method and high order semi-implicit scheme for the phase field crystal equation, SIAM J. Sci. Comput. 38 (1) (2016) A105–A127. J.G. Heywood, R. Rannacher, Finite-element approximation of the nonstationary Navier–Stokes problem. Part IV: error analysis for second-order time discretization, SIAM J. Numer. Anal. 27 (2) (1990) 353–384. Z. Hu, S. Wise, C. Wang, J. Lowengrub, Stable and efficient finite-difference nonlinear-multigrid schemes for the phase field crystal equation, J. Comput. Phys. 228 (15) (2009) 5323–5339. D. Kessler, R.H. Nochetto, A. Schmidt, A posteriori error control for the Allen–Cahn problem: circumventing Gronwall’s inequality, ESAIM: Math. Model. Numer. Anal. 38 (1) (2004) 129–142. S. Li, H. Zheng, W.J. Layton, A decoupling method with different subdomain time steps for the nonstationary Stokes–Darcy model, Numer. Methods Partial Differ. Equ. 29 (2) (2013) 549–583. M.L. Minion, Semi-implicit projection methods for incompressible flow based on spectral deferred corrections, Appl. Numer. Math. 48 (3) (2004) 369–387. N. Provatas, J.A. Dantzig, B. Athreya, P. Chan, P. Stefanovic, N. Goldenfeld, K.R. Elder, Using the phase-field crystal method in the multi-scale modeling of microstructure evolution, JOM 59 (7) (2007) 83–90. J. Shen, T. Tang, L.-L. Wang, Spectral Methods: Algorithms, Analysis and Applications, Springer Series in Computational Mathematics, vol. 41, Springer, Berlin, Heidelberg, 2011. J. Shen, C. Wang, X. Wang, S.M. Wise, Second-order convex splitting schemes for gradient flows with Ehrlich–Schwoebel type energy: application to thin film epitaxy, SIAM J. Numer. Anal. 50 (1) (2012) 105–125. J. Shen, J. Xu, J. Yang, A new class of efficient and robust energy stable schemes for gradient flows, arXiv e-prints, 2017. J. Shen, J. Xu, J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, J. Comput. Phys. 353 (2018) 407–416. J. Shen, X. Yang, Numerical approximations of Allen–Cahn and Cahn–Hilliard equations, Discrete Contin. Dyn. Syst. 28 (2010) 1669–1691. J. Shin, H.G. Lee, J.-Y. Lee, First and second order numerical methods based on a new convex splitting for phase-field crystal equation, J. Comput. Phys. 327 (2016) 519–542. J. Swift, P.C. Hohenberg, Hydrodynamic fluctuations at the convective instability, Phys. Rev. A 15 (1977) 319–328. P. Vignal, L. Dalcin, D. Brown, N. Collier, V. Calo, An energy-stable convex splitting for the phase-field crystal equation, Comput. Struct. 158 (2015) 355–368. S.M. Wise, C. Wang, J.S. Lowengrub, An energy-stable and convergent finite-difference scheme for the phase field crystal equation, SIAM J. Numer. Anal. 47 (3) (2009) 2269–2288. C. Xu, T. Tang, Stability analysis of large time-stepping methods for epitaxial growth models, SIAM J. Numer. Anal. 44 (4) (2006) 1759–1779. X. Yang, D. Han, Linearly first- and second-order, unconditionally energy stable schemes for the phase field crystal model, J. Comput. Phys. 330 (2017) 1116–1134. X. Yang, L. Ju, Efficient linear schemes with unconditional energy stability for the phase field elastic bending energy model, Comput. Methods Appl. Mech. Eng. 315 (2017) 691–712. Z. Zhang, Y. Ma, Z. Qiao, An adaptive time-stepping strategy for solving the phase field crystal model, J. Comput. Phys. 249 (2013) 204–215.