Mesh-robustness of the variable steps BDF2 method for the Cahn-Hilliard model
MMesh-robustness of the variable steps BDF2 methodfor the Cahn-Hilliard model
Hong-lin Liao ∗ Bingquan Ji † Lin Wang ‡ Zhimin Zhang § November 19, 2020
Abstract
The two-step backward differential formula (BDF2) implicit method with unequal time-steps is investigated for the Cahn-Hilliard model by focusing on the numerical influencesof time-step variations. The suggested method is proved to preserve a modified energydissipation law at the discrete levels if the adjoint time-step ratios fulfill a new step-ratiorestriction 0 < r k := τ k /τ k − ≤ r user ( r user can be chosen by the user such that r user < . L norm error estimate is established at the first time under the updated step-ratio restriction0 < r k ≤ r user . The time-stepping scheme is mesh-robustly convergent in the sense that theconvergence constant (prefactor) in the error estimate is independent of the adjoint time-stepratios. On the basis of ample tests on random time meshes, a useful adaptive time-steppingstrategy is applied to efficiently capture the multi-scale behaviors and to accelerate the long-time simulation approaching the steady state. Keywords : Cahn-Hilliard model; adaptive BDF2 method; discrete energy dissipation law;orthogonal convolution kernels, discrete convolution embedding inequality; error estimate
AMS subject classifications.
The Cahn-Hilliard (CH) model is an efficient approach to describe the coarsening dynamics ofa binary alloy system [4] and has been applied in other fields including image inpainting [3] andtumor growth [10]. Consider a free energy functional of Ginzburg–Landau type, E [Φ] = (cid:90) Ω (cid:104) (cid:15) |∇ Φ | + F (Φ) (cid:105) d x with F (Φ) := 14 (Φ − (1.1) ∗ ORCID 0000-0003-0777-6832; Department of Mathematics, Nanjing University of Aeronautics and Astronau-tics, Nanjing 211106, P. R. China. Hong-lin Liao ([email protected] and [email protected]) is supported by agrant 12071216 from National Natural Science Foundation of China. † Department of Mathematics, Nanjing University of Aeronautics and Astronautics, 211101, P. R. China.Bingquan Ji ([email protected]). ‡ Beijing Computational Science Research Center, Beijing 100193, P. R. China. E-mail: [email protected]. § Beijing Computational Science Research Center, Beijing 100193, P. R. China; and Department of Mathe-matics, Wayne State University, Detroit, MI 48202, USA. E-mails: [email protected] and [email protected] author is supported in part by the NSFC grant 11871092 and NSAF grant U1930402. a r X i v : . [ m a t h . NA ] F e b here x ∈ Ω ⊆ R and 0 < (cid:15) < H − gradient flow associatedwith the free energy functional E [ φ ], ∂ t Φ = κ ∆ µ with µ := δEδ Φ = F (cid:48) (Φ) − (cid:15) ∆Φ , (1.2)where the parameter κ is the mobility related to the characteristic relaxation time of systemand µ is the chemical potential. Assume that Φ is periodic over the domain Ω. By applying theintegration by parts, one can find the volume conservation, (cid:0) Φ( t ) , (cid:1) = (cid:0) Φ( t ) , (cid:1) , (1.3)and the following energy dissipation law,d E d t = (cid:0) δEδ Φ , ∂ t Φ (cid:1) = ( µ, ∆ µ ) = − (cid:107)∇ µ (cid:107) ≤ , (1.4)where ( u, v ) := (cid:82) Ω uv d x , and the associated L norm (cid:107) v (cid:107) = (cid:112) ( v, v ) for all u, v ∈ L (Ω).The main aim of this paper is to present a rigorous stability and convergence analysis of theBDF2 implicit method with variable time-steps for simulating the CH model (1.2). Considerthe nonuniform time levels 0 = t < t < · · · < t N = T with the time-step sizes τ k := t k − t k − for 1 ≤ k ≤ N , and denote the maximum time-step size τ := max ≤ k ≤ N τ k . Let the adjointtime-step ratio r k := τ k /τ k − for 2 ≤ k ≤ N . Our analysis will focus on the influence of non-uniform time grids (with the associated time-step ratios) on the numerical solution by carefullyevaluating the stability and convergence.This is motivated by the following facts: • The BDF2 method is A-stable and L-stable such that it would be more suitable thanCrank-Nicolson type schemes for solving the stiff dissipative problems, see e.g. [5, 6, 29]. • The nonuniform grid and adaptive time-stepping techniques [13, 18–20, 24] are powerfulin capturing the multi-scale behaviors and accelerating the long-time simulations of phasefield models including the CH model. • The convergence theory of variable-steps BDF2 scheme remains incomplete for nonlinearparabolic equations. Actually, the required step-ratio constraint for the L norm stabilityare severer than the classical zero-stability condition r k < √
2, given by Grigorieff [12].Always, they contain some undesirable pre-factors C r exp( C r Γ n ) or C r exp( C r t n ), see e.g.[2,8,11,17,28], where Γ n may be unbounded when certain time-step variations appear and C r may be infinity as the step-ratios approach the zero-stability limit 1 + √ discrete orthogonal convolution (DOC) kernels was suggested to verify that, if the adjoint time-step ratios r k satisfy a newzero-stability condition [22], that is, 0 < r k < (cid:0) √ (cid:1) / ≈ . S0 . 0 < r k ≤ r user ( < . ≤ k ≤ N ,where the value of r user can be chosen in adaptive time-stepping computations by the user suchthat r user < . r user = 2 , S0 , we will present an L norm error estimate with an improved prefactor, see Theorem 4.1, C φ exp (cid:0) c (cid:15) t n − (cid:1) . Here and hereafter, any subscripted C , such as C u and C φ , denotes a generic positive constant,not necessarily the same at different occurrences; while, any subscripted c , such as c (cid:15) , c Ω , c p , c z and so on, denotes a fixed constant. The appeared constants may be dependent on the givendata (typically, the interface width parameter (cid:15) ) and the solution but are always independentof the spatial lengths, the time t n , the step sizes τ n and the step ratios r n . It is interesting toemphasize that, under the step-ratio constraint S0 , the involved constants are bounded evenwhen the step-ratios r n approach the user limit r user such that the variable-steps BDF2 schemeis mesh-robustly convergent.To the best of our knowledge, this is the first time such an optimal L norm error estimate ofvariable-steps BDF2 method is established for the Cahn-Hiliard (and Allen-Cahn) type models.As a closely related work, the variable-steps BDF2 scheme for the Allen-Chan equation was alsoinvestigated in [20] by using the discrete complementary convolution kernels. The BDF2 schemewas proved to preserve the maximum bound principle if the step-ratios satisfy the classicalzero-stability condition r k < √
2. The maximum norm error estimate with a prefactor − η exp( t n − η ) was obtained, where the parameter η → r k → √
2. It is to mentionthat, under the constraint S0 , one can follow the present analysis to obtain a new L norm errorestimate that is robustly stable to the variations of time-steps.Given a grid function { v k } Nk =0 , put (cid:79) τ v k := v k − v k − , ∂ τ v k := (cid:79) τ v k /τ k for k ≥
1. Taking v n = v ( t n ), we view the variable-steps BDF2 formula as a discrete convolution summation D v n := n (cid:88) k =1 b ( n ) n − k (cid:79) τ v k for n ≥ , (1.5)where the discrete convolution kernels b ( n ) n − k are defined for n ≥ b ( n )0 := 1 + 2 r n τ n (1 + r n ) , b ( n )1 := − r n τ n (1 + r n ) and b ( n ) j := 0 for j ≥ . (1.6)Without losing the generality, assume that an accurate solution φ is available. We consider thestability and convergence of the BDF2 scheme for solving the CH equation (1.2): D φ n = κ ∆ h µ n with µ n := F (cid:48) ( φ n ) − (cid:15) ∆ h φ n for 2 ≤ n ≤ N (1.7)subject to the periodic boundary conditions. The spatial operators are approximated by theFourier pseudo-spectral method, as described in the next section.3he unique solvability of the BDF2 scheme (1.7) is established in Theorem 2.1 by using thefact that the solution of nonlinear scheme (1.7) is equivalent to the minimization of a convexfunctional. Lemma 2.1 shows that the BDF2 convolution kernels b ( n ) n − k are positive definiteprovided the adjacent time-step ratios r k satisfy the sufficient condition S0 . Then we verifyin Theorem 2.2 that the adaptive BDF2 time-stepping method (1.7) has a modified energydissipation law at the discrete levels under a proper step-size constraint.We are to emphasize that the solution estimates in section 2 are based on the standard form(1.7) of the BDF2 implicit scheme, but in the subsequent L norm error analysis we will use anequivalent convolution form with a class of discrete orthogonal convolution (DOC) kernels. TheDOC kernels { θ ( n ) n − k } nk =2 are defined by (this definition is slightly different from those in [18,19,22]since we do not introduce the discrete kernel b (1)0 for the first-level solver) θ ( n )0 := 1 b ( n )0 for n ≥ θ ( n ) n − k := − b ( k )0 n (cid:88) j = k +1 θ ( n ) n − j b ( j ) j − k for n ≥ k + 1 ≥ . (1.8)One has the following discrete orthogonal identity n (cid:88) j = k θ ( n ) n − j b ( j ) j − k ≡ δ nk for 2 ≤ k ≤ n , (1.9)where δ nk is the Kronecker delta symbol. By exchanging the summation order and using theidentity (1.9), it is not difficult to check that n (cid:88) j =2 θ ( n ) n − j D v j = n (cid:88) j =2 θ ( n ) n − j b ( j ) j − (cid:79) τ v + n (cid:88) j =2 θ ( n ) n − j j (cid:88) (cid:96) =2 b ( j ) j − (cid:96) (cid:79) τ v (cid:96) = θ ( n ) n − b (2)1 (cid:79) τ v + (cid:79) τ v n for n ≥
2. (1.10)Acting the DOC kernels θ ( m ) m − n on the first equation in (1.7) and summing n from n = 2 to m ,we apply (1.10) to find the equivalent convolution form (replacing m by n ) (cid:79) τ φ n = − θ ( n ) n − b (2)1 (cid:79) τ φ + κ n (cid:88) j =2 θ ( n ) n − j ∆ h µ j for 2 ≤ n ≤ N . (1.11)Note that, by following the proof of [21, Lemma 2.1], we have m (cid:88) j = k b ( m ) m − j θ ( j ) j − k ≡ δ mk for 2 ≤ k ≤ m . (1.12)With the help of this mutually orthogonal identity, one can recover the original form (1.7) byacting the BDF2 kernels b ( m ) m − n on the new formulation (1.11). In this sense, the DOC kernelsdefine a reversible discrete transform between (1.7) and the convolution form (1.11).To perform the L norm error analysis, section 3 presents some properties of the DOC kernels θ ( n ) n − k and some new convolution embedding inequalities with respect to the DOC kernels, seeLemmas 3.1–3.9. By making use of the H norm solution bound obtained in Lemma 2.2, weestablish an optimal L norm error estimate in section 4. Numerical tests and comparisons arepresented in section 5 to validate the accuracy and effectiveness of the BDF2 method (1.7),especially when coupled with an adaptive stepping strategy.4 Solvability and energy dissipation law
We use the same spatial notations in [18]. Set the space domain Ω = (0 , L ) and consider theuniform length h x = h y = h := L/M in each direction for an even positive integer M . LetΩ h := (cid:8) x h = ( ih, jh ) | ≤ i, j ≤ M (cid:9) and put ¯Ω h := Ω h ∪ ∂ Ω. Denote the space of L -periodicgrid functions V h := { v | v = ( v h ) is L -periodic for x h ∈ ¯Ω h } . For a periodic function v ( x ) on ¯Ω, let P M : L (Ω) → F M be the standard L projectionoperator onto the space F M , consisting of all trigonometric polynomials of degree up to M/ I M : L (Ω) → F M be the trigonometric interpolation operator [26], i.e.,( P M v ) ( x ) = M/ − (cid:88) (cid:96),m = − M/ (cid:98) v (cid:96),m e (cid:96),m ( x ) , ( I M v ) ( x ) = M/ − (cid:88) (cid:96),m = − M/ (cid:101) v (cid:96),m e (cid:96),m ( x ) , where the complex exponential basis function e (cid:96),m ( x ) := e i ν ( (cid:96)x + my ) with ν = 2 π/L . The coef-ficients (cid:98) v (cid:96),m refer to the standard Fourier coefficients of function v ( x ), and the pseudo-spectralcoefficients (cid:101) v (cid:96),m are determined such that ( I M v ) ( x h ) = v h .The Fourier pseudo-spectral first and second order derivatives of v h are given by D x v h := M/ − (cid:88) (cid:96),m = − M/ (i ν(cid:96) ) (cid:101) v (cid:96),m e (cid:96),m ( x h ) , D x v h := M/ − (cid:88) (cid:96),m = − M/ (i ν(cid:96) ) (cid:101) v (cid:96),m e (cid:96),m ( x h ) . The differentiation operators D y and D y can be defined in the similar fashion. In turn, we candefine the discrete gradient and Laplacian in the point-wise sense, respectively, by ∇ h v h := ( D x v h , D y v h ) T and ∆ h v h := ∇ h · ( ∇ h v h ) = (cid:0) D x + D y (cid:1) v h For any grid functions v, w ∈ V h , define the discrete inner product (cid:104) v, w (cid:105) := h (cid:80) x h ∈ Ω h v h w h ,and the associated L norm (cid:107) v (cid:107) := (cid:107) v (cid:107) l = (cid:112) (cid:104) v, v (cid:105) . Also, we will use the discrete l q norm (cid:107) v (cid:107) l q := q (cid:113) h (cid:80) x h ∈ Ω h | v h | q and the H seminorm (cid:13)(cid:13) ∇ h v (cid:13)(cid:13) := (cid:113) h (cid:80) x h ∈ Ω h |∇ h v h | . It is easyto check the discrete Green’s formulas, (cid:104)− ∆ h v, w (cid:105) = (cid:104)∇ h v, ∇ h w (cid:105) and (cid:10) ∆ h v, w (cid:11) = (cid:104) ∆ h v, ∆ h w (cid:105) ,see [26] for more details. Also we have the following discrete embedding inequality simulatingthe Sobolev embedding H (Ω) (cid:44) → L (Ω), (cid:13)(cid:13) v (cid:13)(cid:13) l ≤ c Ω (cid:0)(cid:13)(cid:13) v (cid:13)(cid:13) + (cid:13)(cid:13) ∇ h v (cid:13)(cid:13)(cid:1) for any v ∈ V h . (2.1)For the underlying volume-conservative problem, it is also to define a mean-zero functionspace ˚ V h := (cid:8) v ∈ V h | (cid:104) v, (cid:105) = 0 (cid:9) ⊂ V h . As usual, following the arguments in [26], one canintroduce an discrete version of inverse Laplacian operator ( − ∆ h ) − γ as follows. For a gridfunction v ∈ ˚ V h , define( − ∆ h ) − γ v h := M/ − (cid:88) (cid:96), m = − M/ (cid:96), m ) (cid:54) = (cid:0) ν (cid:0) (cid:96) + m (cid:1)(cid:1) − γ (cid:101) v (cid:96),m e (cid:96),m ( x h ) , and an H − inner product (cid:104) v, w (cid:105) − := (cid:10) ( − ∆ h ) − v, w (cid:11) . The associated H − norm (cid:107)·(cid:107) − can bedefined by (cid:107) v (cid:107) − := (cid:113) (cid:104) v, v (cid:105) − . We have the following Poincar´e type inequality with the usual5oincar´e constant c p , (cid:13)(cid:13) v (cid:13)(cid:13) − ≤ c p (cid:13)(cid:13) v (cid:13)(cid:13) , and the generalized H¨older inequality, (cid:13)(cid:13) v (cid:13)(cid:13) ≤ (cid:13)(cid:13) ∇ h v (cid:13)(cid:13)(cid:13)(cid:13) v (cid:13)(cid:13) − for any v ∈ ˚ V h . (2.2)Also the discrete embedding inequality (2.1) can be simplified as ( c z = c Ω + c Ω c p ) (cid:13)(cid:13) v (cid:13)(cid:13) l ≤ c z (cid:13)(cid:13) ∇ h v (cid:13)(cid:13) for any v ∈ ˚ V h . (2.3) To focus on the numerical analysis of the BDF2 solution, it is to assume that A1 . A starting scheme is properly chosen to compute an accurate first-level solution φ suchthat it preserves the initial volume, (cid:10) φ , (cid:11) = (cid:10) φ , (cid:11) = (cid:10) P M Φ , (cid:11) , and also preservescertain (maybe, modified) energy dissipation law. There exists a positive constant c ,depended on the domain Ω, the interface parameter (cid:15) and the initial value φ , such that (cid:15) (cid:13)(cid:13) ∇ h φ (cid:13)(cid:13) + (cid:10) F ( φ ) , (cid:11) + τ κ (cid:13)(cid:13) ∂ τ φ (cid:13)(cid:13) − ≤ c . Remark 1.
Assumption A1 can be satisfied by many of first-level solvers. The BDF1 schemewould be suited for computing a second-order solution φ ; however, a very small initial step τ would not be suggested here since it arrives at a large step-ratio r and eventually affects theaccuracy of solution in the whole simulation, see numerical results in [23].The Crank-Nicolson scheme at the first time-level can generate a second-order differencequotient ∂ τ φ ; but a very small initial step τ would not be suggested either because it would beprone to generate nonphysical oscillations. To control possibly initial oscillations, we suggesta special step-ratio r = √ / in the implementation of our scheme (1.7) . Actually, by taking φ γ := φ , φ := φ , τ ∗ := τ + τ and γ := τ /τ ∗ with r = 1 /γ − , the first two steps of (1.7) are equivalent to the following TR-BDF2 method φ γ − φ γτ ∗ = κ h µ γ + κ h µ , − γ (1 − γ ) τ ∗ φ − γ (1 − γ ) τ ∗ φ γ + 1 − γγτ ∗ φ = κ ∆ h µ , which was shown to be L-stable for γ = 2 − √ , see [15, 27]. Generally, the TR-BDF2 methodis defined as a one-step method resulting from the composition of the trapezoidal rule in thefirst substep, followed by BDF2 formula in the second substep. This combination is empiricallyjustified under the rationale of combining the good accuracy of the trapezoidal rule with thestability and damping of fast modes guaranteed by BDF2 scheme.The user can also choose other self-starting L-stable schemes, such as the second-order singly-diagonal implicit Runge-Kutta method [1, 14] for α = (2 − √ / , φ α − φ ατ = κ ∆ h µ α , φ − φ τ = ακ ∆ h µ + (1 − α ) κ ∆ h µ α . Under the assumption A1 , the solution φ n of BDF2 scheme (1.7) preserves the volume, (cid:10) φ n , (cid:11) = (cid:10) φ , (cid:11) for n ≥
2. Actually, taking the inner product of (1.7) by 1 and applying the6ummation by parts, one can check that (cid:10) D φ j , (cid:11) = 0 for j ≥
2. Multiplying both sides of thisequality by the DOC kernels θ ( n ) n − j and summing the index j from j = 2 to n , we get n (cid:88) j =2 θ ( n ) n − j (cid:10) D φ j , (cid:11) = 0 for n ≥ . It leads to (cid:10) (cid:79) τ φ n , (cid:11) = 0 directly by taking v j = φ j in the equality (1.10). Simple inductionyields the volume conversation law, (cid:10) φ n , (cid:11) = (cid:10) φ n − , (cid:11) = · · · = (cid:10) φ , (cid:11) for n ≥ Theorem 2.1. If A1 holds and the time-step size τ n ≤ (cid:15) (1+2 r n ) κ (1+ r n ) , the variable-steps BDF2time-stepping scheme (1.7) is uniquely solvable.Proof. For any fixed time-level indexes n ≥
2, consider the following energy functional G on thespace V ∗ h := (cid:8) z ∈ V h | (cid:10) z, (cid:11) = (cid:10) φ n − , (cid:11)(cid:9) ,G [ z ] := 12 b ( n )0 (cid:13)(cid:13) z − φ n − (cid:13)(cid:13) − + b ( n )1 (cid:10) (cid:79) τ φ n − , z − φ n − (cid:11) − + κ(cid:15) (cid:13)(cid:13) ∇ h z (cid:13)(cid:13) + κ (cid:10) F ( z ) , (cid:11) . (2.4)Under the time-step condition τ n ≤ (cid:15) (1+2 r n ) κ (1+ r n ) or b ( n )0 (cid:15) ≥ κ , the functional G is strictly convexsince, for any λ ∈ R and any ψ ∈ ˚ V h ,d G d λ [ z + λψ ] (cid:12)(cid:12)(cid:12) λ =0 = b ( n )0 (cid:13)(cid:13) ψ (cid:13)(cid:13) − + κ(cid:15) (cid:13)(cid:13) ∇ h ψ (cid:13)(cid:13) + 3 κ (cid:13)(cid:13) zψ (cid:13)(cid:13) − κ (cid:13)(cid:13) ψ (cid:13)(cid:13) ≥ (cid:13)(cid:13) ∇ h ψ (cid:13)(cid:13)(cid:13)(cid:13) ψ (cid:13)(cid:13) − − κ (cid:13)(cid:13) ψ (cid:13)(cid:13) + 3 κ (cid:13)(cid:13) zψ (cid:13)(cid:13) > , where the generalized H¨older inequality (2.2) has been applied in the last step. Thus thefunctional G has a unique minimizer, denoted by φ n , if and only if it solves the equation0 = d G d λ [ z + λψ ] (cid:12)(cid:12)(cid:12) λ =0 = (cid:10) b ( n )0 ( z − φ n − ) + b ( n )1 (cid:79) τ φ n − , ψ (cid:11) − + κ (cid:10) − (cid:15) ∆ h z + F (cid:48) ( z ) , ψ (cid:11) = (cid:68) b ( n )0 ( z − φ n − ) + b ( n )1 (cid:79) τ φ n − − κ ∆ h (cid:0) F (cid:48) ( z ) − (cid:15) ∆ h z (cid:1) , ψ (cid:69) − . This equation holds for any ψ ∈ ˚ V h if and only if the unique minimizer φ n ∈ V ∗ h solves b ( n )0 ( φ n − φ n − ) + b ( n )1 (cid:79) τ φ n − − κ ∆ h (cid:0) F (cid:48) ( φ n ) − (cid:15) ∆ h φ n (cid:1) = 0 , which is just the BDF2 scheme (1.7). It completes the proof.The proof of Theorem 2.1 also says that the solution of the BDF2 scheme (1.7) is equivalentto the minimization of a convex functional G [ z ] under the condition τ n ≤ (cid:15) (1+2 r n ) κ (1+ r n ) . We seethat the BDF2 implicit scheme is also convex according to Xu et al. [29]. In our previous work [22, Lemma 2.1], the BDF2 kernels b ( n ) n − k are shown to be positive definiteif the adjacent time-step ratios 0 < r k < √ . The following result shows that this sufficientcondition can be further improved in the theoretical manner. This improvement is inspired7y [18, LemmaA.1] to find a lower bound for the eigenvalues of the step-scaled matrix (cid:101) B , seeLemma 3.2 below. For simplicity, we denote R L ( z, s ) := 2 + 4 z − z / z − s / s for 0 ≤ z, s ≤ r ∗ , (2.5)where r ∗ ≈ .
864 is the positive root of the equation 1 + 2 r ∗ − r / ∗ = 0. According to the proofof [18, LemmaA.1], R L (cid:0) z, s (cid:1) is increasing in (0 ,
1) and decreasing in (1 , r ∗ ) with respect to z .Also, R L ( z, s ) is decreasing with respect to s such that R L ( z, s ) > min { R L (0 , r ∗ ) , R L (cid:0) r ∗ , r ∗ (cid:1) } = 2(1 + 2 r ∗ − r / ∗ )1 + z = 0 for 0 < z, s < r ∗ . Lemma 2.1.
Let < r k < . for ≤ k ≤ N . For any real sequence { w k } nk =1 with n entries,it holds that w k k (cid:88) j =1 b ( k ) k − j w j ≥ r / k +1 r k +1 w k τ k − r / k r k w k − τ k − + R L ( r k , r k +1 ) w k τ k for k ≥ .So the discrete convolution kernels b ( n ) n − k are positive definite in the sense that n (cid:88) k =2 w k k (cid:88) j =2 b ( k ) k − j w j ≥ n (cid:88) k =2 R L ( r k , r k +1 ) w k τ k for n ≥ . Proof.
Applying the inequality − ab ≥ − a − b , we take u k := w k / √ τ k to find2 w k k (cid:88) j =1 b ( k ) k − j w j = 2 τ k b ( k )0 u k + 2 √ τ k τ k − b ( k )1 u k u k − ≥ r k r k u k − r / k r k (cid:0) u k + u k − (cid:1) = r / k +1 r k +1 w k τ k − r / k r k w k − τ k − + R L ( r k , r k +1 ) w k τ k for k ≥ k = 2 to n , it is straightforward to obtain the claimed secondinequality and complete the proof. Remark 2.
This lemma updates the sufficient condition of [22, Lemma 2.1]. Thus by followingthe discussions in [22, Remark 3 and Remark 5], one can verify that the variable-step BDF2method is A-stable and zero-stable if the adjoint time-step ratios r k satisfy the updated sufficientcondition: < r k < . for ≤ k ≤ N . Let E [ φ k ] be the discrete version of free energy functional (1.1), given by E [ φ k ] := (cid:15) (cid:13)(cid:13) ∇ h φ k (cid:13)(cid:13) + (cid:10) F ( φ k ) , (cid:11) for k ≥
1. (2.6)Next theorem shows that the numerical scheme (1.7) preserves a modified energy dissipationproperty at the discrete levels, and it is mesh-robustly stable in an energy norm if S0 holds.8 heorem 2.2. Let S0 holds. If the time-step sizes are properly small such that τ n ≤ (cid:15) κ min (cid:110) r n r n , R L ( r n , r n +1 ) (cid:111) for n ≥ , (2.7) the BDF2 scheme (1.7) preserves the following energy dissipation law E [ φ n ] ≤ E [ φ n − ] ≤ E [ φ ] for n ≥ ,where the modified discrete energy E [ φ n ] is defined by E [ φ k ] := E [ φ k ] + √ r k +1 τ k +1 κ (1 + r k +1 ) (cid:13)(cid:13) ∂ τ φ k (cid:13)(cid:13) − for k ≥ . (2.8) Proof.
The first condition of (2.7) ensures the unique solvability. We will establish the energydissipation law under the second condition of (2.7). The volume conversation implies (cid:79) τ φ n ∈ ˚ V h for n ≥
1. Then we make the inner product of (1.7) by ( − ∆ h ) − (cid:79) τ φ n /κ and obtain1 κ (cid:10) D φ n , (cid:79) τ φ n (cid:11) − − (cid:15) (cid:10) ∆ h φ n , (cid:79) τ φ n (cid:11) + (cid:10) F (cid:48) ( φ n ) , (cid:79) τ φ n (cid:11) = 0 . (2.9)With the help of the summation by parts and 2 a ( a − b ) = a − b + ( a − b ) , the second termat the left hand side of (2.9) reads (cid:15) (cid:10) ∇ h φ n , (cid:79) τ ∇ h φ n (cid:11) = (cid:15) (cid:13)(cid:13) ∇ h φ n (cid:13)(cid:13) − (cid:15) (cid:13)(cid:13) ∇ h φ n − (cid:13)(cid:13) + (cid:15) (cid:13)(cid:13) ∇ h (cid:79) τ φ n (cid:13)(cid:13) . It is easy to check the following identity4 (cid:0) a − a (cid:1) ( a − b ) = (cid:0) a − (cid:1) − (cid:0) b − (cid:1) − a − b ) + 2 a ( a − b ) + (cid:0) a − b (cid:1) . Then the third term in (2.9) can be bounded by (cid:10) F (cid:48) ( φ n ) , (cid:79) τ φ n (cid:11) ≥ (cid:10) F ( φ n ) , (cid:11) − (cid:10) F ( φ n − ) , (cid:11) − (cid:13)(cid:13) (cid:79) τ φ n (cid:13)(cid:13) . Recalling the generalized H¨older inequality (2.2), one gets12 (cid:13)(cid:13) (cid:79) τ φ n (cid:13)(cid:13) ≤ (cid:13)(cid:13) ∇ h (cid:79) τ φ n (cid:13)(cid:13)(cid:13)(cid:13) (cid:79) τ φ n (cid:13)(cid:13) − ≤ (cid:15) (cid:13)(cid:13) ∇ h (cid:79) τ φ n (cid:13)(cid:13) + 18 (cid:15) (cid:13)(cid:13) (cid:79) τ φ n (cid:13)(cid:13) − . Thus it follows from (2.9) that1 κ (cid:10) D φ n , (cid:79) τ φ n (cid:11) − − (cid:15) (cid:13)(cid:13) (cid:79) τ φ n (cid:13)(cid:13) − + E [ φ n ] ≤ E [ φ n − ] for n ≥
2. (2.10)The second condition of (2.7) gives R L ( r n , r n +1 ) ≥ κτ n / (4 (cid:15) ) . Taking w j = (cid:79) τ φ j in the firstinequality of Lemma 2.1, it is not difficult to get1 κ (cid:10) D φ n , (cid:79) τ φ n (cid:11) − ≥ √ r n +1 τ n +1 κ (1 + r n +1 ) (cid:13)(cid:13) ∂ τ φ n (cid:13)(cid:13) − − √ r n τ n κ (1 + r n ) (cid:13)(cid:13) ∂ τ φ n − (cid:13)(cid:13) − + 18 (cid:15) (cid:13)(cid:13) (cid:79) τ φ n (cid:13)(cid:13) − . Combining it with (2.10) yields E [ φ n ] ≤ E [ φ n − ] for n ≥
2. It completes the proof.9 emark 3.
It is seen that this time-step constraint (2.7) requires τ n = O ( (cid:15) /κ ) and is essen-tially determined by the interface width parameter (cid:15) . In practical simulations, this conditionis acceptable since τ n = O ( (cid:15) /κ ) is always necessary in the L norm stability or convergenceanalysis, cf. [7–9, 29]. Lemma 2.2.
Let S0 and A1 hold. If the step sizes τ n fulfill (2.7) , the solution of BDF2 implicittime-stepping scheme (1.7) is bounded in the sense that (cid:13)(cid:13) φ n (cid:13)(cid:13) + (cid:13)(cid:13) ∇ h φ n (cid:13)(cid:13) ≤ c := (cid:112) (cid:15) − c + (2 + (cid:15) ) | Ω h | for n ≥ ,where c is dependent on the domain Ω , the interface parameter (cid:15) and the starting value φ , butindependent of the time t n , the time-step sizes τ n and the time-step ratios r n .Proof. Under the assumption A1 , the definition (2.8) of E [ φ n ] gives E [ φ ] ≤ E [ φ ] + τ κ (cid:13)(cid:13) ∂ τ φ (cid:13)(cid:13) − ≤ c Thus the discrete energy dissipation law in Theorem 2.2 implies c ≥ E [ φ n ] ≥ E [ φ n ] . Remindingthe inequality (cid:107) φ n (cid:107) l ≥ (cid:15) ) (cid:107) φ n (cid:107) − (1 + (cid:15) ) | Ω h | , due to the simple fact ( a − − (cid:15) ) ≥ E [ φ n ] to get4 c ≥ (cid:15) (cid:13)(cid:13) ∇ h φ n (cid:13)(cid:13) + 4 (cid:10) F ( φ n ) , (cid:11) ≥ (cid:15) (cid:13)(cid:13) ∇ h φ n (cid:13)(cid:13) + 2 (cid:15) (cid:107) φ n (cid:107) − (cid:15) (2 + (cid:15) ) | Ω h | , and then (cid:0)(cid:13)(cid:13) φ n (cid:13)(cid:13) + (cid:13)(cid:13) ∇ h φ n (cid:13)(cid:13)(cid:1) ≤ (cid:107) φ n (cid:107) + 2 (cid:13)(cid:13) ∇ h φ n (cid:13)(cid:13) ≤ (cid:15) − c + (2 + (cid:15) ) | Ω h | for n ≥ Our error analysis is closely related to the convolution form (1.11), so we need some detailproperties and discrete convolution inequalities with respect to the DOC kernels θ ( n ) n − j . It is toemphasize that the positive constants m , m and m involved in this section are independent ofthe time t n , time-step sizes τ n and the step ratios r n . Actually, they would take different valuesfor different choices of step ratios r n , but are bounded with respect to the changes of step ratios,even when r n approaches the user limit r user . Following the proofs of [22, Lemma 2.2, Corollary 2.1 and Lemma 2.3], we can obtain somesimple properties of the DOC kernels.
Lemma 3.1. If S0 holds, the DOC kernels θ ( n ) n − j defined in (1.8) satisfy:(I) The discrete kernels θ ( n ) n − j are positive definite;(II) The discrete kernels θ ( n ) n − j are positive and θ ( n ) n − j = 1 b ( j )0 n (cid:89) i = j +1 r i r i for ≤ j ≤ n ; III) n (cid:88) j =2 θ ( n ) n − j ≤ τ n such that n (cid:88) k =2 k (cid:88) j =2 θ ( k ) k − j ≤ t n for n ≥ . We introduce the following two ( n − × ( n −
1) matrices B := b (2)0 b (3)1 b (3)0 . . . . . . b ( n )1 b ( n )0 and Θ := θ (2)0 θ (3)1 θ (3)0 ... ... . . . θ ( n ) n − θ ( n ) n − · · · θ ( n )0 , where the discrete kernels b ( n ) n − k and θ ( n ) n − k are defined by (1.6) and (1.8), respectively. It followsfrom the discrete orthogonal identity (1.9) thatΘ = B − . (3.1)If the step ratios condition S0 holds, Lemma 2.1 shows that the real symmetric matrix B := B + B T (3.2)is positive definite, that is, w T B w = 2 n (cid:88) k =2 w k k (cid:88) j =2 b ( k ) k − j w j ≥ n (cid:88) k =1 R L ( r k , r k +1 ) τ k ( w k ) , where the function R L ( z, s ) is defined by (2.5) and the vector w := ( w , w , · · · , w n ) T . Accordingto Lemma 3.1 (I), the following symmetric matrixΘ := Θ + Θ T = B − + ( B − ) T = ( B − ) T ( B + B T ) B − = ( B − ) T BB − (3.3)is also positive definite in the sense of w T Θ w = 2 (cid:80) n,kk,j θ ( k ) k − j w j w k >
0. Here and hereafter, wedenote (cid:80) n,kk,j := (cid:80) nk =2 (cid:80) kj =2 for the simplicity of presentation. To facilitate the proofs in what follows, we are to define the following step-scaled matrix (cid:101) B := Λ τ B Λ τ = ˜ b (2)0 ˜ b (3)1 ˜ b (3)0 . . . . . .˜ b ( n )1 ˜ b ( n )0 ( n − × ( n − , (3.4)where the the diagonal matrix Λ τ := diag (cid:0) √ τ , √ τ , · · · , √ τ n (cid:1) so that the step-scaled discretekernels ˜ b ( k )0 and ˜ b ( k )1 are given by˜ b ( k )0 = 1 + 2 r k r k and ˜ b ( k )1 = − r / k r k for 2 ≤ k ≤ n. (3.5)11oreover, we will use the following real symmetric matrix, (cid:101) B := (cid:101) B + (cid:101) B T = Λ τ B Λ τ . (3.6)The following two lemmas present some eigenvalue estimates of (cid:101) B and (cid:101) B T (cid:101) B . To avoidpossible confusions, we define the vector norm (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) · (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) by (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) := √ u T u for any real vector u andthe associated matrix norm (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) := (cid:113) λ max (cid:0) A T A (cid:1) . Lemma 3.2. If S0 holds, there exists a positive constant m such that λ min (cid:0) (cid:101) B (cid:1) ≥ m > .Proof. This proof can be followed from [18, Lemma A.1]. We include the main ingredient forthe completeness. Applying the Gerschgorin’s circle theorem to the matrix (cid:101) B , one has λ min (cid:0) (cid:101) B (cid:1) ≥ min ≤ k ≤ n R L ( r k , r k +1 ) > R L ( r user , r user ) = 2(1 + 2 r user − r / )1 + r user > , where R L (cid:0) z, s (cid:1) is defined by (2.5). It completes the proof by taking m = r user − r / )1+ r user . Lemma 3.3. If S0 holds, there exists a positive constant m such that λ max (cid:0) (cid:101) B T (cid:101) B (cid:1) ≤ m .Proof. This proof can be followed from [18, Lemma A.2]. We include the main ingredient forthe completeness. By writing out the tri-diagonal matrix (cid:101) B T (cid:101) B and applying the Gerschgorin’scircle theorem, one can find λ max (cid:0) (cid:101) B T (cid:101) B (cid:1) ≤ max ≤ k ≤ n R U ( r k , r k +1 ) < R U ( r user , r user ) , where the function R U ( z, s ) is defined by R U ( z, s ) := (1 + 2 z )(1 + 2 z + z / )(1 + z ) + s / (1 + 2 s + s / )(1 + s ) for 0 ≤ z, s < r user . An upper bound is then obtained by taking m = R U ( r user , r user ).By the above two lemmas, we can bound the minimum eigenvalue of Θ. Lemma 3.4. If S0 holds, the real symmetric matrix Θ in (3.3) satisfies v T Θ v ≥ m m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Λ τ v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) for any vector v . Proof.
Lemma 3.2 says that real symmetric matrix (cid:101) B ia positive definite. There exists a non-singular upper triangular matrix (cid:101) U such that (cid:101) B = (cid:101) U T (cid:101) U . By using (3.3) and (3.6), one gets v T Θ v = v T ( B − ) T BB − v = v T ( B − ) T Λ − τ (cid:101) B Λ − τ B − v = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:101) U Λ − τ B − v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . Thus it follows that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Λ τ v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Λ τ B Λ τ (cid:101) U − (cid:101) U Λ − τ B − v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:101) B (cid:101) U − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:101) U Λ − τ B − v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:101) B (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:101) U − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v T Θ v = λ max (cid:0) (cid:101) B T (cid:101) B (cid:1) λ max (cid:0) (cid:101) B − (cid:1) v T Θ v . Thus Lemmas 3.2 and 3.3 yield the claimed inequality.12o evaluate the maximum eigenvalue of Θ, consider the inverse matrix of the matrix (cid:101) B , (cid:101) Θ := (cid:101) B − = Λ − τ Θ Λ − τ = ˜ θ (2)0 ˜ θ (3)1 ˜ θ (3)0 ... ... . . .˜ θ ( n ) n − ˜ θ ( n ) n − · · · ˜ θ ( n )0 , (3.7)where the step-scaled DOC kernels ˜ θ ( k ) k − j follow from Lemma 3.1 (II),˜ θ ( k ) k − j := 1 √ τ k τ j θ ( k ) k − j = 1 + r j r j k (cid:89) i = j +1 r / i r i for 2 ≤ j ≤ k ≤ n. (3.8) Lemma 3.5. If S0 holds, then there exists a positive constant m such that v T Θ v ≤ m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Λ τ v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) for any vector v . Proof.
Let (cid:101)
Θ = (cid:101) Θ + (cid:101) Θ T . Since 0 < x / x < m ∗ := r / r user < x ∈ [0 , r user ], one can applythe formula (3.8) to get R n,k := k (cid:88) j =2 ˜ θ ( k ) k − j + n (cid:88) j = k ˜ θ ( j ) j − k ≤ k (cid:88) j =2 m k − j ∗ + n (cid:88) j = k m j − k ∗ < − m ∗ for 2 ≤ k ≤ n .One has λ max (cid:0) (cid:101) Θ (cid:1) ≤ max ≤ k ≤ n R n,k < m := − m ∗ by the Gerschgorin’s circle theorem. Itimplies w T (cid:101) Θ w ≤ m (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) w (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) for any w and the choice w := Λ τ v completes the proof. The following two lemmas describe the Young-type convolution inequality.
Lemma 3.6. If S0 holds, then for any real sequences { v k } nk =2 and { w k } nk =2 , n,k (cid:88) k,j θ ( k ) k − j w k v j ≤ ε n,k (cid:88) k,j θ ( k ) k − j v k v j + 12 m ε n (cid:88) k =2 τ k ( w k ) for ∀ ε > .Proof. Let w := ( w , w , · · · , w n ) T . A similar proof of [18, Lemma A.3] gives n,k (cid:88) k,j θ ( k ) k − j v j w k ≤ ε n,k (cid:88) k,j θ ( k ) k − j v j v k + 12 ε w T B − w for any ε > . From the proof Lemma 3.4, we have B − = Λ τ (cid:101) U − (cid:0) Λ τ (cid:101) U − (cid:1) T and then w T B − w = w T Λ τ (cid:101) U − (cid:0) Λ τ (cid:101) U − (cid:1) T w = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:0) (cid:101) U − (cid:1) T Λ τ w (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:0) (cid:101) U − (cid:1) T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Λ τ w (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = λ max (cid:0) ( (cid:101) B ) − (cid:1) w T Λ τ w ≤ m − n (cid:88) k =2 τ k ( w k ) , where Lemma 3.2 has been used. It completes the proof.13 emma 3.7. If S0 holds, then for any real sequences { v k } nk =2 and { w k } nk =2 , n,k (cid:88) k,j θ ( k ) k − j w k v j ≤ ε n (cid:88) k =2 τ k ( v k ) + m m ε n (cid:88) k =2 τ k ( w k ) for ∀ ε > . Proof.
For fixed time index n , taking ε := 2 ε / m in Lemma 3.6 yields n,k (cid:88) k,j θ ( k ) k − j w k v j ≤ ε m n,k (cid:88) k,j θ ( k ) k − j v k v j + m m ε n (cid:88) k =2 τ k ( w k ) ≤ ε n (cid:88) k =2 τ k ( v k ) + m m ε n (cid:88) k =2 τ k ( w k ) , where Lemma 3.5 was used in the last inequality. It completes the proof by choosing ε := ε .We now present two discrete embedding-type convolution inequalities by considering threetime-space discrete functions u k , v k and w k (1 ≤ k ≤ n ) in the space V h or its subspace ˚ V h . Lemma 3.8.
Assume that u k , w k ∈ V h , v k ∈ ˚ V h (2 ≤ k ≤ n ) and there exists a constant c u such that (cid:13)(cid:13) u k (cid:13)(cid:13) l ≤ c u for ≤ k ≤ n . If S0 holds, then for any ε > , n,k (cid:88) k,j θ ( k ) k − j (cid:10) u j v j , w k (cid:11) ≤ ε n,k (cid:88) k,j θ ( k ) k − j (cid:10) ∇ h v j , ∇ h v k (cid:11) + c z c u m m m ε n (cid:88) k =2 τ k (cid:13)(cid:13) w k (cid:13)(cid:13) . Proof.
For fixed time index n , taking v j := u jh v jh and ε := ε in Lemma 3.7, we have n,k (cid:88) k,j θ ( k ) k − j (cid:10) u j v j , w k (cid:11) ≤ ε n (cid:88) k =2 τ k (cid:13)(cid:13) u k v k (cid:13)(cid:13) + m m ε n (cid:88) k =2 τ k (cid:13)(cid:13) w k (cid:13)(cid:13) . The well–known H¨older inequality and the discrete embedding inequality (2.3) imply that (cid:13)(cid:13) u k v k (cid:13)(cid:13) ≤ (cid:13)(cid:13) u k (cid:13)(cid:13) l (cid:13)(cid:13) v k (cid:13)(cid:13) l ≤ c z (cid:13)(cid:13) u k (cid:13)(cid:13) l (cid:13)(cid:13) ∇ h v k (cid:13)(cid:13) ≤ c z c u (cid:13)(cid:13) ∇ h v k (cid:13)(cid:13) . We derive that n (cid:88) k =2 τ k (cid:13)(cid:13) u k v k (cid:13)(cid:13) ≤ c z c u n (cid:88) k =2 τ k (cid:13)(cid:13) ∇ h v k (cid:13)(cid:13) . Then it follows that n,k (cid:88) k,j θ ( k ) k − j (cid:10) u j v j , w k (cid:11) ≤ ε c z c u n (cid:88) k =2 τ k (cid:13)(cid:13) ∇ h v k (cid:13)(cid:13) + m m ε n (cid:88) k =2 τ k (cid:13)(cid:13) w k (cid:13)(cid:13) . (3.9)Following the proof of Lemma 3.4, it is not difficult to get (cf. [19]) n (cid:88) k =2 τ k (cid:13)(cid:13) ∇ h v k (cid:13)(cid:13) ≤ m m n,k (cid:88) k,j θ ( k ) k − j (cid:10) ∇ h v j , ∇ h v k (cid:11) . Inserting this inequality into (3.9) and choosing the parameter ε := m ε/ (2 c z c u m ) , we get theclaimed inequality and complete the proof. 14 emma 3.9. Assume that u k ∈ V h , w k ∈ ˚ V h (2 ≤ k ≤ n ) and there exists a constant c u suchthat (cid:13)(cid:13) u k (cid:13)(cid:13) l ≤ c u for ≤ k ≤ n . If S0 holds, then for any ε > , n,k (cid:88) k,j θ ( k ) k − j (cid:10) u j w j , ∆ h w k (cid:11) ≤ ε n,k (cid:88) k,j θ ( k ) k − j (cid:10) ∆ h w j , ∆ h w k (cid:11) + c z c u m m m ε n (cid:88) k =2 τ k (cid:13)(cid:13) w k (cid:13)(cid:13) . Proof.
For fixed time index n , we start the proof from (3.9) by setting w j := ∆ h w j , v j := w j and ε := m m / ( ε m ), that is, n,k (cid:88) k,j θ ( k ) k − j (cid:10) u j w j , ∆ h w k (cid:11) ≤ c z c u m m m ε n (cid:88) k =2 τ k (cid:13)(cid:13) ∇ h w k (cid:13)(cid:13) + m ε m n (cid:88) k =2 τ k (cid:13)(cid:13) ∆ h w k (cid:13)(cid:13) ≤ c z c u m m m ε n (cid:88) k =2 τ k (cid:13)(cid:13) ∇ h w k (cid:13)(cid:13) + ε n,k (cid:88) k,j θ ( k ) k − j (cid:10) ∆ h w j , ∆ h w k (cid:11) , (3.10)where Lemma 3.4 has been used to handle the last term. Furthermore, by using the classicalYoung’s inequality and Lemma 3.4, one gets n (cid:88) k =2 τ k (cid:13)(cid:13) ∇ h w k (cid:13)(cid:13) = n (cid:88) k =2 τ k (cid:10) − ∆ h w k , w k (cid:11) ≤ ε n (cid:88) k =2 τ k (cid:13)(cid:13) ∆ h w k (cid:13)(cid:13) + 12 ε n (cid:88) k =2 τ k (cid:13)(cid:13) w k (cid:13)(cid:13) ≤ m ε m n,k (cid:88) k,j θ ( k ) k − j (cid:10) ∆ h w j , ∆ h w k (cid:11) + 12 ε n (cid:88) k =2 τ k (cid:13)(cid:13) w k (cid:13)(cid:13) . Inserting this inequality into (3.10), we have n,k (cid:88) k,j θ ( k ) k − j (cid:10) u j w j , ∆ h w k (cid:11) ≤ (cid:16) c z c u m m ε m ε + ε (cid:17) n,k (cid:88) k,j θ ( k ) k − j (cid:10) ∆ h w j , ∆ h w k (cid:11) + c z c u m m m ε ε n (cid:88) k =2 τ k (cid:13)(cid:13) w k (cid:13)(cid:13) . Now by choosing ε := ε/ ε := m ε ε/ ( c z c u m m ), we obtain the claimed inequality. L norm error estimate Let ξ j Φ := D Φ( t j ) − ∂ t Φ( t j ) be the local consistency error of BDF2 formula at the time t = t j .We will consider a convolutional consistency error Ξ k Φ defined byΞ k Φ := k (cid:88) j =2 θ ( k ) k − j ξ j Φ = k (cid:88) j =2 θ ( k ) k − j ( D Φ( t j ) − ∂ t Φ( t j )) for k ≥
2. (4.1)By following the proof of [18, Lemma 3.4], one has (cid:12)(cid:12) Ξ k Φ (cid:12)(cid:12) ≤ k (cid:88) j =1 θ ( k ) k − j τ j (cid:90) t j t j − (cid:12)(cid:12) Φ (cid:48)(cid:48)(cid:48) ( s ) (cid:12)(cid:12) d s for k ≥ emma 4.1. If S0 holds, the convolutional consistency error Ξ k Φ in (4.1) satisfies n (cid:88) k =2 (cid:12)(cid:12) Ξ k Φ (cid:12)(cid:12) ≤ t n max ≤ j ≤ n (cid:16) τ j (cid:90) t j t j − (cid:12)(cid:12) Φ (cid:48)(cid:48)(cid:48) ( s ) (cid:12)(cid:12) d s (cid:17) for n ≥ . We use the standard seminorms and norms in the Sobolev space H m (Ω) for m ≥
0. Let C ∞ per (Ω) be a set of infinitely differentiable L -periodic functions defined on Ω, and H mper (Ω) bethe closure of C ∞ per (Ω) in H m (Ω), endowed with the semi-norm | · | H mper and the norm (cid:107)·(cid:107) H mper .For simplicity, denote | · | H m := | · | H mper , (cid:107)·(cid:107) H m := (cid:107)·(cid:107) H mper , and (cid:107)·(cid:107) L := (cid:107)·(cid:107) H . Next lemmalists some approximations, cf. [25, 26], of the L -projection operator P M and trigonometricinterpolation operator I M defined in subsection 2.1. Lemma 4.2.
For any u ∈ H qper (Ω) and ≤ s ≤ q , it holds that (cid:107) P M u − u (cid:107) H s ≤ C u h q − s | u | H q , (cid:107) P M u (cid:107) H s ≤ C u (cid:107) u (cid:107) H s ; (4.2) and, in addition if q > / , (cid:107) I M u − u (cid:107) H s ≤ C u h q − s | u | H q , (cid:107) I M u (cid:107) H s ≤ C u (cid:107) u (cid:107) H s . (4.3) Note that, the energy dissipation law (1.4) of CH model (1.2) shows that E [Φ n ] ≤ E [Φ( t )]. Fromthe formulation (1.1), it is easy to check that (cid:13)(cid:13) Φ n (cid:13)(cid:13) H can be bounded by a time-independentconstant. Let Φ nM := (cid:0) P M Φ (cid:1) ( · , t n ) be the L -projection of exact solution at time t = t n . Theprojection estimate (4.2) in Lemma 4.2 yields (cid:13)(cid:13) Φ nM (cid:13)(cid:13) + (cid:13)(cid:13) ∇ h Φ nM (cid:13)(cid:13) ≤ (cid:13)(cid:13) P M Φ n (cid:13)(cid:13) H ≤ c for 1 ≤ n ≤ N , (4.4)where c is dependent on the domain Ω and initial data Φ( t ), but independent of the time t n .We are in the position to prove the L norm convergence of the adaptive BDF2 scheme (1.7). Theorem 4.1.
Assume that the CH problem (1.2) has a solution Φ ∈ C (cid:0) [0 , T ]; H m +4 per (cid:1) forsome integer m ≥ . Suppose further that the step-ratios condition S0 and the time-step sizeconstraint (2.7) hold such that the BDF2 scheme (1.7) is unique solvable and energy stable. Ifthe maximum step size τ ≤ /c (cid:15) , the solution φ n is robustly convergent in the L norm, (cid:13)(cid:13) Φ n − φ n (cid:13)(cid:13) ≤ C φ exp (cid:0) c (cid:15) t n − (cid:1)(cid:16) (cid:13)(cid:13) Φ M − φ (cid:13)(cid:13) + τ (cid:13)(cid:13) ∂ τ (Φ M − φ ) (cid:13)(cid:13) + t n h m + t n τ max A substitution of the pro-jection solution Φ M and differentiation operator ∆ h into the original equation (1.2) yields thesemi-discrete system ∂ t Φ M = κ ∆ h µ M + ζ P with µ M = F (cid:48) (Φ M ) − (cid:15) ∆ h Φ M , (4.6)where ζ P ( x h , t ) represents the spatial consistency error arising from the projection of exactsolution, that is, ζ P := ∂ t Φ M − ∂ t Φ + κ (∆ µ − ∆ h µ M ) for x h ∈ Ω h . (4.7)Following the proof of [18, Theorem 3.1], and using Lemma 4.2, it is not difficult to obtainthat (cid:107) ζ P (cid:107) ≤ C φ h m and (cid:107) ζ P ( t j ) (cid:107) ≤ C φ h m for j ≥ 2. Then Lemma 3.1 (III) yields n (cid:88) k =2 (cid:13)(cid:13) Υ kP (cid:13)(cid:13) ≤ C φ h m n (cid:88) k =2 k (cid:88) j =2 θ ( k ) k − j ≤ C φ t n h m where Υ kP := k (cid:88) j =2 θ ( k ) k − j ζ P ( t j ) for k ≥ 2. (4.8) Stage 2: L norm error of fully discrete system From the projection equation (4.6), onecan apply the BDF2 formula to obtain the following approximation equation D Φ nM = κ ∆ h µ nM + ζ nP + ξ n Φ with µ nM = F (cid:48) (Φ nM ) − (cid:15) ∆ h Φ nM , (4.9)where ξ n Φ is the local consistency error of BDF2 formula and ζ nP := ζ P ( t n ) is defined by (4.7).Subtracting the full discrete scheme (1.7) from the approximation equation (4.9), we have thefollowing error system D e n = κ ∆ h (cid:0) − (cid:15) ∆ h e n + f nφ e n (cid:1) + ζ nP + ξ n Φ for 2 ≤ n ≤ N , (4.10)17here the nonlinear term f nφ := (Φ nM ) + Φ nM φ n + ( φ n ) − 1. Thanks to the solution estimatesin Lemma 2.2 and (4.4), one applies the discrete embedding inequality (2.1) to find that (cid:13)(cid:13) f nφ (cid:13)(cid:13) l ≤ (cid:13)(cid:13) Φ nM (cid:13)(cid:13) l + (cid:13)(cid:13) Φ nM (cid:13)(cid:13) l (cid:13)(cid:13) φ n (cid:13)(cid:13) l + (cid:13)(cid:13) φ n (cid:13)(cid:13) l + (cid:12)(cid:12) Ω h (cid:12)(cid:12) / ≤ c ( c + c c + c ) + (cid:12)(cid:12) Ω h (cid:12)(cid:12) / = c . (4.11)Multiplying both sides of equation (4.10) by the DOC kernels θ ( k ) k − n , and summing up n from n = 2 to k , we apply the equality (1.10) with v j = e j to obtain (cid:79) τ e k = − θ ( k ) k − b (2)1 (cid:79) τ e + κ k (cid:88) j =2 θ ( k ) k − j ∆ h (cid:0) − (cid:15) ∆ h e j + f jφ e j (cid:1) + Υ kP + Ξ k Φ (4.12)for 2 ≤ k ≤ N , where Ξ k Φ and Υ kP are defined by (4.1) and (4.8), respectively. Making the innerproduct of (4.12) with 2 e k , and summing k from 2 to n , we obtain (cid:13)(cid:13) e n (cid:13)(cid:13) ≤ (cid:13)(cid:13) e (cid:13)(cid:13) − n (cid:88) k =2 θ ( k ) k − b (2)1 (cid:13)(cid:13) e k (cid:13)(cid:13)(cid:13)(cid:13) (cid:79) τ e (cid:13)(cid:13) + J n + 2 n (cid:88) k =2 (cid:10) Υ kP + Ξ k Φ , e k (cid:11) (4.13)for 2 ≤ n ≤ N , where J n is defined by J n := 2 κ n,k (cid:88) k,j θ ( k ) k − j (cid:10) − (cid:15) ∆ h e j + f jφ e j , ∆ h e k (cid:11) . (4.14)Taking u j := f jφ (with the upper bound c u := c ), w j := e j and ε = (cid:15) in Lemma 3.9, one appliesthe solution bound (4.11) to obtain2 κ n,k (cid:88) k,j θ ( k ) k − j (cid:10) f jφ e j , ∆ h e k (cid:11) ≤ κ(cid:15) n,k (cid:88) k,j θ ( k ) k − j (cid:10) ∆ h e j , ∆ h e k (cid:11) + 2 κc z c m m m (cid:15) n (cid:88) k =2 τ k (cid:13)(cid:13) e k (cid:13)(cid:13) . Then the term J n in (4.14) can be bounded by J n ≤ c (cid:15) n (cid:88) k =2 τ k (cid:13)(cid:13) e k (cid:13)(cid:13) . Therefore, it follows from (4.13) that (cid:13)(cid:13) e n (cid:13)(cid:13) ≤ (cid:13)(cid:13) e (cid:13)(cid:13) − n (cid:88) k =2 θ ( k ) k − b (2)1 (cid:13)(cid:13) e k (cid:13)(cid:13)(cid:13)(cid:13) (cid:79) τ e (cid:13)(cid:13) + c (cid:15) n (cid:88) k =2 τ k (cid:13)(cid:13) e k (cid:13)(cid:13) + 2 n (cid:88) k =2 (cid:13)(cid:13) e k (cid:13)(cid:13)(cid:13)(cid:13) Υ kP + Ξ k Φ (cid:13)(cid:13) for 2 ≤ n ≤ N . Choosing some integer n (1 ≤ n ≤ n ) such that (cid:13)(cid:13) e n (cid:13)(cid:13) = max ≤ k ≤ n (cid:13)(cid:13) e k (cid:13)(cid:13) .Taking n := n in the above inequality, one can obtain (cid:13)(cid:13) e n (cid:13)(cid:13) ≤ (cid:13)(cid:13) e (cid:13)(cid:13) − (cid:13)(cid:13) ∂ τ e (cid:13)(cid:13) n (cid:88) k =2 θ ( k ) k − b (2)1 τ + c (cid:15) n (cid:88) k =2 τ k (cid:13)(cid:13) e k (cid:13)(cid:13) + 2 n (cid:88) k =2 (cid:13)(cid:13) Υ kP + Ξ k Φ (cid:13)(cid:13) . 18y using Lemma 3.1(II), one has − θ ( k ) k − b (2)1 τ = τ k (cid:89) i =2 r i r i = τ k k (cid:89) i =2 r i r i ≤ τ k k − for 2 ≤ k ≤ N , such that − n (cid:88) k =2 θ ( k ) k − b (2)1 τ ≤ τ n (cid:88) k =2 k − ≤ τ for 2 ≤ n ≤ N . Thus one gets (cid:13)(cid:13) e n (cid:13)(cid:13) ≤ (cid:13)(cid:13) e n (cid:13)(cid:13) ≤ (cid:13)(cid:13) e (cid:13)(cid:13) + 2 τ (cid:13)(cid:13) ∂ τ e (cid:13)(cid:13) + c (cid:15) n (cid:88) k =2 τ k (cid:13)(cid:13) e k (cid:13)(cid:13) + 2 n (cid:88) k =2 (cid:13)(cid:13) Υ kP + Ξ k Φ (cid:13)(cid:13) . Under the maximum step constraint τ ≤ /c (cid:15) , we have (cid:13)(cid:13) e n (cid:13)(cid:13) ≤ (cid:13)(cid:13) e (cid:13)(cid:13) + 4 τ (cid:13)(cid:13) ∂ τ e (cid:13)(cid:13) + c (cid:15) n − (cid:88) k =2 τ k (cid:13)(cid:13) e k (cid:13)(cid:13) + 4 n (cid:88) k =2 (cid:13)(cid:13) Υ kP + Ξ k Φ (cid:13)(cid:13) . The discrete Gr¨onwall inequality [22, Lemma 3.1] yields the following estimate (cid:13)(cid:13) e n (cid:13)(cid:13) ≤ (cid:0) c (cid:15) t n − (cid:1)(cid:16)(cid:13)(cid:13) e (cid:13)(cid:13) + 2 τ (cid:13)(cid:13) ∂ τ e (cid:13)(cid:13) + 2 n (cid:88) k =2 (cid:13)(cid:13) Υ kP (cid:13)(cid:13) + 2 n (cid:88) k =2 (cid:13)(cid:13) Ξ k Φ (cid:13)(cid:13)(cid:17) for 2 ≤ n ≤ N .Furthermore, Lemma 4.1 together with (cid:13)(cid:13) ∂ t Φ (cid:13)(cid:13) = (cid:13)(cid:13) I M ∂ t Φ (cid:13)(cid:13) L ≤ C φ (cid:13)(cid:13) ∂ t Φ (cid:13)(cid:13) L , due to Lemma 4.2,gives the bound of the global temporal error term (cid:80) nk =2 (cid:13)(cid:13) Ξ k Φ (cid:13)(cid:13) . Therefore by applying the errorestimate (4.8) and the triangle inequality (4.5), we complete the proof. We run the variable-steps BDF2 scheme (1.7) for the CH equation (1.2). The TR-BDF2 methodis always employed to obtain the first-level solution. A simple fixed-point iteration with thetermination error 10 − is employed to solve the nonlinear algebra equations at each time level. Example 1. To facilitate the robustness test of the variable-steps BDF2 method (1.7) , we con-sider an exact solution Φ( x , t ) = cos( t ) sin( x ) sin( y ) with the model parameters κ = 1 and (cid:15) = 0 . by adding a corresponding exterior force on the right hand side of the CH model (1.2) . In the following examinations, the computational domain (0 , π ) is discretized by using 128 spatial meshes. Then the problem is solved until time T = 1 on random time meshes. To be moreprecise, we take the time step sizes τ k := T σ k /S for 1 ≤ k ≤ N , where σ k ∈ (0 , 1) is the uniformlydistributed random number and S = (cid:80) Nk =1 σ k . Since the spectral accuracy in space is standard,we only test the time accuracy with the numerical error e ( N ) := max ≤ n ≤ N (cid:107) Φ( t n ) − φ n (cid:107) ineach run. The numerical order of convergence is estimated byOrder := log ( e ( N ) /e (2 N ))log ( τ ( N ) /τ (2 N )) , N τ e ( N ) Order max r k N 40 3.96e-02 1.82e-04 1.69 17.27 380 2.44e-02 4.60e-05 2.84 46.22 5160 1.29e-02 1.16e-05 2.16 167.41 16320 6.28e-03 2.90e-06 1.93 264.04 29640 3.05e-03 7.27e-07 1.92 1584.01 62where τ ( N ) denotes the maximum time-step size for total N subintervals.The numerical results obtained using a set of random meshes are tabulated in Table 1. Inaddition to the discrete L numerical error between the exact solution and the numerical solution,the maximum time-step size τ , the maximum step ratio max r k and the number (denote by N )of time levels with the step ratios r k ≥ . 864 are also recorded, respectively.As observed, the variable-steps BDF2 method (1.7) still achieves the second-order accuracyon arbitrary nonuniform meshes even though some step ratios lager than the update zero-stabilitylimit r ∗ ≈ . < r k < . 864 is still a sufficient condition for second-order convergence.In the following experiments, we always set r user = 4, which is large enough for practical adaptivetime-stepping simulations. (a) time-step size τ = 10 − (b) time-step size τ = 2 × − (c) time-step size τ = 10 − Figure 1: Solution curves by BDF2, CN and CNCS methods at T = 0 . Example 2. We next simulate the coarsening dynamics of the CH equation (1.2) . Precisely, theinitial condition is taken as Φ ( x ) = rand ( x ) , where rand ( x ) generates random numbers between − . to . uniformly. Here, the mobility coefficient κ = 0 . and the interfacial thickness (cid:15) = 0 . are taken in the following numerical simulations. Always, the spatial domain (0 , π ) is discretized by using spatial meshes. To further benchmark the BDF2 scheme with the random initial data generated in Example 2,we run several numerical tests to explore the numerical behaviors near the initial time. We also20 a) time-step size τ = 10 − (b) time-step size τ = 2 × − (c) time-step size τ = 10 − Figure 2: Original energy curves by BDF2, CN and CNCS methods until T = 0 . ∂ τ φ n = κ ∆ h µ n − with µ n − = 12 (cid:2) ( φ n ) + ( φ n − ) (cid:3) φ n − − φ n − − ε ∆ h φ n − , and the second-order Crank-Nicolson convex-splitting (denoted by CNCS) method [9], ∂ τ φ n = κ ∆ h ˆ µ n − with ˆ µ n − = 12 (cid:2) ( φ n ) + ( φ n − ) (cid:3) φ n − − ˇ φ n − − ε ∆ h ˆ φ n − , where φ n − := ( φ n + φ n − ) / 2, ˆ φ n − := (cid:0) φ n + φ n − (cid:1) / φ n − = (cid:0) φ n − − φ n − (cid:1) / 2. Sincethe CNCS method requires two initialization steps, a first-order convex-splitting scheme [8] isused here to obtain the first-level solution. (a) Original energy (b) Time step sizes Figure 3: Energy curves and adaptive time-step sizes for different parameters β .The random initial data initiates a fast coarsening dynamics at the beginning time. Weuse a random initial profile to test the effectiveness of various numerical methods with differenttime step sizes. The numerical solution curves are summarized in Figure 1, where the referencesolution is obtained by using the BDF2 method with a uniform time-step size τ = 10 − . Weobserve that solutions of CN and CNCS methods tend to generate non-physical oscillations whensome large time steps are used. In contrast, the BDF2 solution is more robust and accuratethan the CN and CNCS schemes with the same time step size. It seems that the BDF2 methodis more suitable than Crank-Nicolson type schemes when large time-step sizes are adopted.21able 2: CPU time (in seconds) and total time steps comparisons.Strategies τ = 10 − β = 10 β = 10 β = 10 CPU time 337.736 13.870 26.568 62.797Time levels 30000 493 1190 3554 In this subsection, we simulate the coarsening dynamics by using the variable-steps BDF2method (1.7) with the random initial condition. In what follows, to capture the multiple timescales accurately and to improve the computational efficiency for long-time simulations, the timesteps are selected by using the following adaptive time-stepping strategy [16], τ ada = max (cid:40) τ min , τ max (cid:113) β (cid:13)(cid:13) ∂ τ φ n (cid:13)(cid:13) (cid:41) so that τ n +1 = min (cid:8) τ ada , r user τ n (cid:9) , (5.1)hwere (cid:13)(cid:13) ∂ τ φ n (cid:13)(cid:13) denotes the change rate between two consecutive time step numerical solutions, β > τ max and τ min are the predetermined maximum and minimumtime steps, respectively. (a) time t = 10 (b) time t = 50 (c) time t = 100(d) time t = 200 (e) time t = 300 (f) time t = 500 Figure 4: The profile of numerical solution φ at different time for the CH model.We take τ min = 10 − and τ max = 10 − in the adaptive time-stepping algorithm (5.1), and runthe BDF2 method (1.7) until time T = 30. The reference solution is obtained by applying a small22 a) Energy scaling (b) Volume difference (c) Adaptive step sizes Figure 5: Numerical results show original energy, volume and adaptive time steps of the CHequation during the coarsening dynamics.time step τ = 10 − . As seen in Figure 3, we use three different user parameters β = 10 , and10 to compute the discrete original energy and the corresponding adaptive time-steps. One canobserve that the discrete energy curves using the adaptive stepping algorithm are comparableto the the reference one. On the other hand, the adjustments of time-steps are closely reliedon the user parameter β . As expected, a large β leads to small time-step sizes, and a small β generates large step sizes. The CPU time (in seconds) and the adaptive time levels recordedin Table 2 show the effectiveness and efficiency of the adaptive time-stepping algorithm, whichmakes the long-time dynamics simulations practical.We next perform the coarsening dynamic simulations by using the above adaptive time-stepping strategy with the setting β = 10 until time T = 500. The evolution of microstructurefor the CH model due to the phase separation at different time are summarized in Figure 4. Asseen, the microstructure is relatively fine and consists of many precipitations at early time. Thecoarsening, dissolution, merging processes are also observed. The time evolutions of originalenergy, volume and the adaptive step sizes are summarized in Figure 5. The subplot (a) ofFigure 5 demonstrates a very good agreement with the expected scaling law, i.e., the energydecrease as O ( t − ). References [1] R. Alexander , Diagonally implicit Runge-Kutta methods for stiff ODEs , SIAM J. Numer.Anal., 14(6) (1977), pp. 1006–1021.[2] J. Becker , A second order backward difference method with variable steps for a parabolicproblem , BIT, 38(4) (1998), pp. 644–662.[3] A. Bertozzi, S. Esedoglu, and A. Gillette , Inpainting of binary images using theCahn-Hilliard equation , IEEE Trans. Image Process., 16 (2007), pp. 285–291.[4] J. Cahn and J. Hilliard , Free energy of a nonuniform system I. interfacial free energy ,J. Chem. Phys., 28 (1958), pp. 258–267.[5] K. Cheng, W. Feng, C. Wang, and S. Wise , An energy stable fourth order finitedifference scheme for the Cahn-Hilliard equation , J. Comput. Appl. Math., 362 (2019), pp.574–595. 236] K. Cheng, C. Wang, and S. Wise , An energy stable BDF2 Fourier pseudo-spectralnumerical scheme for the square phase field crystal equation , Comm. Comput. Phys., 26(2019), pp. 1335–1364.[7] W. Chen, C. Wang and X. Wang , A linear iteration algorithm for a second-order energystable scheme for a thin film model without slope selection , J Sci. Comput., 59 (2014), pp.574–601.[8] W. Chen, X. Wang, Y. Yan and Z. Zhang , A second order BDF numerical schemewith variable steps for the Cahn–Hilliard equation , SIAM J. Numer. Anal., 57 (1) (2019),pp. 495–525.[9] K. Cheng, C. Wang, S. Wise, and X. Yue , A second-order, weakly energy-stable pseudo-spectral scheme for the Cahn-Hilliard equation and its solution by the homogeneous lineariteration method , J. Sci. Comput., 69 (2016), pp. 1083–1114.[10] V. Cristini, X. Li, J. Lowengrub, and S. Wise. Nonlinear simulations of solid tumor growthusing a mixture model: invasion and branching. J. Math. Biol. , 58 (2009), pp. 723–763.[11] E. Emmrich , Stability and error of the variable two-step BDF for semilinear parabolicproblems , J. Appl. Math. & Computing, 19 (2005), pp. 33–55.[12] R.D. Grigorieff , Stability of multistep-methods on variable grids , Numer. Math., 42(1983), pp. 359–377.[13] H. Gomez and T. Hughes , Provably unconditionally stable, second-order time-accurate,mixed variational methods for phase-field models , J. Comput. Phys., 230 (2011), pp. 5310–5327.[14] E. Hairer, S.P. Nørsett and G. Wanner , Solving Ordinary Differential EquationsI: Nonstiff Problems , Volume 8 of Springer Series in Computational Mathematics, SecondEdition, Springer-Verlag, 1992.[15] M.E. Hosea and L.F. Shampine , Analysis and implementation of TR-BDF2 , Appl. Nu-mer. Math., 20 (1996), pp. 21-37.[16] J. Huang, C. Yang, and Y. Wei , Parallel energy-stable solver for a coupled Allen–Cahnand Cahn–Hilliard system , SIAM J. Sci. Comput., 42(5) (2020), pp. C294–C312.[17] M.-N. Le Roux , Variable step size multistep methods for parabolic problems , SIAM J.Numer. Anal., 19 (4) (1982), pp. 725–741.[18] H.-L. Liao, B. Ji and L. Zhang , An adaptive BDF2 implicit time-stepping method forthe phase field crystal model , IMA J. Numer. Anal., 2020, doi:10.1093/imanum/draa075.[19] H.-L. Liao, X. Song, T. Tang and T. Zhou , Analysis of the second order BDF schemewith variable steps for the molecular beam epitaxial model without slope selection , Sci. ChinaMath., 2020, doi:10.1007/s11425-020-1817-4.[20] H.-L. Liao, T. Tang and T. Zhou , On energy stable, maximum-principle preserving,second order BDF scheme with variable steps for the Allen-Cahn equation , SIAM J. Numer.Anal., 58(4) (2020), pp. 2294-2314. 2421] H.-L. Liao, T. Tang and T. Zhou , Positive definiteness of real quadratic forms result-ing from variable-step approximations of convolution operators , arXiv:2011.13383v1, 2020,submitted.[22] H.-L. Liao and Z. Zhang , Analysis of adaptive BDF2 scheme for diffusion equations ,Math. Comp., 2020, DOI: 10.1090/mcom/3585.[23] H. Nishikawa , On large start-up error of BDF2 , J. Comput. Phys., 392 (2019), pp. 456–461.[24] Z. Qiao, Z. Zhang and T. Tang , An adaptive time-stepping strategy for the molecularbeam epitaxy models , SIAM J. Sci. Comput., 33 (2011), pp. 1395–1414.[25] J. Shen, C. Wang, X. Wang and S.M. Wise , Second-order convex splitting schemes forgradient flows with Ehrlich–Schwoebel type energy: application to thin film epitaxy , SIAMJ. Numer. Anal., 50(1) (2012), pp. 105–125.[26] J. Shen, T. Tang, and L. Wang , Spectral methods: Algorithms, analysis and applica-tions , Springer-Verlag, Berlin Heidelberg, 2011.[27] G. Tumolo and L. Bonaventura , A semi-implicit, semi- Lagrangian, DG frameworkfor adaptive numerical weather prediction , Quarterly Journal of the Royal MeteorologicalSociety, DOI: 10.1002/qj.2544, 2015.[28] W. Wang, Y. Chen and H. Fang On the variable two-step IMEX BDF method forparabolic integro-differential equations with nonsmooth initial data arising in finance , SIAMJ. Numer. Anal., 57(3) (2019), pp. 1289–1317.[29] J. Xu, Y.K. Li, S.N. Wu and A.Bousequet , On the stability and accuracy of partiallyand fully implicit schemes for phase field modeling , Comput. Methods Appl. Mech. Engrg.,345 (2019), pp. 826-853.[30]