An efficient two-grid fourth-order compact difference scheme with variable-step BDF2 method for the semilinear parabolic equation

Due to the lack of corresponding analysis on appropriate mapping operator between two grids, high-order two-grid difference algorithms are rarely studied. In this paper, we firstly discuss the boundedness of a local bi-cubic Lagrange interpolation operator. And then, taking the semilinear parabolic equation as an example, we first construct a variable-step high-order nonlinear difference algorithm using compact difference technique in space and the second-order backward differentiation formula (BDF2) with variable temporal stepsize in time. With the help of discrete orthogonal convolution (DOC) kernels and a cut-off numerical technique, the unique solvability and corresponding error estimates of the high-order nonlinear difference scheme are established under assumptions that the temporal stepsize ratio satisfies rk<4.8645 and the maximum temporal stepsize satisfies tau = o(h^1/2 ). Then, an efficient two-grid high-order difference algorithm is developed by combining a small-scale variable-step high-order nonlinear difference algorithm on the coarse grid and a large-scale variable-step high-order linearized difference algorithm on the fine grid, in which the constructed piecewise bi-cubic Lagrange interpolation mapping operator is adopted to project the coarse-grid solution to the fine grid. Under the same temporal stepsize ratio restriction rk<4.8645 and a weaker maximum temporal stepsize condition tau = o(H^1.2 ), optimal fourth-order in space and second-order in time error estimates of the two-grid difference scheme is established if the coarse-fine grid stepsizes satisfy H = O(h^4/7). Finally, several numerical experiments are carried out to demonstrate the effectiveness and efficiency of the proposed scheme.


Introduction
It is well known that the analytical solutions of nonlinear parabolic PDEs arising from a variety of physical and engineering applications are not available in most cases.Thus, numerous efforts have been devoted to the development of efficient numerical schemes, see [2,8,12,16,19,33].Generally speaking, fully-implicit numerical schemes are usually proved to be unconditionally stable.Unfortunately, at each time step, one has to solve a system of nonlinear equations [16,18], in which an extra iterative process must be imposed, and this in turn may cause severe computational costs.Instead, a very popular and widely-used approach is the so-called implicit-explicit scheme, which treats the linear term implicitly and the nonlinear term explicitly.However, if the corresponding globally continuous condition of the nonlinear term (e.g., | f ′′ (w)| ≤ K f for w ∈ R) cannot be imposed or the boundedness of numerical solution in L ∞ norm cannot be obtained, this method usually suffers from a very restrictive temporal stepsize condition caused by use of inverse inequality for convergence, e.g., τ = O(h d/2p ), where h is the spatial mesh size, d is the space dimension and p is the accuracy of the time discretization.Therefore, such restrictions may lead to the use of a small temporal stepsize, and thus much computational time may be consumed.Here we refer to [2,3,41,48] for an incomplete list of references.
Another efficient and powerful strategy is the two-grid method which is proposed by Xu et al. [33,45].The basic idea of this kind of method is to reduce the solution of the large-scale nonlinear problem on the fine grid to a small-scale nonlinear problem on the coarse grid and a large-scale linear problem on the fine grid.Hence, basically it includes two solution steps: First, one solve a nonlinear system on the coarse grid to obtain a rough approximation, and then solve a linearized system resulting from the rough solution to derive a corrected solution on the fine grid.Up to now, this technique has been widely applied to numerically solve many types of nonlinear PDEs, e.g., [8,15,44] for parabolic equations, [32,36] for Darcy-Forchheimer equations, [11] for Navier-Stokes equations, [10,22,23] for time-fractional equations and [5,46] for other nonlinear equations.Compared to traditional implicit-explicit scheme, the main advantages of the two-grid scheme are twofold: (i) if the boundedness of numerical solution in L ∞ norm is not obtained, the temporal stepsize restriction of the method under the local Lipschitz continuous condition on nonlinear term usually is only related to H instead of h, which is much more weaker, see [8,46]; (ii) the two-grid method which treats the nonlinearity on the coarse grid and solving the linear system on the fine grid [33], is much more stable and accurate than the implicit-explicit one when solving the nonlinear PDEs whose solutions change rapidly with respect to time.In this paper, we will investigate the numerical stability of these two methods for general semilinear parabolic PDEs by carrying out representative numerical examples.
Unlike the finite element method, which generates pointwise solution in space and thus is easy to develop two-grid finite element method [7,11,45], the solution yielded by finite difference method is only on discrete grids, and therefore, an appropriate accuracy-preserving mapping operator from the coarse-grid function space to the fine-grid function space is required to construct and analyze the two-grid difference method, e.g., piecewise linear/bilinear interpolation for second-order two-grid difference schemes [8,23,36,46].However, due to the lack of corresponding analysis on appropriate high-order mapping operator, the work about high-order two-grid difference method is meager and, in fact, numerical analysis is also lack.This motivates us to develop high-order two-grid difference scheme for semilinear parabolic PDEs with general boundary condition, e.g., Dirichlet or periodic boundary condition, by introducing and analyzing appropriate mapping operator.
For many time-dependent PDEs, e.g., Allen-Cahn equations [29], whose solutions admiting multiple time scales, adaptive temporal stepsize strategies [24,34] are heuristic and available methods to improve accuracy or efficiency.Due to its strong stability, variable-step BDF2 method is practically valuable for stiff or differential-algebraic problems [13,37].However, compared to those one-step methods, such as the backward Euler and Crank-Nicolson schemes, the numerical analysis of nonuniform BDF2 method would be challenging.In particular, for a linear parabolic problem, [4] proved that, if 0 < r k := τ k /τ k−1 ≤ 1.868 with τ k := t k − t k−1 the kth temporal stepsize, the variable-step BDF2 scheme is zero-stable and second-order convergence containing a prefactor exp(CΓ n ), where Γ n := n−2 k [r k − r k+2 ] + with [x] + the positive part of x.Recently, by introducing a generalized discrete Grönwall inequality, Chen et al. [6] circumvented such a prefactor in error analysis under a little stronger step-ratio restriction 0 < r k ≤ 1.53.In [43], the authors developed an implicit-explicit BDF2 method with variable stepsize for the parabolic partial integro-differential equations and proved its stability and convergence with 0 < r k ≤ 1.91.The authors in [29] considered the fully-implicit BDF2 scheme for the Allen-Cahn equation and established the maximum-norm stability under r k < 1 + √ 2 by developing a novel kernel recombination and complementary technique.To analyze the variable-step BDF2 scheme for linear reaction-diffusion equations, Liao and Zhang [30] introduced a new concept, namely, discrete orthogonal convolution (DOC) kernels, and they improved the unconditional stability in the L 2 norm to r k ≤ 3.561.Subsequently, with the help of DOC kernels and corresponding convolution inequalities, there is a great progress on the stability and error estimates of variable-step BDF2 method for nonlinear PDEs under r k < 3.561 [26,27,39] and the further improved step-ratio restriction r k < 4.8645 [9,20,25], respectively.Among all the variable-step BDF2 methods for nonlinear PDEs in the literature mentioned above, they treat the nonlinear terms fully or partially implicit, in which a nonlinear iteration must be implemented at each time step.Very recently, Zhao et al. [47] presented a linearized variable-step BDF2 scheme for solving nonlinear parabolic equation, and they proved the unconditional error estimate under r k < 4.8645 and the maximum temporal stepsize τ ≤ C 1 √ N by adopting the error splitting approach.After then, Li et al. extended this method to solve a nonlinear Ginzburg-Landau equation [42] and coupled Ginzburg-Landau equations [21] under the same conditions.In [31], a positivity-preserving and energy stable BDF2 scheme with variable stepsize was developed for the Cahn-Hilliard equation with nonlinear logarithmic potential, and convergence analysis in L 2 norm was established under τ ≤ Ch.For solving gradient flow problems, Hou and Qiao [14] proposed an unconditionally energy stable implicit-explicit BDF2 scheme with variable temporal stepsize using the SAV method, and derived the error estimates under the mild restriction on the adjacent temporal stepsize ratio r k < 4.8645.Our goal is to construct and analyze an efficient high-order two-grid difference algorithm with nonuniform BDF2 method for the nonlinear parabolic equation (1.1)- (1.3).Compared to the existing literature, our contributions are mainly threefold: • An efficient variable-step two-grid fourth-order compact difference method is proposed, by using compact difference scheme, two-grid method, variable-step BDF2 formula as well as a developed piecewise bi-cubic Lagrange interpolation operator.
• Under the local continuous condition imposed on the nonlinear term, see (3.2), by adopting the DOC kernels and a cut-off technique, we rigorously prove the unique solvability and error estimate for the nonlinear compact difference scheme, under the temporal stepsize ratio restriction r k < 4.8645 and the maximum temporal stepsize condition τ = o(h ).Furthermore, by discussing the boundedness of the proposed piecewise bi-cubic Lagrange interpolation operator with periodic and Dirichlet boundary conditions, optimal error estimate of the two-grid compact difference scheme is established, under a weaker maximum temporal stepsize condition τ = o(H 1 2 ) and coarse-fine-grid condition H = O(h 4/7 ).
• Several numerical experiments are presented to illustrate the effectiveness and efficiency of the proposed variable-step (adaptive) two-grid compact difference method, and comparisons of computational efficiency and stability with other schemes are also tested.
The remainder of this paper is organized as follows.In Section 2, we introduce and analyze the high-order mapping operator between the coarse-grid function space and fine-grid function space.In Section 3, we first propose a nonlinear compact difference scheme with variable-step BDF2 method for the semilinear parabolic equation subject to Dirichlet boundary condition, and then rigorously prove the unique solvability and convergence, based on which we construct an efficient variable-step two-grid compact difference scheme in Section 4, and optimal-order error analysis under a weaker maximum temporal stepsize condition and a coarse-fine-grid condition is derived.Moreover, in Section 5, the developed methods and techniques are extended to the context of periodic boundary condition.Several numerical experiments are presented to demonstrate the accuracy and efficiency of the proposed method in Section 6.Finally, some concluding remarks are drawn in the last section.

High-order mapping operator between two grids
In this section, we first propose and analyze a high-order mapping operator between two grids based on Lagrange interpolation, which plays a significant role in the construction and numerical analysis of high-order two-grid difference method in the subsequent sections.

Some notations
Given two positive integers N H x and N H y , we define a uniform coarse grid ( x and 0 ≤ j ≤ N H y , with corresponding coarse mesh sizes y and define a uniform fine grid (x h i , y h j ) := (ih x , jh y ) for 0 ≤ i ≤ N h x and 0 ≤ j ≤ N h y , with corresponding fine mesh sizes h x := H x /M x and h y := H y /M y .Denote H := max{H x , H y } and h := max{h x , h y }.
Let ωκ := (i, j) | 0 ≤ i ≤ N κ x , 0 ≤ j ≤ N κ y , ω κ := ωκ ∩ Ω and ∂ω κ := ωκ ∩ ∂Ω denote the sets of spatial grids, where κ = H or h.Accordingly, we define the following discrete spaces of grid functions For any w, q ∈ V κ , we introduce the following notations Similarly, the notations d κ,y w i, j+ 1   2   , d 2 κ,y w i, j and A κ,y w i, j can be defined.Furthermore, we denote x and A κ := A κ,x A κ,y .Besides, we also introduce the discrete inner products and corresponding discrete L 2 and L ∞ norms It is easy to check that A κ w κ ≤ w κ for any w ∈ V κ .Moreover, two well-known and useful lemmas are listed below.

Piecewise bi-cubic Lagrange interpolation
An important tool used in the construction of high-order two-grid method is the local highorder Lagrange interpolation from coarse-grid space to fine-grid space.We shall present and discuss its properties in this subsection.
We first define the one-dimensional piecewise cubic Lagrange interpolation along x-direction.
x − 1, we use φ x i,s (x) 3 s=0 to represent the cubic Lagrange interpolation basis functions.For 1 ≤ i ≤ N H x − 2, φ x i,s (x) is defined as (2.1) For i = 0, i.e., x ∈ (x H 0 , x H 1 ), we define φ x 0,s (x) := φ x 1,s (x); and for i x −2,s (x).Then, for any continuous function w(x), the piecewise cubic Lagrange interpolation operator Π H,x along x-direction is defined as ( where w i = w(x H i ) for 0 ≤ i ≤ N H x .Similarly, we can define the cubic Lagrange interpolation basis functions φ y j,s (y) 3 s=0 and corresponding piecewise cubic Lagrange interpolation operator Π H,y along y-direction.Therefore, the piecewise bi-cubic Lagrange interpolation operator Π H can be defined as the tensor product of the one-dimensional piecewise cubic Lagrange interpolation operators in two directions, that is, Π H := Π H,y Π H,x .When no confusion caused, below we denote (Π H w) i, j = Π H w i, j for (i, j) ∈ ωh and w ∈ V H .

Lemma 2.3 ([35]
). Assume that w ∈ W 4,∞ (Ω), there exists positive constants C 1 and C 2 , independent of H x and H y , such that Next, the bounds for φ x i,s (x) 3 s=0 are given in the following lemma in order to support the proof of boundedness conclusions for the operator Π H . Lemma 2.4.The cubic Lagrange interpolation basis functions φ x i,s (x) 3 s=0 are bounded, i.e., for Note that its first derivative which, together with the definition of φ x i,0 (x), leads to the first conclusion.The remaining conclusions can be similarly proved.
✷ Next, we apply Lemma 2.4 to establish the boundedness results in the discrete L 2 and L ∞ norms for the piecewise bi-cubic Lagrange interpolation operator Π H . Lemma 2.5.For any w ∈ V 0 H , the following estimate holds Proof.Denote ξ i, j = Π H,x w i, j , and then (2.4) We use the identity 4 i=1 a i 2 ≤ 4 4 i=1 a 2 i and the definition (2.2) of cubic Lagrange interpolation operator Π H,y to obtain which, together with Lemma 2.4, yields Analogous to the process (2.5)-(2.6),we can easily derive (2.7) 9) which further implies For fixed 1 ≤ j ≤ N H y − 1, we replace {ξ, y, j} with {w, x, i} in (2.4)-(2.10) to similarly obtain Consequently, we have 12) The proof is completed.✷ Lemma 2.6.For any w ∈ V 0 H , the following estimate holds y , then the triangle inequality and Lemma 2.4 give us and Thus, with a similar treatment to the above estimates, we derive which, together with (2.13)-(2.15),completes the proof.✷

A nonlinear variable-step compact difference scheme
In this section, combined with the variable-step BDF2 method, we are committed to establishing a nonlinear compact difference scheme for the semilinear parabolic equation (1.1)-(1.3)enclosed with Dirichlet boundary condition.Meanwhile, we shall develop the corresponding error estimates by imposing the following regularity assumption By the Sobolev embedding theorem, the above assumption implies that the solution of problem (1.1)-(1.3) is bounded, i.e., there exist two constants m and M such that Furthermore, for a fixed small δ > 0, we denote In this paper, we only assume that the nonlinear term f ∈ C 2 (B δ ) holds locally on B δ , that is,

Nonlinear compact difference scheme
To construct variable temporal stepsize schemes, we consider a nonuniform temporal mesh partition 0 Denote the maximum temporal stepsize τ := max 1≤k≤N τ k and adjacent temporal stepsize ratio In particular, when n = 1, we use the BDF1 (i.e., backward Euler) formula D 1 w 1 = 1 τ 1 ∇ τ w 1 for the first time level discretization.Let r 1 = 0, we then rewrite the above variable-step BDF formula as a unified discrete convolution summation where the discrete convolution kernels b (n) n−k are defined by b (1)  0 := 1/τ 1 and Next, we introduce the discrete orthogonal convolution (DOC) kernels of b (n) n−k by [9,25,30] where By exchanging the summation order and using definition (3.5), it is easy to check that Let U n i, j := u(x h i , y h i , t n ) be the exact nodal solutions.Then, we apply the compact difference operator A h in space and variable-step BDF method in time to get the following equation Under the regularity condition (3.1), by using the Taylor expansion and well-known Bramble-Hilbert Lemma, it is easy to see Let u n i, j be the finite difference approximations to U n i, j , then we drop the local truncation errors in (3.7) to obtain the nonlinear compact difference scheme In the following, we aim to prove the unique solvability and error analysis for the nonlinear variable temporal stepsize compact difference scheme (3.10).For this purpose, we introduce an auxiliary solution ūn = {ū n i, j } satisfying the following scheme subject to the same initial and boundary conditions ( It is obvious that f (w) is globally Lipschitz continuous on R with Lipschitz constant K f .Below we shall show the unique solvability and error analysis for the nonlinear auxiliary scheme (3.11) instead of (3.10).Several useful lemmas are presented below which will be used for our theoretical analysis.

Unique solvability
It is noticed that in [29], the authors proposed a nonlinear variable-step BDF2 scheme for the Allen-Cahn equation with a polynomial type double-well nonlinear potential, and they proved its unique solvability by showing that the solution of the nonuniform BDF2 scheme is equivalent to a minimization problem with strictly convex energy functional.In this subsection, the unique solvability of the auxiliary BDF2 scheme (3.11) for semilinear parabolic equations with general nonlinearity will be discussed by the Browder's fixed point theorem (see e.g.[1]).By homogenization treatment, it suffices to consider the corresponding homogeneous case.
Theorem 3.4.The auxiliary nonlinear compact difference scheme (3.11) Proof.Note that (3.11) can be equivalently rewritten as where G n i, j is defined as Denote the mapping Then, taking the inner product of T ūn with ūn in the sense of (•, •) h gives us Next, we estimate (3.16) term-by-term.For S 2 , utilizing Cauchy-Schwarz inequality yields Moreover, by summation by parts and homogeneous boundary conditions, and noting that A h,x and A h,y are self-adjoint and positive definite operators, we have Besides, we use Cauchy-Schwarz inequality to get which further, together with the globally Lipschitz continuous of f and Lemma 2.1, leads to Then inserting the above estimates into (3.16)gives and subtract these two equations to get (3.17) In (3.17), we let n = m, and then, multiplying it by the DOC kernels θ (k) k−m and summing m from 1 to k yields where we have used the orthogonal identity (3.5).Furthermore, taking the inner product of the above equation with 2v k , and summing the resulting equation from k = 1 to n gives where v 0 = 0 has been used.Due to A h,x and A h,y are self-adjoint and positive definite operators, there exist η x and η y such that A h,x = η 2 x and A h,y = η 2 y .Thus, by summation by parts and homogeneous boundary conditions, the second left-hand side term of (3.19) could be rewritten as where the positive semi-definiteness of the DOC kernels θ (k) k−m (cf.Lemma 3.1) has been used in the last inequality.
Using Cauchy-Schwarz inequality and Lemma 2.1 we see which, together with the global Lipschitz continuous property of f , gives . Then, the above inequality yields which, eliminating v n * A,h from both sides and using Lemma 3.2, implies Then, for τ n ≤ 1/(4 √ 3K f ), an application of the discrete Grönwall inequality in Lemma 3.3 gives which proves the uniqueness of the solutions to scheme (3.11).✷

Error analysis
Denote ēn := U n − ūn and e n := U n − u n for 0 ≤ n ≤ N. It is easy to see that ē0 = e 0 = 0. Next, we shall first give an error estimates in the discrete L 2 norm between the exact solution and the numerical solution yielded by the auxiliary BDF2 scheme (3.11)., the following estimate holds for the auxiliary nonlinear compact difference scheme (3.11) where Proof.Subtracting (3.11) from (3.7), we can get the following error equation for 1 ≤ i ≤ N h x − 1 and 1 ≤ j ≤ N h y − 1. Analogous to the proof of (3.17)-(3.18), it follows from the orthogonality (3.5) of the DOC kernels {θ (k)  k−m } that which, by taking the inner product with 2ē k and summing the resulting equality from k = 1 to n, gives With a similar treatment to (3.20)-(3.21), the above inequality yields Similarly, if ēn * A,h = max 1≤k≤n ēn A,h , using Lemma 3.2 we have (3.28) For the last term of (3.28), we exchange the order of summation and utilize Lemma 3.2 to obtain which, together with the triangle inequality and (3.8)-(3.9),gives Now, insert the above inequality into (3.28),for τ ≤ 1/(4 Then, applying the discrete Grönwall inequality to (3.30) and using Lemma 2.1 yields which implies the theorem.✷ Finally, we would like to give the L 2 norm error estimates between the exact solution and the numerical solution yielded by the variable-step BDF2 compact scheme (3.10).Note that by Theorem 3.6 and Lemma 2.2, we see ), and for τ, h sufficiently small, it holds that ūn i, j ∈ B δ , which further implies f ( ūn i, j ) = f ( ūn i, j ) and thus in this context ūn i, j ≡ u n i, j .We summarize these conclusions in the following theorem.Theorem 3.7.Under the conditions in Theorem 3.6 and if the maximum temporal stepsize satisfies τ = o(h 1 2 ), the nonlinear compact difference scheme (3.10) admits a unique solution satisfying 4. An efficient variable-step two-grid compact difference scheme In order to solve the nonlinear system (1.1)-(1.3)efficiently, we shall propose a two-grid compact finite difference scheme based on the variable-step BDF2 method and the piecewise bi-cubic Lagrange interpolation developed in Section 2 as follows.
Step 1.On the coarse grid, solve a small-scale nonlinear compact difference scheme to find a rough solution u n H = {u n H,i, j } by subject to the initial and boundary conditions (1.2)-(1.3)defined on coarse grid.
Step 2. On the fine grid, solve a large-scale linearized compact difference scheme to produce a corrected solution u n h = {u n h,i, j } based on the rough solution u n H in Step 1 by subject to the initial and boundary conditions (1.2)-(1.3),where F n i, j represents a Newton linearization from coarse grid to fine grid defined as On the coarse grid, denote e n H = U n − u n H for 0 ≤ n ≤ N. Analogous to Theorem 3.7, we can immediately reach the following conclusions.), the nonlinear compact difference scheme (4.1) defined on the coarse grid admits a unique solution satisfying Based on this theorem and Lemmas 2.5-2.6, the following corollary can be derived.
Corollary 4.2.Under the conditions in Theorem 4.1, the numerical solution u n H of the nonlinear compact difference scheme (4.1) satisfies and consequently the interpolation solution Π H u n H is bounded on the fine grid and satisfies where Proof.We perform the splitting where the first right-hand side term could be bounded as which completes the proof of estimate (4.4).Analogous to the proof of (4.4), by Lemmas 2.2, 2.3, 2.6 and Theorem 4.1, we can conclude that Furthermore, an application of the triangle inequality yields which leads to a much worse restriction τ = o(h ).On the other hand, in practical computation the condition τ = o(H 1 2 ) is not too severe, as the coarse grid size H is large enough compared with the fine grid size h, and a mild choice τ = O(H 2 ) suggested by Theorem 4.1 can naturally satisfy the restrictive condition.
At last, we shall give an error estimate for e n h = U n − u n h on the fine grid for the linearized compact difference scheme (4.2).Theorem 4.4.Under the conditions in Theorem 4.1, the following error estimate holds for the solution of the two-grid compact difference scheme (4.1)-(4.2) where Proof.For the linearized scheme (4.2) on the fine grid, we can get a very similar error equation which, together with a similar treatment to (3.24)-(3.26),leads to (4.8) Compared with the estimate (3.26) in Theorem 3.6, the only difference lies in the treatment of the nonlinear term.We apply Taylor expansion of f for some constant µ m i, j between U m i, j and [Π H u m H ] i, j .Then, subtract F m i, j in (4.3) from this equation, we have . Furthermore, we apply Lemma 2.1, Cauchy-Schwarz inequality and Corollary 4.2 to derive (4.9) The other terms in (4.8) could be estimated similarly as in the proof of Theorem 3.6.Then, substituting (4.9) into (4.8)gives us As was done before, choosing an integer n * (1 Therefore, for τ ≤ 1/(4 √ 3K f ), we use Lemma 3.2 and estimate (3.29) to obtain from which an application of the discrete Grönwall inequality in Lemma 3.3 and Lemma 2.1 yield the following estimate The proof is completed.✷

Extension to periodic boundary condition
In this section, we extend the ideas and derivations in previous sections to the semilinear parabolic equation (1.1)-(1.3)with periodic boundary condition.Firstly, we denote the following spaces of grid functions on grids ωκ Furthermore, for any grid functions w, q ∈ V p κ , the discrete inner product and corresponding norms are redefined as (w, q) κ = κ x κ y N κ x i=1 N κ y j=1 w i, j q i, j , w κ = (w, w) κ , w A,κ = (A κ w, w) κ , and a useful lemma is listed below.Lemma 5.1 ( [28]).For any w ∈ V p κ , we have 2 3 w κ ≤ w A,κ ≤ w κ .In the context of periodic boundary case, we shall still adopt the piecewise bi-cubic Lagrange interpolation operator defined in Section 2 to construct high-order two-grid difference scheme.However, a small modification of the proof of Lemmas 2.5-2.6 should be given to show the boundedness conclusions of the piecewise bi-cubic Lagrange interpolation operator under the discrete L 2 and L ∞ norms.Lemma 5.2.For any w ∈ V p H , the following estimate holds Proof.Denote ξ i, j = Π H,x w i, j , and then Π H w 2 h can be rewritten as where I 1 , I 2 are defined in Lemma 2.5 and While, the application of Lemma 2.4 shows This together with (2.6)-(2.7)implies Finally, the claimed result can be derived immediately by following a completely similar process as (2.10)-(2.12).✷ Lemma 5.3.For any w ∈ V p H , the following estimate holds where Proof.Similar as the proof of Lemma 2.6, estimate (2.13) holds for M y ≤ j * ≤ (N H y − 1)M y .Besides, if 1 ≤ j * ≤ M y , we have and similarly, if (N H y − 1)M y ≤ j * ≤ N h y , we have Finally, analogous to (2.16), we see x w i,k holds, which completes the proof.✷ Now, an efficient two-grid fourth-order compact difference scheme for model (1.1)-(1.3) is proposed similarly as follows.
Step 1.On the coarse grid, solve a small-scale nonlinear compact finite difference scheme to find a rough solution u n H by subject to the initial condition (1.2) and periodic boundary condition.
Step 2. On the fine grid, solve a large-scale linearized compact difference scheme to produce a corrected solution u n h based on the rough solution u n H in Step 1 by subject to the initial condition (1.2) and periodic boundary condition.Following the proofs of Corollary 4.2, Theorems 4.1 and 4.4, together with Lemmas 5.1-5.3, the unique solvability and error estimates for the two-grid algorithm (5.3)-(5.4)can be proved very similarly, and we skip the detailed proof here.) and τ ≤ 1 6K f , the two-grid compact difference scheme (5.3)-(5.4)admits a unique solution satisfying where C 12 := 3 exp(6K f T ) max{6C 5 + 3C 6 T + C 9 C 10 K f T, 9C 7 T }.

Numerical examples
In this section, we shall present several numerical experiments to test the effectiveness and efficiency of the variable-step two-grid compact difference scheme.In the computation, a Newtontype iterative procedure with tolerance error 1.0 × 10 −13 is performed to solve the nonlinear algebra systems at each time level.
Figure 1 displays the evolution of these solutions with time, in which it can be clearly observed that the solution in Case I changes most smoothly while the solution in Case III changes most sharply with respect to time.Firstly, we set N h x = N h y and M x = M y = 10 and adjust N and N h x to investigate the spatial convergence of these three methods.Numerical results for Cases I-III are listed in Tables 1-3 respectively, which indicates the fourth-order spatial accuracy for both the nonlinear method and the two-grid method.But the implicit-explicit method is only successfully implemented for Case I and fails for the other two cases, which may be caused by improper treatment of the time dependence and nonlinearity.
Secondly, to test the temporal convergence rates, we fix N h x = N h y = 10N H x = 10N H y = 300 and present the numerical results with respect to N in Table 4 for Case I. We can observe that these three methods all have second-order temporal accuracy as proved.However, when the solution changes dramatically over time (e.g.Case II or III), the implicit-explicit method becomes unstable (|u| → ∞), while both the nonlinear method and two-grid method can generate the desired numerical solutions with the same magnitude accuracy, as seen in Tables 5-6.Furthermore, the numerical results show that if we further refine the temporal grids (e.g., changing N from 256 to 512 in Table 5), the implicit-explicit method may produce a correct result.But as shown in Table 6, its stability requirement for the temporal grid is quite related to the smoothness of the solution with respect to time.It is seen that the implicit-explicit discretization has much more strict stability condition compared to the other two methods when the solution u changes sharply with respect to t, for which a convincing explanation is that approximating the nonlinear term via solutions at previous time levels may leads to inaccuracy in this context.

Accuracy and efficiency tests on variable-step temporal grids
To check the accuracy and efficiency on variable-step temporal grids, we consider model (1.1)-(1.3)on (0, 1) 2 × (0, 1] with c = 1 8π 2 and f (u) = u − u 3 .The linear part g is determined such that the exact solution is u(x, y, t) = sin(t) sin(2πx) sin(2πy).The variable-step temporal grids are generated randomly by where θ k is randomly drawn from the uniform distribution on the interval (1/4.8645, 1) such that the adjacent temporal stepsize ratio r k < 4.8645.
As the implicit-explicit scheme (6.1) may generate wrong results, we only test the nonlinear scheme (3.10) and two-grid scheme (4.1)-(4.2).We firstly test the errors and convergence rates in spatial and temporal directions for both methods with N h x = N h y and M x = M y = 5.The corresponding numerical results are listed in Tables 7-8 respectively, which indicates the fourthorder accuracy in space and second-order accuracy in time as proved in Theorems 3.6 and 4.4.Moreover, we compare the CPU times consumed by the two methods in Table 9 under N h x = N h y and M x = M y , in which we can clearly observe that the proposed two-grid method has significantly improved the computational efficiency, for example, it takes about two and a half hours for the implementation of the nonlinear scheme when N = N h x = 480, while the two-grid scheme consumes only about one hour to desire the same error!3) on (0, 1) 2 × (0, 4] with c = 1 and f (u) = sin u, and the exact solution is chosen as u(x, y, t) = 1 + 20e −40(t−1) 2 + 30e −60(t−4) 2 sin(2πx) sin(2πy).Figure 2 (left) depicts the evolution of solution with respect to time at fixed point (0.25, 0.25), which consists of two peaks and admits multiple time scales.Therefore, the variable-step twogrid scheme based on the adaptive temporal stepsize strategy [17,25] will be adopted to improve the temporal accuracy where ∂ τ u n := ∇ τ u n /τ n and r max = 4.8 which satisfies the restrictions in Theorems 3.6 and 4.4.
Here τ min and τ max are the pre-determined minimum and maximum temporal stepsize and η is a pre-chosen parameter.In this example, we uniformly set τ max = 0.2, η = 500 and gradually reduce τ min to generate grids with distinct stepsizes.We select N h x = N h y = 250 to test the temporal convergence rates yielded by the two-grid scheme on uniform and adaptive temporal grids with M x = M y = 10.The numerical results are presented in Table 10, which shows that the standard two-grid scheme is second-order accurate in time, while the variable-step two-grid scheme based on the adaptive temporal stepsize strategy (6.3) has better convergence and much smaller errors with the same number of grids.In other words, the standard two-grid scheme on uniform temporal grids requires more temporal steps, and of course more CPU times, to generate numerical solutions with the same magnitude accuracy as adaptive method which uses less temporal steps.The reason can be more intuitively observed from Figure 2 (right), which shows that when the solution varies sharply, small temporal stepsizes are adaptively created to capture the fast evolution process, while otherwise large temporal stepsizes are generated to accelerate the time integration.

Application to phase-field Allen-Cahn equation
In this subsection, we consider the following Allen-Cahn equation with a polynomial doublewell potential, subject to periodic boundary condition where ε is the interaction length that describes the thickness of the transition boundary between materials.It is well known that the energy dissipation law [14,38]

E[u](t) ≤ E[u](s), ∀t > s
holds for the Allen-Cahn equation, where E[u](t) represents the Lyapunov energy functional, namely Moreover, it has been observed that the evolution of the energy E[u](t) usually involves both fast and slow stages of change in the long time simulation.Thus, it is highly desirable for numerical methods to preserve the discrete energy dissipation law on the nonuniform temporal grids.Define the discrete energy functional E[u n ] as In this test, we would also compare effectiveness and efficiency of the high-order nonlinear difference scheme (3.10), the two-grid difference scheme (4.1)-(4.2) and the adaptive two-grid difference scheme using a similar adaptive temporal stepsize strategy (6.3) by replacing Example 6.1.In this example, we set ε = 0.02 and apply these three methods to simulate the merging of four bubbles with an initial condition In this simulation, a 384 × 384 uniform mesh is taken to discretize the spatial domain Ω = (−1, 1) 2 and the ratio of coarse-fine grids are setted as M x = M y = 3.We start with the modeling of the solution by the nonlinear scheme and two-grid scheme with a constant temporal stepsize τ = 0.1 until time T = 100.Parameters in the adaptive temporal stepsize strategy (6.3) are selected as τ min = 0.1, τ max = 1 and η = 3200.In Figure 3, it displays a comparison on the evolution of solution snapshots among the three methods, in which the gradually merging and shrinking process of the initial four-drops over time can be clearly observed.As can be seen in the figures, there seems no distinguishable differences among these methods.Next, we investigate the efficiency of the two-grid method on uniform grid and on adaptive grid for long time modeling.As seen in Figure 4 (left), the evolution of the free energy with respect to time for these three methods coincide, which consists very well with the corresponding results in [29].Moreover, Table 11 indicates that the adaptive two-grid scheme has significant advantage in computational efficiency over the other two schemes.For example, it takes about 10 hours for the implementation of the nonlinear scheme up to T = 100, while the two-grid method using uniform temporal grid consumes about 5 hours.What is even more amazing is that the developed adaptive two-grid method using variable-step temporal grid takes only about one hour!In fact, the total number of adaptive temporal steps is only 156, while it takes 1000 steps for the uniform grid.Finally, the adaptive temporal stepsize curve of the adaptive two-grid method is plotted in Figure 4 (right), which also demonstrates the superiority of the variable-step two-grid compact difference scheme.Example 6.2.In this example, we consider the coarsening process governed by the Allen-Cahn equation with the model parameter ε = 0.01 and computational domain Ω = (0, 1) 2 .Here we choose a random initial condition u 0 (x, y) = −0.05+ 0.1 × rand(x, y).
This simulation is performed under N h x = N h y = 384 and M x = M y = 3. Due to the fact that the initial values are randomly given, we provide the startup values on coarse grid for the twogrid scheme by implementing the nonlinear algorithm up to T = 0.5.In Figure 5, it displays the 29      evolution of the coarsening dynamic for nonlinear and two-grid compact schemes with different time strategies.It is observed that the nonlinear scheme with large uniform temporal stepsize τ = 1 yields inaccurate solution u, while the adaptive two-grid scheme gives the correct coarsening pattern which is consistent with the results obtained by the nonlinear and two-grid methods with small uniform temporal stepsize τ = 0.01.In Figure 6, we depict the evolution of discrete energies and temporal stepsizes with respect to time, which shows that the energy dissipation of the adaptive two-grid method agrees very well with the the nonlinear and two-grid methods using small uniform temporal stepsize.Moreover, the efficiency of the proposed two-grid scheme using the adaptive temporal stepsize strategy can also be seen from Figure 6 (right) and Table 12.For example, the adaptive two-grid method using variable-step temporal grid costs only 33 minutes for time marching to T = 100, while the two-grid method using uniform temporal stepsize τ = 0.01 consumes more than 8 hours, even worse the implementation of the nonlinear scheme runs nearly 19 hours.

Concluding remarks
High-order two-grid difference scheme for nonlinear PDEs are rarely studied in existing literature due to, e.g., the lack of the appropriate accuracy-preserving mapping operator.To address this issue, we introduce a piecewise bi-cubic Lagrange interpolation operator between two grids, and discuss its boundedness under L 2 and L ∞ norms.Moreover, to effectively solve the nonlinear PDEs whose solutions may admit multiple time scales, variable-step temporal discretization methods, in particular, the variable-step multistep methods are natrually and valuable to improve accuracy for stiff problems.However, its numerical analysis is much more challenging than the single-step methods.As an illustration, combined with the variable-step BDF2 scheme, an efficient high-order two-grid difference method is developed for the semilinear parabolic equation with Dirichlet or periodic boundary conditions.The unique solvability of the nonlinear problem on coarse grid is shown by Browder's fixed point theorem.Moreover, with the help of DOC kernels and the boundedness of the high-order mapping operator, optimal-order error estimates for the two-grid method on both coarse and fine grids are rigourously proved under r k := τ k /τ k−1 < 4.8645 and the maximum temporal stepsize condition τ = o(H 1 2 ), where the cut-off technique is used to reduce the regularity requirement on the nonlinear term f to the local Lipschitz continuous condition.Several numerical examples are carried out to confirm the theoretical findings.

Theorem 4 . 1 .
Under the conditions in Theorem 3.6 and if the maximum temporal stepsize satisfies τ = o(H 1 2 3, and due to the linear property of the piecewise bi-cubic Lagrange interpolation operator with respect ro the interpolated function, Lemma 2.5 and Theorem 4.1 give us

Figure 1 :
Figure 1: Evolution of exact solution with time at point (0.25, 0.25)

Figure 2 :
Figure 2: Evolution of solution with time (left) at (0.25, 0.25) and temporal stepsize (right) of the two-grid scheme and adaptive two-grid scheme

Figure 3 :
Figure 3: Solution snapshots of Allen-Cahn equation at t = 1, 10, 50, 100 (from left to right) yielded by nonlinear scheme, two-grid scheme and adaptive two-grid scheme

Figure 4 :
Figure 4: Evolutions of energy (left) and time steps (right) for the nonlinear scheme, two-grid scheme and adaptive two-grid scheme until time T = 30

Figure 5 :
Figure 5: Solution snapshots of coarsening dynamics for Allen-Cahn equation at t = 1, 20, 50, 80, 100 (from left to right) yielded by nonlinear scheme, two-grid scheme and adaptive two-grid scheme

Figure 6 :
Figure 6: Evolutions of energy (left) and time steps (right) for the nonlinear scheme, two-grid scheme and adaptive two-grid scheme until time T = 100 Consequently, Browder's fixed point theorem shows there exists a ūn ∈ V 0 h such that T ūn = 0 and ūn The solution of the auxiliary nonlinear compact difference scheme (3.11) is unique if the adjacent temporal stepsize ratios r k satisfy 0 < r k < 4.8645 and the maximum temporal stepsize τ ≤ 1 The argument is by contradiction.Suppose that there are two solutions ūn 1 and ūn 2 satisfying(3.11)andthe same initial and boundary conditions (1.2)-(1.3),i.e., Corollary 4.2 declares that the interpolation solution Π H u n H is bounded under the discrete L ∞ norm on the fine grid, if the maximum temporal stepsize satisfies τ = o(H 1 2 ).This restriction is compatible with previous work, see Refs.[8, Theorem 4.2] and [46, Theorem 4.2].Notice that the boundedness of Π H in the sense of L ∞ norm in Lemma 2.6 plays an important role in the proof of (4.5)-(4.6).Actually, if Lemma 2.6 does not hold, we have which completes the proof of(4.6).✷ Remark 4.3.

Table 10 :
Temporal convergence of two-grid scheme and adaptive two-grid scheme

Table 11 :
CPU times and the total number of temporal steps for three schemes

Table 12 :
CPU times and the total number of temporal steps for three schemes