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

Bingyin Zhang<EMAIL_ADDRESS>Hongfei Fu<EMAIL_ADDRESS>School
of Mathematical Sciences, Ocean University of China, Qingdao, Shandong 266100,
China

###### Abstract

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
$r_{k}:=\tau_{k}/\tau_{k-1}<4.8645$ and the maximum temporal stepsize
satisfies $\tau=o(h^{\frac{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 $r_{k}<4.8645$ and a weaker
maximum temporal stepsize condition $\tau=o(H^{\frac{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=\mathcal{O}(h^{4/7})$. Finally, several numerical experiments are carried
out to demonstrate the effectiveness and efficiency of the proposed scheme.

###### keywords:

High-order two-grid difference algorithm, Variable-step BDF2 method, Piecewise
bi-cubic Lagrange interpolation, DOC kernels, Unique solvability, Error
estimate

##  1 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, 16, 33, 8, 12, 19].
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^{\prime\prime}(w)|\leq K_{f}$ for
$w\in\mathbb{R}$) cannot be imposed or the boundedness of numerical solution
in $L^{\infty}$ 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., $\tau=\mathcal{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, 48, 41] 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, [36, 32] for Darcy-
Forchheimer equations, [11] for Navier-Stokes equations, [23, 22, 10] 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^{\infty}$ 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 [45, 7, 11], 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,
36, 23, 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.

In this paper, to illustrate the application of the proposed high-order two-
grid difference method, we focus on the following semilinear parabolic
equation

$\displaystyle u_{t}(x,y,t)-c\Delta
u(x,y,t)=f(u(x,y,t))+g(x,y,t),\quad(x,y)\in\Omega,\ t\in(0,T],$ (1.1)

where $\Omega=(0,L_{x})\times(0,L_{y})$, subject to the initial condition

$\displaystyle u(x,y,0)=u_{0}(x,y),\quad(x,y)\in\Omega,$ (1.2)

and periodic boundary condition or Dirichlet boundary condition

$\displaystyle u(x,y,t)=\psi(x,y,t),\quad(x,y)\in\partial\Omega,\ t\in[0,T],$
(1.3)

where $c>0$ and $\Delta$ is the Laplacian operator
$\partial^{2}_{x}+\partial^{2}_{y}$, $f(u)$ is the nonlinear term, $u_{0}$,
$\psi$ and $g$ are given smooth functions.

For many time-dependent PDEs, e.g., Allen-Cahn equations [29], whose solutions
admiting multiple time scales, adaptive temporal stepsize strategies [34, 24]
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}:=\tau_{k}/\tau_{k-1}\leq 1.868$ with $\tau_{k}:=t_{k}-t_{k-1}$ the
$k$th temporal stepsize, the variable-step BDF2 scheme is zero-stable and
second-order convergence containing a prefactor $\exp(C\Gamma_{n})$, where
$\Gamma_{n}:=\sum^{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}\leq 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}\leq 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+\sqrt{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}\leq 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$ [27, 26, 39] and the further improved step-ratio
restriction $r_{k}<4.8645$ [9, 25, 20], 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
$\tau\leq C\frac{1}{\sqrt{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 $\tau\leq 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:

  * 1.

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.

  * 2.

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 $\tau=o(h^{\frac{1}{2}})$. 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
$\tau=o(H^{\frac{1}{2}})$ and coarse-fine-grid condition
$H=\mathcal{O}(h^{4/7})$.

  * 3.

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.

##  2 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.

###  2.1 Some notations

Given two positive integers $N^{H}_{x}$ and $N^{H}_{y}$, we define a uniform
coarse grid $(x^{H}_{i},y^{H}_{j}):=(iH_{x},jH_{y})$ for $0\leq i\leq
N^{H}_{x}$ and $0\leq j\leq N^{H}_{y}$, with corresponding coarse mesh sizes
$H_{x}:=L_{x}/N^{H}_{x}$ and $H_{y}:=L_{y}/N^{H}_{y}$. Moreover, for fixed
positive integers $M_{x},M_{y}\geq 2$, denote $N^{h}_{x}:=M_{x}N^{H}_{x}$,
$N^{h}_{y}:=M_{y}N^{H}_{y}$ and define a uniform fine grid
$(x^{h}_{i},y^{h}_{j}):=(ih_{x},jh_{y})$ for $0\leq i\leq N^{h}_{x}$ and
$0\leq j\leq 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 $\bar{\omega}_{\kappa}:=\left\\{(i,j)\mid 0\leq i\leq N^{\kappa}_{x},\
0\leq j\leq N^{\kappa}_{y}\right\\}$,
$\omega_{\kappa}:=\bar{\omega}_{\kappa}\cap\Omega$ and
$\partial\omega_{\kappa}:=\bar{\omega}_{\kappa}\cap\partial\Omega$ denote the
sets of spatial grids, where $\kappa=H$ or $h$. Accordingly, we define the
following discrete spaces of grid functions

$\mathcal{V}_{\kappa}=\left\\{v=\\{v_{i,j}\\}\mid(i,j)\in\bar{\omega}_{\kappa}\right\\}\quad{\rm
and}\quad\mathcal{V}^{0}_{\kappa}=\left\\{v\mid v\in\mathcal{V}_{\kappa}\ \
\text{and}\ \ v_{i,j}=0\ {\rm if}\ (i,j)\in\partial\omega_{\kappa}\right\\}.$

For any $w,q\in\mathcal{V}_{\kappa}$, we introduce the following notations

$\begin{array}[]{l}\displaystyle
d_{\kappa,x}w_{i+\frac{1}{2},j}:=\frac{1}{\kappa_{x}}\left(w_{i+1,j}-w_{i,j}\right),\quad\displaystyle
d^{2}_{\kappa,{x}}w_{i,j}:=\frac{1}{\kappa_{x}}\left(\left[d_{\kappa,x}w\right]_{i+\frac{1}{2},j}-\left[d_{\kappa,x}w\right]_{i-\frac{1}{2},j}\right),\end{array}$
$\begin{array}[]{ll}\mathcal{A}_{\kappa,x}w_{i,j}:=\left\\{\begin{split}&w_{i,j}+\frac{\kappa_{x}^{2}}{12}d^{2}_{\kappa,x}w_{i,j}=\frac{1}{12}(w_{i-1,j}+10w_{i,j}+w_{i+1,j}),&\quad
1\leq i\leq N^{\kappa}_{x}-1,\\\ &w_{i,j},&\quad
i=0,N^{\kappa}_{x}.\end{split}\right.\end{array}$

Similarly, the notations $d_{\kappa,y}w_{i,j+\frac{1}{2}}$,
$d^{2}_{\kappa,y}w_{i,j}$ and $\mathcal{A}_{\kappa,y}w_{i,j}$ can be defined.
Furthermore, we denote $\Delta_{\kappa}:=d^{2}_{\kappa,x}+d^{2}_{\kappa,y}$,
$\Lambda_{\kappa}:=\mathcal{A}_{\kappa,x}d^{2}_{\kappa,y}+\mathcal{A}_{\kappa,y}d^{2}_{\kappa,x}$
and $\mathcal{A}_{\kappa}:=\mathcal{A}_{\kappa,x}\mathcal{A}_{\kappa,y}$.

Besides, we also introduce the discrete inner products

$(w,q)_{\kappa}=\kappa_{x}\kappa_{y}\sum^{N^{\kappa}_{x}-1}_{i=1}\sum^{N^{\kappa}_{y}-1}_{j=1}w_{i,j}q_{i,j},$
$(w,q)_{\kappa,x}=\kappa_{x}\kappa_{y}\sum^{N^{\kappa}_{x}-1}_{i=0}\sum^{N^{\kappa}_{y}-1}_{j=1}w_{i+\frac{1}{2},j}q_{i+\frac{1}{2},j},~{}~{}(w,q)_{\kappa,y}=\kappa_{x}\kappa_{y}\sum^{N^{\kappa}_{x}-1}_{i=1}\sum^{N^{\kappa}_{y}-1}_{j=0}w_{i,j+\frac{1}{2}}q_{i,j+\frac{1}{2}},$
$<w,q>_{\kappa}=\kappa_{x}\kappa_{y}\left[\frac{1}{4}\sum_{i=\\{0,N^{x}_{\kappa}\\}}\sum_{j=\\{0,N^{y}_{\kappa}\\}}+\frac{1}{2}\sum^{N^{x}_{\kappa}-1}_{i=1}\sum_{j=\\{0,N^{y}_{\kappa}\\}}+\frac{1}{2}\sum_{i=\\{0,N^{x}_{\kappa}\\}}\sum^{N^{y}_{\kappa}-1}_{j=1}+\sum^{N^{x}_{\kappa}-1}_{i=1}\sum^{N^{y}_{\kappa}-1}_{j=1}\right]w_{i,j}q_{i,j},$

and corresponding discrete $L^{2}$ and $L^{\infty}$ norms

$\|w\|_{\kappa}=\sqrt{(w,w)_{\kappa}},\quad\interleave
w\interleave_{\kappa}=\sqrt{<w,w>_{\kappa}},\quad\|w\|_{\mathcal{A},\kappa}=\sqrt{(\mathcal{A}_{\kappa}w,w)_{\kappa}},\quad\|w\|_{\kappa,\infty}=\max_{(i,j)\in\bar{\omega}_{\kappa}}|w_{i,j}|.$

It is easy to check that $\|\mathcal{A}_{\kappa}w\|_{\kappa}\leq\interleave
w\interleave_{\kappa}$ for any $w\in\mathcal{V}_{\kappa}$. Moreover, two well-
known and useful lemmas are listed below.

######  Lemma 2.1 ([40])

For any $w\in\mathcal{V}^{0}_{\kappa}$, we have
$\frac{1}{3}\|w\|_{\kappa}^{2}\leq\|w\|_{\mathcal{A},\kappa}^{2}\leq\|w\|_{\kappa}^{2}.$

######  Lemma 2.2 ([40])

For any $w\in\mathcal{V}_{\kappa}$, there exists a positive constant $C_{0}$,
independent of $\kappa$, such that

$\|w\|_{\kappa,\infty}\leq C_{0}\ \kappa^{-1}\|w\|_{\kappa}.$

###  2.2 Piecewise bi-cubic Lagrange interpolation

An important tool used in the construction of high-order two-grid method is
the local high-order 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. For each $x\in(x^{H}_{i},x^{H}_{i+1})$ with $0\leq i\leq
N^{H}_{x}-1$, we use $\big{\\{}\phi^{x}_{i,s}(x)\big{\\}}^{3}_{s=0}$ to
represent the cubic Lagrange interpolation basis functions. For $1\leq i\leq
N^{H}_{x}-2$, $\phi^{x}_{i,s}(x)$ is defined as

$\begin{array}[]{l}\phi^{x}_{i,s}(x):=\left\\{\begin{split}&-\frac{(x-x^{H}_{i})(x-x^{H}_{i+1})(x-x^{H}_{i+2})}{6H^{3}_{x}},\qquad
s=0,\\\
&\frac{(x-x^{H}_{i-1})(x-x^{H}_{i+1})(x-x^{H}_{i+2})}{2H^{3}_{x}},\qquad\ \ \
s=1,\\\
&-\frac{(x-x^{H}_{i-1})(x-x^{H}_{i})(x-x^{H}_{i+2})}{2H^{3}_{x}},\qquad\
s=2,\\\
&\frac{(x-x^{H}_{i-1})(x-x^{H}_{i})(x-x^{H}_{i+1})}{6H^{3}_{x}},\qquad\quad\
s=3.\end{split}\right.\end{array}$ (2.1)

For $i=0$, i.e., $x\in(x^{H}_{0},x^{H}_{1})$, we define
$\phi^{x}_{0,s}(x):=\phi^{x}_{1,s}(x)$; and for $i=N^{H}_{x}-1$, i.e.,
$x\in(x^{H}_{N^{H}_{x}-1},x^{H}_{N^{H}_{x}})$, we define
$\phi^{x}_{N^{H}_{x}-1,s}(x):=\phi^{x}_{N^{H}_{x}-2,s}(x)$. Then, for any
continuous function $w(x)$, the piecewise cubic Lagrange interpolation
operator $\Pi_{H,x}$ along $x$-direction is defined as

$\Pi_{H,x}w(x):=\left\\{\begin{aligned} &\sum^{3}_{s=0}w_{s}\
\phi^{x}_{0,s}(x),&\quad x\in(x^{H}_{0},x^{H}_{1}),\ i=0,\\\
&\sum^{3}_{s=0}w_{i-1+s}\ \phi^{x}_{i,s}(x),&\quad
x\in(x^{H}_{i},x^{H}_{i+1}),\ 1\leq i\leq N^{H}_{x}-2,\\\
&\sum^{3}_{s=0}w_{N^{H}_{x}-3+s}\ \phi^{x}_{N^{H}_{x}-1,s}(x),&\quad
x\in(x^{H}_{N^{H}_{x}-1},x^{H}_{N^{H}_{x}}),\
i=N^{H}_{x}-1.\end{aligned}\right.$ (2.2)

where $w_{i}=w(x^{H}_{i})$ for $0\leq i\leq N^{H}_{x}$.

Similarly, we can define the cubic Lagrange interpolation basis functions
$\big{\\{}\phi^{y}_{j,s}(y)\big{\\}}^{3}_{s=0}$ and corresponding piecewise
cubic Lagrange interpolation operator $\Pi_{H,y}$ along $y$-direction.
Therefore, the piecewise bi-cubic Lagrange interpolation operator $\Pi_{H}$
can be defined as the tensor product of the one-dimensional piecewise cubic
Lagrange interpolation operators in two directions, that is,
$\Pi_{H}:=\Pi_{H,y}\Pi_{H,x}$. When no confusion caused, below we denote
$(\Pi_{H}w)_{i,j}=\Pi_{H}w_{i,j}$ for $(i,j)\in\bar{\omega}_{h}$ and
$w\in\mathcal{V}_{H}$.

######  Lemma 2.3 ([35])

Assume that $w\in W^{4,\infty}(\Omega)$, there exists positive constants
$C_{1}$ and $C_{2}$, independent of $H_{x}$ and $H_{y}$, such that

$\|w-\Pi_{H}w\|_{h,\infty}\leq C_{1}H^{4}\quad{\rm
and}\quad\|w-\Pi_{H}w\|_{h}\leq C_{2}H^{4}.$

Next, the bounds for $\big{\\{}\phi^{x}_{i,s}(x)\big{\\}}^{3}_{s=0}$ are given
in the following lemma in order to support the proof of boundedness
conclusions for the operator $\Pi_{H}$.

###### Lemma 2.4

The cubic Lagrange interpolation basis functions
$\big{\\{}\phi^{x}_{i,s}(x)\big{\\}}^{3}_{s=0}$ are bounded, i.e., for $1\leq
i\leq N^{H}_{x}-2$, we have

$\begin{array}[]{l}\left|\phi^{x}_{i,0}(x)\right|\leq\left\\{\begin{split}&1,\qquad\qquad\
\ \ x\in(x^{H}_{i-1},x^{H}_{i}),\\\ &\frac{\sqrt{3}}{27},\qquad\quad\ \ \
x\in(x^{H}_{i},x^{H}_{i+1}),\\\ &\frac{\sqrt{3}}{27},\qquad\quad\
x\in(x^{H}_{i+1},x^{H}_{i+2}),\end{split}\right.\end{array}\qquad\begin{array}[]{l}\left|\phi^{x}_{i,1}(x)\right|\leq\left\\{\begin{split}&\frac{7\sqrt{7}+10}{27},\qquad\quad
x\in(x^{H}_{i-1},x^{H}_{i}),\\\ &1,\qquad\qquad\qquad\ \ \
x\in(x^{H}_{i},x^{H}_{i+1}),\\\ &\frac{7\sqrt{7}-10}{27},\qquad\ \
x\in(x^{H}_{i+1},x^{H}_{i+2}),\end{split}\right.\end{array}$
$\begin{array}[]{l}\left|\phi^{x}_{i,2}(x)\right|\leq\left\\{\begin{split}&\frac{7\sqrt{7}-10}{27},\qquad\quad
x\in(x^{H}_{i-1},x^{H}_{i}),\\\ &1,\qquad\qquad\qquad\ \ \
x\in(x^{H}_{i},x^{H}_{i+1}),\\\ &\frac{7\sqrt{7}+10}{27},\qquad\ \
x\in(x^{H}_{i+1},x^{H}_{i+2}),\end{split}\right.\end{array}\qquad\begin{array}[]{l}\left|\phi^{x}_{i,3}(x)\right|\leq\left\\{\begin{split}&\frac{\sqrt{3}}{27},\qquad\
\quad\ x\in(x^{H}_{i-1},x^{H}_{i}),\\\ &\frac{\sqrt{3}}{27},\qquad\quad\ \
x\in(x^{H}_{i},x^{H}_{i+1}),\\\ &1,\qquad\qquad\
x\in(x^{H}_{i+1},x^{H}_{i+2}).\end{split}\right.\end{array}$

Proof. Consider auxiliary function

$Z(x)=\left(x-x^{H}_{i}\right)\left(x-x^{H}_{i+1}\right)\left(x-x^{H}_{i+2}\right)=(x-iH_{x})(x-(i+1)H_{x})(x-(i+2)H_{x}).$

Note that its first derivative
$Z^{\prime}(x)=3x^{2}-6(i+1)H_{x}x+(3i^{2}+6i+2)H_{x}^{2}$ has two zero-points
$x_{-}=\left(i+1-\frac{\sqrt{3}}{3}\right)H_{x}$ and
$x_{+}=\left(i+1+\frac{\sqrt{3}}{3}\right)H_{x}$. Thus, we have

$\displaystyle\max_{x\in(x^{H}_{i-1},x^{H}_{i})}\left|Z(x)\right|=\max\left\\{Z\left(\left(i-1\right)H_{x}\right),Z\left(iH_{x}\right)\right\\}=6H_{x}^{3},$
$\displaystyle\max_{x\in(x^{H}_{i},x^{H}_{i+1})}\left|Z(x)\right|=\max\left\\{Z\left(iH_{x}\right),Z\left(\left(i+1-\frac{\sqrt{3}}{3}\right)H_{x}\right),Z\left(\left(i+1\right)H_{x}\right)\right\\}=\frac{2\sqrt{3}}{9}H_{x}^{3},$
$\displaystyle\max_{x\in(x^{H}_{i+1},x^{H}_{i+2})}\left|Z(x)\right|=\max\left\\{Z\left(\left(i+1\right)H_{x}\right),Z\left(\left(i+1+\frac{\sqrt{3}}{3}\right)H_{x}\right),Z\left(\left(i+2\right)H_{x}\right)\right\\}=\frac{2\sqrt{3}}{9}H_{x}^{3},$

which, together with the definition of $\phi^{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^{\infty}$ norms for the piecewise bi-cubic Lagrange
interpolation operator $\Pi_{H}$.

###### Lemma 2.5

For any $w\in\mathcal{V}_{H}^{0}$, the following estimate holds

$\displaystyle\|\Pi_{H}w\|_{h}\leq C_{3}\|w\|_{H},$

where $C_{3}=4\left(\frac{3+27^{2}+(10+7\sqrt{7})^{2}}{27^{2}}\right)$.

Proof. Denote $\xi_{i,j}=\Pi_{H,x}w_{i,j}$, and then

$\displaystyle\|\Pi_{H}w\|^{2}_{h}=h_{x}h_{y}\sum^{N^{h}_{x}-1}_{i=1}\sum^{N^{h}_{y}-1}_{j=1}\left(\Pi_{H}w_{i,j}\right)^{2}=h_{x}\sum^{N^{h}_{x}-1}_{i=1}\left(h_{y}\sum^{N^{h}_{y}-1}_{j=1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}\right).$
(2.3)

For fixed $1\leq i\leq N^{h}_{x}-1$, we know

$\displaystyle
h_{y}\sum^{N^{h}_{y}-1}_{j=1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}$
$\displaystyle=h_{y}\sum^{N^{H}_{y}}_{k=1}\sum^{kM_{y}}_{j=(k-1)M_{y}+1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}$
(2.4)
$\displaystyle=h_{y}\sum^{M_{y}}_{j=1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}+h_{y}\sum^{N^{H}_{y}-1}_{k=2}\sum^{kM_{y}}_{j=(k-1)M_{y}+1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}+h_{y}\sum^{N^{h}_{y}-1}_{j=(N^{H}_{y}-1)M_{y}+1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}$
$\displaystyle=:I_{1}+I_{2}+I_{3}.$

We use the identity $\left(\sum^{4}_{i=1}a_{i}\right)^{2}\leq
4\sum^{4}_{i=1}a^{2}_{i}$ and the definition (2.2) of cubic Lagrange
interpolation operator $\Pi_{H,y}$ to obtain

$\displaystyle
I_{1}=h_{y}\sum^{M_{y}}_{j=1}\left(\sum^{3}_{s=0}\phi^{y}_{0,s}(y^{h}_{j})\xi_{i,s}\right)^{2}\leq
4h_{y}\sum^{M_{y}}_{j=1}\sum^{3}_{s=0}\left(\phi^{y}_{0,s}(y^{h}_{j})\right)^{2}\left(\xi_{i,s}\right)^{2},$
(2.5)

which, together with Lemma 2.4, yields

$\displaystyle I_{1}$ $\displaystyle\leq
4h_{y}\sum^{M_{y}}_{j=1}\left(\xi_{i,0}^{2}+\frac{(7\sqrt{7}+10)^{2}}{27^{2}}\xi_{i,1}^{2}+\frac{(7\sqrt{7}-10)^{2}}{27^{2}}\xi_{i,2}^{2}+\frac{3}{27^{2}}\xi_{i,3}^{2}\right)$
(2.6)
$\displaystyle=4H_{y}\left(\xi_{i,0}^{2}+\frac{(7\sqrt{7}+10)^{2}}{27^{2}}\xi_{i,1}^{2}+\frac{(7\sqrt{7}-10)^{2}}{27^{2}}\xi_{i,2}^{2}+\frac{3}{27^{2}}\xi_{i,3}^{2}\right).$

Analogous to the process (2.5)–(2.6), we can easily derive

$\displaystyle I_{2}\leq
4H_{y}\sum^{N^{H}_{y}-1}_{k=2}\left(\frac{3}{27^{2}}\xi_{i,k-2}^{2}+\xi_{i,k-1}^{2}+\xi_{i,k}^{2}+\frac{3}{27^{2}}\xi_{i,k+1}^{2}\right).$
(2.7) $\displaystyle I_{3}\leq
4H_{y}\left(\frac{3}{27^{2}}\xi_{i,N^{H}_{y}-3}^{2}+\frac{(7\sqrt{7}-10)^{2}}{27^{2}}\xi_{i,N^{H}_{y}-2}^{2}+\frac{(7\sqrt{7}+10)^{2}}{27^{2}}\xi_{i,N^{H}_{y}-1}^{2}+\xi_{i,N^{H}_{y}}^{2}\right).$
(2.8)

Inserting (2.6)–(2.8) into (2.4) gives us

$\displaystyle
h_{y}\sum^{N^{h}_{y}-1}_{j=1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}$
$\displaystyle\leq
4H_{y}\left(\frac{3+27^{2}+(7\sqrt{7}+10)^{2}}{27^{2}}\xi_{i,1}^{2}+\frac{3+2\times
27^{2}+(7\sqrt{7}-10)^{2}}{27^{2}}\xi_{i,2}^{2}\right.$ (2.9)
$\displaystyle\qquad\quad+\frac{9+2\times
27^{2}}{27^{2}}\xi_{i,3}^{2}+\frac{6+2\times
27^{2}}{27^{2}}\sum^{N^{H}_{y}-4}_{j=4}\xi_{i,j}^{2}+\frac{9+2\times
27^{2}}{27^{2}}\xi_{i,N^{H}_{y}-3}^{2}$
$\displaystyle\qquad\quad\left.+\frac{3+2\times
27^{2}+(7\sqrt{7}-10)^{2}}{27^{2}}\xi_{i,N^{H}_{y}-2}^{2}+\frac{3+27^{2}+(7\sqrt{7}+10)^{2}}{27^{2}}\xi_{i,N^{H}_{y}-1}^{2}\right)$
$\displaystyle\leq C_{3}H_{y}\sum^{N^{H}_{y}-1}_{j=1}\xi_{i,j}^{2},$

which further implies

$\displaystyle\|\Pi_{H}w\|^{2}_{h}=h_{x}\sum^{N^{h}_{x}-1}_{i=1}\left(h_{y}\sum^{N^{h}_{y}-1}_{j=1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}\right)\leq
C_{3}H_{y}\sum^{N^{H}_{y}-1}_{j=1}\left(h_{x}\sum^{N^{h}_{x}-1}_{i=1}\left(\Pi_{H,x}w_{i,j}\right)^{2}\right).$
(2.10)

For fixed $1\leq j\leq N^{H}_{y}-1$, we replace $\left\\{\xi,y,j\right\\}$
with $\left\\{w,x,i\right\\}$ in (2.4)–(2.10) to similarly obtain

$\displaystyle h_{x}\sum^{N^{h}_{x}-1}_{i=1}\left(\Pi_{H,x}w_{i,j}\right)^{2}$
$\displaystyle\leq C_{3}H_{x}\sum^{N^{H}_{x}-1}_{i=1}w_{i,j}^{2}.$ (2.11)

Consequently, we have

$\displaystyle\|\Pi_{H}w\|^{2}_{h}\leq
C_{3}^{2}H_{x}H_{y}\sum^{N^{H}_{x}-1}_{i=1}\sum^{N^{H}_{y}-1}_{j=1}w_{i,j}^{2}=C_{3}^{2}\|w\|^{2}_{H}.$
(2.12)

The proof is completed.

###### Lemma 2.6

For any $w\in\mathcal{V}_{H}^{0}$, the following estimate holds

$\displaystyle\|\Pi_{H}w\|_{h,\infty}\leq C_{4}\|w\|_{H,\infty},$

where $C_{4}=\left(\frac{54+2\sqrt{3}}{27}\right)^{2}$.

Proof. Suppose $\|\Pi_{H}w\|_{h,\infty}=\left|\Pi_{H}w_{i^{*},j^{*}}\right|$
where $iM_{x}\leq i^{*}\leq(i+1)M_{x}$ and $jM_{y}\leq j^{*}\leq(j+1)M_{y}$.
Denote $\xi_{i^{*},j}=\Pi_{H,x}w_{i^{*},j}$ for $0\leq j\leq N^{H}_{y}$, then
the triangle inequality and Lemma 2.4 give us

$\displaystyle\|\Pi_{H}w\|_{h,\infty}$
$\displaystyle=\left|\Pi_{H,y}\xi_{i^{*},j^{*}}\right|=\left|\sum^{3}_{s=0}\phi^{y}_{j,s}(y^{h}_{j^{*}})\xi_{i^{*},j-1+s}\right|$
(2.13)
$\displaystyle\leq\frac{\sqrt{3}}{27}\left|\xi_{i^{*},j-1}\right|+\left|\xi_{i^{*},j}\right|+\left|\xi_{i^{*},j+1}\right|+\frac{\sqrt{3}}{27}\left|\xi_{i^{*},j+2}\right|$
$\displaystyle\leq\frac{54+2\sqrt{3}}{27}\max_{1\leq k\leq
N^{H}_{y}-1}\left|\xi_{i^{*},k}\right|,\qquad\text{ if }\ M_{y}\leq
j^{*}\leq(N^{H}_{y}-1)M_{y},$

and

$\displaystyle\|\Pi_{H}w\|_{h,\infty}$
$\displaystyle\leq\frac{7\sqrt{7}+10}{27}\left|\xi_{i^{*},1}\right|+\frac{7\sqrt{7}-10}{27}\left|\xi_{i^{*},2}\right|+\frac{\sqrt{3}}{27}\left|\xi_{i^{*},3}\right|$
(2.14) $\displaystyle\leq\frac{14\sqrt{7}+\sqrt{3}}{27}\max_{1\leq k\leq
N^{H}_{y}-1}\left|\xi_{i^{*},k}\right|,\qquad\text{ if }\ 1\leq j^{*}<M_{y},$
$\displaystyle\|\Pi_{H}w\|_{h,\infty}\leq\frac{14\sqrt{7}+\sqrt{3}}{27}\max_{1\leq
k\leq N^{H}_{y}-1}\left|\xi_{i^{*},k}\right|,\qquad\text{ if }\
(N^{H}_{y}-1)M_{y}<j^{*}\leq N^{h}_{y}-1.$ (2.15)

Next, it is necessary to estimate $\max_{1\leq k\leq
N^{H}_{y}-1}\left|\xi_{i^{*},k}\right|$. For fixed index $1\leq k\leq
N^{H}_{y}-1$, let $\left|\xi_{i^{**},k}\right|=\max_{1\leq\ell\leq
N^{h}_{x}-1}\left|\xi_{\ell,k}\right|$ where $iM_{x}\leq
i^{**}\leq(i+1)M_{x}$. Thus, with a similar treatment to the above estimates,
we derive

$\displaystyle\left|\xi_{i^{*},k}\right|=\left|\Pi_{H,x}w_{i^{*},k}\right|\leq\left|\Pi_{H,x}w_{i^{**},k}\right|\leq\frac{54+2\sqrt{3}}{27}\max_{1\leq\ell\leq
N^{H}_{x}-1}\left|w_{\ell,k}\right|,$ (2.16)

which, together with (2.13)–(2.15), completes the proof.

##  3 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

$\displaystyle u\in C(0,T;H^{6}(\Omega))\cap C^{3}(0,T;L^{2}(\Omega)).$ (3.1)

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 $u\in B:=[m,M]$. Furthermore, for a fixed small
$\delta>0$, we denote $B_{\delta}:=[m-\delta,M+\delta]$. In this paper, we
only assume that the nonlinear term $f\in C^{2}(B_{\delta})$ holds locally on
$B_{\delta}$, that is,

$|f(w)|+|f^{\prime}(w)|+|f^{\prime\prime}(w)|\leq K_{f}\qquad\text{for }\ w\in
B_{\delta}.$ (3.2)

###  3.1 Nonlinear compact difference scheme

To construct variable temporal stepsize schemes, we consider a nonuniform
temporal mesh partition $0=t_{0}<t_{1}<\cdots<t_{N}=T$ with temporal stepsize
$\tau_{k}:=t_{k}-t_{k-1}$ for $1\leq k\leq N$. Denote the maximum temporal
stepsize $\tau:=\max_{1\leq k\leq N}\tau_{k}$ and adjacent temporal stepsize
ratio $r_{k}:=\tau_{k}/\tau_{k-1}$ $(k\geq 2)$. Set $w^{k}=w(t_{k})$ and
$\nabla_{\tau}w^{k}=w^{k}-w^{k-1}$, then the variable-step BDF2 formula is
defined by

$\mathcal{D}_{2}w^{n}:=\frac{1+2r_{n}}{\tau_{n}(1+r_{n})}\nabla_{\tau}w^{n}-\frac{r_{n}^{2}}{\tau_{n}(1+r_{n})}\nabla_{\tau}w^{n-1}\quad\text{for
}\ n\geq 2.$

In particular, when $n=1$, we use the BDF1 (i.e., backward Euler) formula
$\mathcal{D}_{1}w^{1}=\frac{1}{\tau_{1}}\nabla_{\tau}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

$\displaystyle\mathcal{D}_{2}w^{n}=\sum^{n}_{k=1}b^{(n)}_{n-k}\nabla_{\tau}w^{k}\quad\text{for
}\ n\geq 1,$ (3.3)

where the discrete convolution kernels $b^{(n)}_{n-k}$ are defined by
$b^{(1)}_{0}:=1/\tau_{1}$ and

$\displaystyle b^{(n)}_{0}:=\frac{1+2r_{n}}{\tau_{n}(1+r_{n})},\quad
b^{(n)}_{1}:=-\frac{r_{n}^{2}}{\tau_{n}(1+r_{n})},\quad\text{and}\quad
b^{(n)}_{j}:=0\ \ \text{for}\ 2\leq j\leq n.$ (3.4)

Next, we introduce the discrete orthogonal convolution (DOC) kernels of
$\left\\{b^{(n)}_{n-k}\right\\}$ by [30, 25, 9]

$\displaystyle\sum^{n}_{m=k}\theta^{(n)}_{n-m}b^{(m)}_{m-k}=\delta_{nk},\quad
1\leq k\leq n,\ 1\leq n\leq N,$ (3.5)

where $\delta_{nk}=1$ if $n=k$ and $\delta_{nk}=0$ if $n\neq k$. By exchanging
the summation order and using definition (3.5), it is easy to check that

$\displaystyle\sum_{m=1}^{n}\theta_{n-m}^{(n)}\mathcal{D}_{2}w^{m}=\sum_{k=1}^{n}\nabla_{\tau}w^{k}\sum_{m=k}^{n}\theta_{n-m}^{(n)}b_{m-k}^{(m)}=\nabla_{\tau}w^{n},\quad
1\leq n\leq N.$ (3.6)

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 $\mathcal{A}_{h}$ in space and
variable-step BDF method in time to get the following equation

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}U^{n}_{i,j}-c\Lambda_{h}U^{n}_{i,j}=\mathcal{A}_{h}f(U^{n}_{i,j})+\mathcal{A}_{h}g^{n}_{i,j}+R^{n}_{i,j},\quad
1\leq n\leq N,$ (3.7)

for $1\leq i\leq N^{h}_{x}-1$, $1\leq j\leq N^{h}_{y}-1$. Here
$R^{n}_{i,j}:=(R_{t}^{n})_{i,j}+(R_{s}^{n})_{i,j}$ is the local truncation
error at point $(x^{h}_{i},y^{h}_{i},t_{n})$, where

$(R_{t}^{n})_{i,j}=\mathcal{A}_{h}\left(u_{t}(x^{h}_{i},y^{h}_{j},t_{n})-\mathcal{D}_{2}U^{n}_{i,j}\right)\quad\text{
and }\quad(R_{s}^{n})_{i,j}=-c\mathcal{A}_{h}\Delta
u(x^{h}_{i},y^{h}_{i},t_{n})+c\Lambda_{h}U^{n}_{i,j}.$

Under the regularity condition (3.1), by using the Taylor expansion and well-
known Bramble-Hilbert Lemma, it is easy to see

$\displaystyle\|R_{t}^{1}\|_{h}\leq C_{5}\tau_{1},\quad\|R_{t}^{n}\|_{h}\leq
C_{6}\tau_{n}\tau\quad\text{for}\ 2\leq n\leq N,$ (3.8)
$\displaystyle\|R_{s}^{n}\|_{h}\leq C_{7}h^{4}\quad\text{for}\ 1\leq n\leq N.$
(3.9)

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

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}u^{n}_{i,j}-c\Lambda_{h}u^{n}_{i,j}=\mathcal{A}_{h}f(u^{n}_{i,j})+\mathcal{A}_{h}g^{n}_{i,j},\quad(i,j)\in\omega_{h},~{}1\leq
n\leq N.$ (3.10)

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
$\bar{u}^{n}=\\{\bar{u}^{n}_{i,j}\\}$ satisfying the following scheme

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}\bar{u}^{n}_{i,j}-c\Lambda_{h}\bar{u}^{n}_{i,j}=\mathcal{A}_{h}\bar{f}(\bar{u}^{n}_{i,j})+\mathcal{A}_{h}g^{n}_{i,j},\quad(i,j)\in\omega_{h},~{}1\leq
n\leq N,$ (3.11)

subject to the same initial and boundary conditions (1.2)–(1.3). Here
$\bar{f}(w)$ represents a cut-off function of $f$ on $B_{\delta}$ such that

$\begin{array}[]{l}\bar{f}(w):=\left\\{\begin{split}&f(m-\delta),\qquad\
w<m-\delta,\\\ &f(w),\qquad\qquad w\in B_{\delta},\\\ &f(M+\delta),\qquad\
w>M+\delta.\end{split}\right.\end{array}$ (3.12)

It is obvious that $\bar{f}(w)$ is globally Lipschitz continuous on
$\mathbb{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.

######  Lemma 3.7 ([9, 25])

Assume that the adjacent temporal stepsize ratios $r_{k}$ satisfy
$0<r_{k}<4.8645$. The DOC kernels $\\{\theta^{(n)}_{n-k}\\}$ defined in (3.5)
are positive semi-definite, i.e., for any real sequence
$\\{w_{k}\\}^{n}_{k=1}$, it holds that

$\displaystyle\sum_{k=1}^{n}w_{k}\sum_{m=1}^{k}\theta_{k-m}^{(k)}w_{m}\geq
0\quad{\rm for}\ n\geq 1.$

######  Lemma 3.8 ([30, 9])

The DOC kernels $\\{\theta^{(n)}_{n-k}\\}$ defined in (3.5) have the following
properties

$\displaystyle\theta_{n-m}^{(n)}>0\quad{\rm for}\ 1\leq m\leq
n,\quad\sum_{m=1}^{n}\theta_{n-m}^{(n)}=\tau_{n}\quad{\rm for}\ n\geq 1,$
$\displaystyle\sum^{n}_{m=1}\sum^{n}_{k=m}\theta^{(k)}_{k-m}=t_{n},\quad{\rm
and}\quad\sum^{n}_{k=m}\theta^{(k)}_{k-m}\leq 2\tau\quad{\rm for}\ 1\leq m\leq
n.$

######  Lemma 3.9 ([30])

Let $\lambda\geq 0$, and the time sequences $\\{\xi^{k}\\}_{k=0}^{N}$ and
$\\{w^{k}\\}_{k=1}^{N}$ be nonnegative. If

$\displaystyle
w^{n}\leq\lambda\sum^{n-1}_{k=1}\tau_{k}w^{k}+\sum^{n}_{k=0}\xi^{k}\quad{\rm
for}\ 1\leq n\leq N,$ (3.13)

then it holds that

$\displaystyle w^{n}\leq\exp(\lambda t_{n-1})\sum^{n}_{k=0}\xi^{k}\quad{\rm
for}\ 1\leq n\leq N.$ (3.14)

###  3.2  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.10

The auxiliary nonlinear compact difference scheme (3.11) is solvable if the
maximum temporal stepsize satisfies $\tau\leq\frac{1}{4K_{f}}$.

Proof. Note that (3.11) can be equivalently rewritten as

$\displaystyle
b^{(n)}_{0}\mathcal{A}_{h}\bar{u}^{n}_{i,j}-G^{n}_{i,j}-c\Lambda_{h}\bar{u}^{n}_{i,j}-\mathcal{A}_{h}\bar{f}(\bar{u}^{n}_{i,j})-\mathcal{A}_{h}g^{n}_{i,j}=0,\quad(i,j)\in\omega_{h},$
(3.15)

where $G^{n}_{i,j}$ is defined as

$\displaystyle
G^{n}_{i,j}:=(b^{(n)}_{0}-b^{(n)}_{1})\mathcal{A}_{h}\bar{u}^{n-1}_{i,j}+b^{(n)}_{1}\mathcal{A}_{h}\bar{u}^{n-2}_{i,j}.$

Denote the mapping

$[T\bar{u}]^{n}_{i,j}:=b^{(n)}_{0}\mathcal{A}_{h}\bar{u}^{n}_{i,j}-G^{n}_{i,j}-c\Lambda_{h}\bar{u}^{n}_{i,j}-\mathcal{A}_{h}\bar{f}(\bar{u}^{n}_{i,j})-\mathcal{A}_{h}g^{n}_{i,j},\quad(i,j)\in\omega_{h}.$

Then, taking the inner product of $T\bar{u}^{n}$ with $\bar{u}^{n}$ in the
sense of $(\cdot,\cdot)_{h}$ gives us

$\displaystyle\left(T\bar{u}^{n},\bar{u}^{n}\right)_{h}$
$\displaystyle=b^{(n)}_{0}\left(\mathcal{A}_{h}\bar{u}^{n},\bar{u}^{n}\right)_{h}-\left(G^{n},\bar{u}^{n}\right)_{h}-c\left(\Lambda_{h}\bar{u}^{n},\bar{u}^{n}\right)_{h}$
(3.16)
$\displaystyle\quad-\left(\mathcal{A}_{h}\bar{f}(\bar{u}^{n}),\bar{u}^{n}\right)_{h}-\left(\mathcal{A}_{h}g^{n},\bar{u}^{n}\right)_{h}:=\sum^{5}_{k=1}S_{k}.$

Next, we estimate (3.16) term-by-term. For $S_{2}$, utilizing Cauchy-Schwarz
inequality yields

$\displaystyle
S_{2}\geq-\left((b^{(n)}_{0}-b^{(n)}_{1})\|\bar{u}^{n-1}\|_{\mathcal{A},h}-b^{(n)}_{1}\|\bar{u}^{n-2}\|_{\mathcal{A},h}\right)\|\bar{u}^{n}\|_{\mathcal{A},h}.$

Moreover, by summation by parts and homogeneous boundary conditions, and
noting that $\mathcal{A}_{h,x}$ and $\mathcal{A}_{h,y}$ are self-adjoint and
positive definite operators, we have

$\displaystyle
S_{3}=\left(\mathcal{A}_{h,y}d_{h,x}\bar{u}^{n},d_{h,x}\bar{u}^{n}\right)_{h,x}+\left(\mathcal{A}_{h,x}d_{h,y}\bar{u}^{n},d_{h,y}\bar{u}^{n}\right)_{h,y}\geq
0.$

Besides, we use Cauchy-Schwarz inequality to get

$\displaystyle
S_{4}+S_{5}=-\left(\mathcal{A}_{h}\left(\bar{f}(\bar{u}^{n})+g^{n}\right),\bar{u}^{n}\right)_{h}\geq-\interleave\bar{f}(\bar{u}^{n})+g^{n}\interleave_{h}\|\bar{u}^{n}\|_{h},$

which further, together with the globally Lipschitz continuous of $\bar{f}$
and Lemma 2.1, leads to

$\displaystyle S_{4}+S_{5}$ $\displaystyle\geq-
K_{f}\left(\|\bar{u}^{n}\|_{h}+\|u^{0}\|_{h}\right)\|\bar{u}^{n}\|_{h}-\interleave\bar{f}(u^{0})+g^{n}\interleave_{h}\|\bar{u}^{n}\|_{h}$
$\displaystyle\geq-3K_{f}\|\bar{u}^{n}\|^{2}_{\mathcal{A},h}-\sqrt{3}\left(K_{f}\|u^{0}\|_{h}+\interleave\bar{f}(u^{0})+g^{n}\interleave_{h}\right)\|\bar{u}^{n}\|_{\mathcal{A},h}.$

Denote

$\displaystyle
Q:=(b^{(n)}_{0}-b^{(n)}_{1})\|\bar{u}^{n-1}\|_{\mathcal{A},h}-b^{(n)}_{1}\|\bar{u}^{n-2}\|_{\mathcal{A},h}+\sqrt{3}\left(K_{f}\|u^{0}\|_{h}+\interleave\bar{f}(u^{0})+g^{n}\interleave_{h}\right).$

Then inserting the above estimates into (3.16) gives

$\displaystyle\left(T\bar{u}^{n},\bar{u}^{n}\right)_{h}\geq\left(b^{(n)}_{0}-3K_{f}\right)\|\bar{u}^{n}\|^{2}_{\mathcal{A},h}-Q\|\bar{u}^{n}\|_{\mathcal{A},h}.$

Thus, if $\tau_{n}\leq\frac{1}{4K_{f}}$, we conclude that

$\displaystyle\left(T\bar{u}^{n},\bar{u}^{n}\right)_{h}\geq\left(K_{f}\|\bar{u}^{n}\|_{\mathcal{A},h}-Q\right)\|\bar{u}^{n}\|_{\mathcal{A},h}\geq
0,$

when $\|\bar{u}^{n}\|_{\mathcal{A},h}=Q/K_{f}$. Consequently, Browder’s fixed
point theorem shows there exists a $\bar{u}^{n}\in\mathcal{V}^{0}_{h}$ such
that $T\bar{u}^{n}=0$ and
$\|\bar{u}^{n}\|_{h}\leq\sqrt{3}\|\bar{u}^{n}\|_{\mathcal{A},h}\leq\sqrt{3}Q/K_{f}$,
which implies the solvability of scheme (3.11).

###### Theorem 3.11

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
$\tau\leq\frac{1}{4\sqrt{3}K_{f}}$.

Proof. The argument is by contradiction. Suppose that there are two solutions
$\bar{u}_{1}^{n}$ and $\bar{u}_{2}^{n}$ satisfying (3.11) and the same initial
and boundary conditions (1.2)–(1.3), i.e.,

$\displaystyle\mathcal{D}_{2}[\mathcal{A}_{h}\bar{u}_{1}]^{n}_{i,j}-c[\Lambda_{h}\bar{u}_{1}]^{n}_{i,j}=\mathcal{A}_{h}\bar{f}([\bar{u}_{1}]^{n}_{i,j})+\mathcal{A}_{h}g^{n}_{i,j},$
$\displaystyle\mathcal{D}_{2}[\mathcal{A}_{h}\bar{u}_{2}]^{n}_{i,j}-c[\Lambda_{h}\bar{u}_{2}]^{n}_{i,j}=\mathcal{A}_{h}\bar{f}([\bar{u}_{2}]^{n}_{i,j})+\mathcal{A}_{h}g^{n}_{i,j},$

for $1\leq i\leq N^{h}_{x}-1$, $1\leq j\leq N^{h}_{y}-1$ and $1\leq n\leq N$.
Denote $v^{n}:=\bar{u}_{1}^{n}-\bar{u}_{2}^{n}$ and subtract these two
equations to get

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}v^{n}_{i,j}-c\Lambda_{h}v^{n}_{i,j}=\mathcal{A}_{h}\bar{f}([\bar{u}_{1}]^{n}_{i,j})-\mathcal{A}_{h}\bar{f}([\bar{u}_{2}]^{n}_{i,j}).$
(3.17)

In (3.17), we let $n=m$, and then, multiplying it by the DOC kernels
$\theta^{(k)}_{k-m}$ and summing $m$ from $1$ to $k$ yields

$\displaystyle\nabla_{\tau}\mathcal{A}_{h}v^{k}_{i,j}-c\sum^{k}_{m=1}\theta^{(k)}_{k-m}\Lambda_{h}v^{m}_{i,j}=\sum^{k}_{m=1}\theta^{(k)}_{k-m}\mathcal{A}_{h}\left(\bar{f}([\bar{u}_{1}]^{m}_{i,j})-\bar{f}([\bar{u}_{2}]^{m}_{i,j})\right),$
(3.18)

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

$\displaystyle\|v^{n}\|_{\mathcal{A},h}^{2}-2c\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}(\Lambda_{h}v^{m},v^{k})_{h}$
$\displaystyle\leq
2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left(\mathcal{A}_{h}\bar{f}(\bar{u}_{1}^{n})-\mathcal{A}_{h}\bar{f}(\bar{u}_{2}^{n}),v^{k}\right)_{h},$
(3.19)

where $v^{0}=0$ has been used.

Due to $\mathcal{A}_{h,x}$ and $\mathcal{A}_{h,y}$ are self-adjoint and
positive definite operators, there exist $\eta_{x}$ and $\eta_{y}$ such that
$\mathcal{A}_{h,x}=\eta^{2}_{x}$ and $\mathcal{A}_{h,y}=\eta^{2}_{y}$. Thus,
by summation by parts and homogeneous boundary conditions, the second left-
hand side term of (3.19) could be rewritten as

$\displaystyle\quad-2c\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}(\Lambda_{h}v^{m},v^{k})_{h}$
(3.20)
$\displaystyle=2c\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left((\mathcal{A}_{h,x}d_{h,y}v^{m},d_{h,y}v^{k})_{h}+(\mathcal{A}_{h,y}d_{h,x}v^{m},d_{h,x}v^{k})_{h}\right)$
$\displaystyle=2c\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left((\eta_{x}d_{h,y}v^{m},\eta_{x}d_{h,y}v^{k})_{h}+(\eta_{y}d_{h,x}v^{m},\eta_{y}d_{h,x}v^{k})_{h}\right)\geq
0,$

where the positive semi-definiteness of the DOC kernels $\theta^{(k)}_{k-m}$
(cf. Lemma 3.7) has been used in the last inequality.

Using Cauchy-Schwarz inequality and Lemma 2.1 we see

$\displaystyle
2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left(\mathcal{A}_{h}\bar{f}(\bar{u}_{1}^{n})-\mathcal{A}_{h}\bar{f}(\bar{u}_{2}^{n}),v^{k}\right)_{h}$
$\displaystyle\leq
2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|\bar{f}(\bar{u}_{1}^{n})-\bar{f}(\bar{u}_{2}^{n})\|_{\mathcal{A},h}\|v^{k}\|_{\mathcal{A},h}$
$\displaystyle\leq
2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|\bar{f}(\bar{u}_{1}^{n})-\bar{f}(\bar{u}_{2}^{n})\|_{h}\|v^{k}\|_{\mathcal{A},h},$

which, together with the global Lipschitz continuous property of $\bar{f}$,
gives

$\displaystyle
2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left(\mathcal{A}_{h}\bar{f}(\bar{u}_{1}^{n})-\mathcal{A}_{h}\bar{f}(\bar{u}_{2}^{n}),v^{k}\right)$
$\displaystyle\leq
2K_{f}\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|v^{m}\|_{h}\|v^{k}\|_{\mathcal{A},h}$
(3.21) $\displaystyle\leq
2\sqrt{3}K_{f}\sum^{n}_{k=1}\|v^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|v^{m}\|_{\mathcal{A},h}.$

Now, inserting (3.20)–(3.21) into (3.19), we obtain

$\displaystyle\|v^{n}\|_{\mathcal{A},h}^{2}\leq
2\sqrt{3}K_{f}\sum^{n}_{k=1}\|v^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|v^{m}\|_{\mathcal{A},h}.$

Choosing $n^{*}$ $(1\leq n^{*}\leq n)$ such that
$\|v^{n^{*}}\|_{\mathcal{A},h}=\max_{1\leq k\leq n}\|v^{n}\|_{\mathcal{A},h}$.
Then, the above inequality yields

$\displaystyle\|v^{n^{*}}\|_{\mathcal{A},h}^{2}\leq
2\sqrt{3}K_{f}\sum^{n^{*}}_{k=1}\|v^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|v^{n^{*}}\|_{\mathcal{A},h},$

which, eliminating $\|v^{n^{*}}\|_{\mathcal{A},h}$ from both sides and using
Lemma 3.8, implies

$\displaystyle\|v^{n}\|_{\mathcal{A},h}\leq\|v^{n^{*}}\|_{\mathcal{A},h}\leq
2\sqrt{3}K_{f}\sum^{n^{*}}_{k=1}\|v^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\leq
2\sqrt{3}K_{f}\sum^{n}_{k=1}\tau_{k}\|v^{k}\|_{\mathcal{A},h}.$ (3.22)

Then, for $\tau_{n}\leq 1/(4\sqrt{3}K_{f})$, an application of the discrete
Gr${\rm\ddot{o}}$nwall inequality in Lemma 3.9 gives

$\displaystyle\|v^{n}\|_{\mathcal{A},h}\leq
4\sqrt{3}K_{f}\sum^{n-1}_{k=1}\tau_{k}\|v^{k}\|_{\mathcal{A},h}\leq 0,$ (3.23)

which proves the uniqueness of the solutions to scheme (3.11).

###  3.3  Error analysis

Denote $\bar{e}^{n}:=U^{n}-\bar{u}^{n}$ and $e^{n}:=U^{n}-u^{n}$ for $0\leq
n\leq N$. It is easy to see that $\bar{e}^{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).

###### Theorem 3.12

Assume that the solution of (1.1)–(1.3) satisfy the regularity assumption
condition (3.1) and the adjacent temporal stepsize ratios $r_{k}$ satisfy
$0<r_{k}<4.8645$. If the maximum temporal stepsize
$\tau\leq\frac{1}{4\sqrt{3}K_{f}}$, the following estimate holds for the
auxiliary nonlinear compact difference scheme (3.11) that

$\|U^{n}-\bar{u}^{n}\|_{h}\leq C_{8}\left(\tau^{2}+h^{4}\right)\quad{\rm for}\
1\leq n\leq N,$

where
$C_{8}:=12\exp\big{(}4\sqrt{3}K_{f}T\big{)}\max\\{2C_{5}+C_{6}T,3C_{7}T\\}$.

Proof. Subtracting (3.11) from (3.7), we can get the following error equation

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}\bar{e}^{n}_{i,j}-c\Lambda_{h}\bar{e}^{n}_{i,j}=\mathcal{A}_{h}\bar{f}(U^{n}_{i,j})-\mathcal{A}_{h}\bar{f}(\bar{u}^{n}_{i,j})+R^{n}_{i,j},\quad
1\leq n\leq N.$ (3.24)

for $1\leq i\leq N^{h}_{x}-1$ and $1\leq j\leq N^{h}_{y}-1$. Analogous to the
proof of (3.17)–(3.18), it follows from the orthogonality (3.5) of the DOC
kernels $\\{\theta^{(k)}_{k-m}\\}$ that

$\displaystyle\nabla_{\tau}\mathcal{A}_{h}\bar{e}^{k}_{i,j}-c\sum^{k}_{m=1}\theta^{(k)}_{k-m}\Lambda_{h}\bar{e}^{m}_{i,j}=\sum^{k}_{m=1}\theta^{(k)}_{k-m}\mathcal{A}_{h}\left(\bar{f}(U^{m}_{i,j})-\bar{f}(\bar{u}^{m}_{i,j})\right)+\sum^{k}_{m=1}\theta^{(k)}_{k-m}R^{m}_{i,j},$
(3.25)

which, by taking the inner product with $2\bar{e}^{k}$ and summing the
resulting equality from $k=1$ to $n$, gives

$\displaystyle\|\bar{e}^{n}\|_{\mathcal{A},h}^{2}-2c\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}(\Lambda_{h}\bar{e}^{m},\bar{e}^{k})_{h}$
$\displaystyle\leq
2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left(\mathcal{A}_{h}f(U^{m})-\mathcal{A}_{h}f(\bar{u}^{m}),\bar{e}^{k}\right)_{h}$
(3.26)
$\displaystyle\qquad+2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left(R^{m},\bar{e}^{k}\right)_{h}.$

With a similar treatment to (3.20)–(3.21), the above inequality yields

$\displaystyle\|\bar{e}^{n}\|_{\mathcal{A},h}^{2}\leq
2\sqrt{3}K_{f}\sum^{n}_{k=1}\|\bar{e}^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|\bar{e}^{m}\|_{\mathcal{A},h}+2\sqrt{3}\sum^{n}_{k=1}\|\bar{e}^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|R^{m}\|_{h}.$
(3.27)

Similarly, if $\|\bar{e}^{n^{*}}\|_{\mathcal{A},h}=\max_{1\leq k\leq
n}\|\bar{e}^{n}\|_{\mathcal{A},h}$, using Lemma 3.8 we have

$\displaystyle\|\bar{e}^{n}\|_{\mathcal{A},h}\leq\|\bar{e}^{n^{*}}\|_{\mathcal{A},h}$
$\displaystyle\leq
2\sqrt{3}K_{f}\sum^{n^{*}}_{k=1}\|\bar{e}^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}+2\sqrt{3}\sum^{n^{*}}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|R^{m}\|_{h}$
(3.28) $\displaystyle\leq
2\sqrt{3}K_{f}\sum^{n}_{k=1}\tau_{k}\|\bar{e}^{k}\|_{\mathcal{A},h}+2\sqrt{3}\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|R^{m}\|_{h}.$

For the last term of (3.28), we exchange the order of summation and utilize
Lemma 3.8 to obtain

$\displaystyle\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|R^{m}\|_{h}$
$\displaystyle=\sum^{n}_{m=1}\|R^{m}\|_{h}\sum^{n}_{k=m}\theta^{(k)}_{k-m}=\|R^{1}\|_{h}\sum^{n}_{k=1}\theta^{(k)}_{k-1}+\sum^{n}_{m=2}\|R^{m}\|_{h}\sum^{n}_{k=m}\theta^{(k)}_{k-m}$
$\displaystyle\leq 2\tau\|R^{1}\|_{h}+\max_{2\leq m\leq
n}\|R^{m}\|_{h}\sum^{n}_{m=2}\sum^{n}_{k=m}\theta^{(k)}_{k-m}\leq
2\tau\|R^{1}\|_{h}+t_{n}\max_{2\leq m\leq n}\|R^{m}\|_{h},$

which, together with the triangle inequality and (3.8)–(3.9), gives

$\displaystyle\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|R^{m}\|_{h}\leq\big{(}2C_{5}+C_{6}t_{n}\big{)}\
\tau^{2}+C_{7}\ \big{(}2\tau+t_{n}\big{)}h^{4}.$ (3.29)

Now, insert the above inequality into (3.28), for $\tau\leq
1/(4\sqrt{3}K_{f})$, we have

$\displaystyle\|\bar{e}^{n}\|_{\mathcal{A},h}\leq
4\sqrt{3}K_{f}\sum^{n-1}_{k=1}\tau_{k}\|\bar{e}^{k}\|_{\mathcal{A},h}+(8\sqrt{3}C_{5}+4\sqrt{3}C_{6}T)\tau^{2}+12\sqrt{3}C_{7}Th^{4}.$
(3.30)

Then, applying the discrete Grönwall inequality to (3.30) and using Lemma 2.1
yields

$\displaystyle\|\bar{e}^{n}\|_{h}\leq\sqrt{3}\|\bar{e}^{n}\|_{\mathcal{A},h}\leq
12\exp(4\sqrt{3}K_{f}T)\left((2C_{5}+C_{6}T)\tau^{2}+3C_{7}Th^{4}\right)\leq
C_{8}\left(\tau^{2}+h^{4}\right),$

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.12 and Lemma 2.2, we see

$\displaystyle\|\bar{e}^{n}\|_{h,\infty}\leq
C_{0}h^{-1}\|\bar{e}^{n}\|_{h}\leq C_{0}C_{8}h^{-1}(\tau^{2}+h^{4}).$

Consequently, if $\tau=o(h^{\frac{1}{2}})$, and for $\tau$, $h$ sufficiently
small, it holds that $\bar{u}^{n}_{i,j}\in B_{\delta}$, which further implies
$\bar{f}(\bar{u}^{n}_{i,j})=f(\bar{u}^{n}_{i,j})$ and thus in this context
$\bar{u}^{n}_{i,j}\equiv u^{n}_{i,j}$. We summarize these conclusions in the
following theorem.

###### Theorem 3.13

Under the conditions in Theorem 3.12 and if the maximum temporal stepsize
satisfies $\tau=o(h^{\frac{1}{2}})$, the nonlinear compact difference scheme
(3.10) admits a unique solution satisfying

$\displaystyle\|U^{n}-u^{n}\|_{h}\leq
C_{8}\left(\tau^{2}+h^{4}\right)\quad{\rm for}\ 1\leq n\leq N.$

##  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

$\mathcal{D}_{2}\mathcal{A}_{H}u^{n}_{H,i,j}-c\Lambda_{H}u^{n}_{H,i,j}=\mathcal{A}_{H}f(u^{n}_{H,i,j})+\mathcal{A}_{H}g^{n}_{i,j},\quad(i,j)\in\omega_{H},$
(4.1)

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

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}u^{n}_{h,i,j}-c\Lambda_{h}u^{n}_{h,i,j}=\mathcal{A}_{h}F^{n}_{i,j}+\mathcal{A}_{h}g^{n}_{i,j},\quad(i,j)\in\omega_{h},$
(4.2)

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

$\displaystyle
F^{n}_{i,j}:=f\left(\left[\Pi_{H}u_{H}^{n}\right]_{i,j}\right)+f^{\prime}\left(\left[\Pi_{H}u_{H}^{n}\right]_{i,j}\right)\left(u^{n}_{h,i,j}-\left[\Pi_{H}u_{H}^{n}\right]_{i,j}\right).$
(4.3)

On the coarse grid, denote $e_{H}^{n}=U^{n}-u_{H}^{n}$ for $0\leq n\leq N$.
Analogous to Theorem 3.13, we can immediately reach the following conclusions.

###### Theorem 4.14

Under the conditions in Theorem 3.12 and if the maximum temporal stepsize
satisfies $\tau=o(H^{\frac{1}{2}})$, the nonlinear compact difference scheme
(4.1) defined on the coarse grid admits a unique solution satisfying

$\displaystyle\|U^{n}-u_{H}^{n}\|_{H}\leq
C_{8}\left(\tau^{2}+H^{4}\right)\quad{\rm for}\ 1\leq n\leq N.$

Based on this theorem and Lemmas 2.5–2.6, the following corollary can be
derived.

###### Corollary 4.15

Under the conditions in Theorem 4.14, the numerical solution $u^{n}_{H}$ of
the nonlinear compact difference scheme (4.1) satisfies

$\displaystyle\|U^{n}-\Pi_{H}u^{n}_{H}\|_{h}\leq
C_{9}\left(\tau^{2}+H^{4}\right)\quad{\rm for}\ 1\leq n\leq N,$ (4.4)
$\displaystyle\|U^{n}-\Pi_{H}u^{n}_{H}\|_{h,\infty}\leq
C_{10}\left(H^{-1}\tau^{2}+H^{3}\right)\quad{\rm for}\ 1\leq n\leq N,$ (4.5)

and consequently the interpolation solution $\Pi_{H}u_{H}^{n}$ is bounded on
the fine grid and satisfies

$\Pi_{H}u_{H}^{n}\in B_{\delta},$ (4.6)

where $C_{9}=C_{2}+C_{3}C_{8}$, $C_{10}:=C_{1}+C_{0}C_{4}C_{8}$.

Proof. We perform the splitting

$\|U^{n}-\Pi_{H}u^{n}_{H}\|_{h}\leq\|U^{n}-\Pi_{H}U^{n}\|_{h}+\|\Pi_{H}U^{n}-\Pi_{H}u^{n}_{H}\|_{h},$

where the first right-hand side term could be bounded as
$\|U^{n}-\Pi_{H}U^{n}\|_{h}\leq C_{2}H^{4}$ via Lemma 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.14 give us

$\|\Pi_{H}U^{n}-\Pi_{H}u^{n}_{H}\|_{h}=\|\Pi_{H}e^{n}_{H}\|_{h}\leq
C_{3}\|e^{n}_{H}\|_{H}\leq C_{3}C_{8}(\tau^{2}+H^{4}),$

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.14, we
can conclude that

$\displaystyle\|U^{n}-\Pi_{H}u^{n}_{H}\|_{h,\infty}$
$\displaystyle\leq\|U^{n}-\Pi_{H}U^{n}\|_{h,\infty}+\|\Pi_{H}U^{n}-\Pi_{H}u^{n}_{H}\|_{h,\infty}$
$\displaystyle\leq C_{1}H^{4}+C_{4}\|e^{n}_{H}\|_{H,\infty}\leq
C_{1}H^{4}+C_{0}C_{4}H^{-1}\|e^{n}_{H}\|_{H}$
$\displaystyle\leq(C_{1}+C_{0}C_{4}C_{8})(H^{-1}\tau^{2}+H^{3}).$

Furthermore, an application of the triangle inequality yields

$\|\Pi_{H}u_{H}^{n}\|_{h,\infty}\leq\|U^{n}\|_{h,\infty}+\|U^{n}-\Pi_{H}U^{n}_{H}\|_{h,\infty}\leq\|U^{n}\|_{h,\infty}+C_{10}(H^{-1}\tau^{2}+H^{3}),$

which completes the proof of (4.6).

###### Remark 4.16

Corollary 4.15 declares that the interpolation solution $\Pi_{H}u^{n}_{H}$ is
bounded under the discrete $L^{\infty}$ norm on the fine grid, if the maximum
temporal stepsize satisfies $\tau=o(H^{\frac{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 $\Pi_{H}$ in the sense of $L^{\infty}$
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

$\displaystyle\|\Pi_{H}U^{n}-\Pi_{H}u^{n}_{H}\|_{h,\infty}\leq
C_{0}h^{-1}\|\Pi_{H}e^{n}_{H}\|_{h}\leq
C_{0}C_{3}C_{8}h^{-1}(\tau^{2}+H^{4}),$

which leads to a much worse restriction $\tau=o(h^{\frac{1}{2}})$. On the
other hand, in practical computation the condition $\tau=o(H^{\frac{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 $\tau=\mathcal{O}(H^{2})$ suggested
by Theorem 4.14 can naturally satisfy the restrictive condition.

At last, we shall give an error estimate for $e_{h}^{n}=U^{n}-u_{h}^{n}$ on
the fine grid for the linearized compact difference scheme (4.2).

###### Theorem 4.17

Under the conditions in Theorem 4.14, the following error estimate holds for
the solution of the two-grid compact difference scheme (4.1)–(4.2)

$\|U^{n}-u_{h}^{n}\|_{h}\leq C_{11}\left(\tau^{2}+h^{4}+H^{7}\right)\quad{\rm
for}\ 1\leq n\leq N,$

where
$C_{11}:=2\sqrt{3}\exp\big{(}4\sqrt{3}K_{f}T\big{)}\max\\{4\sqrt{3}C_{5}+2\sqrt{3}C_{6}T+C_{9}C_{10}K_{f}T,6\sqrt{3}C_{7}T\\}$.

Proof. For the linearized scheme (4.2) on the fine grid, we can get a very
similar error equation

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}e^{n}_{h,i,j}-c\Lambda_{h}e^{n}_{h,i,j}=\mathcal{A}_{h}f(U^{n}_{i,j})-\mathcal{A}_{h}F^{n}_{i,j}+R^{n}_{i,j},\quad(i,j)\in\omega_{h},$
(4.7)

which, together with a similar treatment to (3.24)–(3.26), leads to

$\displaystyle\|e_{h}^{n}\|_{\mathcal{A},h}^{2}-2c\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}(\Lambda_{h}e_{h}^{m},e_{h}^{k})$
$\displaystyle\leq
2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left(\mathcal{A}_{h}f(U^{m})-\mathcal{A}_{h}F^{m},e_{h}^{k}\right)$
(4.8)
$\displaystyle\quad+2\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\left(R^{m},e_{h}^{k}\right).$

Compared with the estimate (3.26) in Theorem 3.12, the only difference lies in
the treatment of the nonlinear term. We apply Taylor expansion of $f(U^{m})$
at $\Pi_{H}u_{H}^{m}$ to obtain

$\displaystyle f(U_{i,j}^{m})$
$\displaystyle=f\left(\left[\Pi_{H}u_{H}^{m}\right]_{i,j}\right)+f^{\prime}\left(\left[\Pi_{H}u_{H}^{m}\right]_{i,j}\right)\left(U^{m}_{i,j}-\left[\Pi_{H}u_{H}^{m}\right]_{i,j}\right)+\frac{1}{2}f^{\prime\prime}(\mu_{i,j}^{m})\left(U^{m}_{i,j}-\left[\Pi_{H}u_{H}\right]^{m}_{i,j}\right)^{2},$

for some constant $\mu^{m}_{i,j}$ between $U^{m}_{i,j}$ and
$[\Pi_{H}u_{H}^{m}]_{i,j}$. Then, subtract $F^{m}_{i,j}$ in (4.3) from this
equation, we have

$\displaystyle
f(U_{i,j}^{m})-F^{m}_{i,j}=f^{\prime}\left([\Pi_{H}u_{H}^{m}]_{i,j}\right)\left(U^{m}_{i,j}-u_{h,i,j}^{m}\right)+\frac{1}{2}f^{\prime\prime}(\mu_{i,j}^{m})\left(U^{m}_{i,j}-\left[\Pi_{H}u_{H}^{m}\right]_{i,j}\right)^{2}.$

By Corollary 4.15, if $\tau=o(H^{\frac{1}{2}})$, we have $\Pi_{H}u_{H}^{m}\in
B_{\delta}$ and $\mu^{m}\in B_{\delta}$. Thus, assumption (3.2) implies that
$\big{|}f^{\prime}([\Pi_{H}u_{H}^{m}]_{i,j})\big{|}+\big{|}f^{\prime\prime}(\mu_{i,j}^{m})\big{|}\leq
K_{f}$. Furthermore, we apply Lemma 2.1, Cauchy-Schwarz inequality and
Corollary 4.15 to derive

$\displaystyle\quad\left(\mathcal{A}_{h}f(U^{m})-\mathcal{A}_{h}F^{m},e_{h}^{k}\right)\leq\|f(U^{m})-F^{m}\|_{h}\|e_{h}^{k}\|_{\mathcal{A},h}$
(4.9) $\displaystyle\leq
K_{f}\|e_{h}^{m}\|_{h}\|e_{h}^{k}\|_{\mathcal{A},h}+\frac{1}{2}K_{f}\|U^{m}-\Pi_{H}u_{H}^{m}\|_{h,\infty}\|U^{m}-\Pi_{H}u_{H}^{m}\|_{h}\|e_{h}^{k}\|_{\mathcal{A},h}$
$\displaystyle\leq\sqrt{3}K_{f}\|e_{h}^{m}\|_{\mathcal{A},h}\|e_{h}^{k}\|_{\mathcal{A},h}+\frac{1}{2}C_{9}C_{10}K_{f}\big{(}\tau^{2}+H^{7}\big{)}\|e^{k}_{h}\|_{\mathcal{A},h}.$

The other terms in (4.8) could be estimated similarly as in the proof of
Theorem 3.12. Then, substituting (4.9) into (4.8) gives us

$\displaystyle\|e_{h}^{n}\|_{\mathcal{A},h}^{2}$ $\displaystyle\leq
2\sqrt{3}K_{f}\sum^{n}_{k=1}\|e_{h}^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|e_{h}^{m}\|_{\mathcal{A},h}+C_{9}C_{10}K_{f}\big{(}\tau^{2}+H^{7}\big{)}\sum^{n}_{k=1}\|e^{k}_{h}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}$
$\displaystyle\qquad+2\sqrt{3}\sum^{n}_{k=1}\|e_{h}^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|R^{m}\|_{h}.$

As was done before, choosing an integer $n^{*}$ $(1\leq n^{*}\leq n)$ such
that $\|e_{h}^{n^{*}}\|_{\mathcal{A},h}=\max_{1\leq k\leq
n}\|e_{h}^{n}\|_{\mathcal{A},h}$, we have

$\displaystyle\|e_{h}^{n}\|_{\mathcal{A},h}$ $\displaystyle\leq
2\sqrt{3}K_{f}\sum^{n}_{k=1}\|e_{h}^{k}\|_{\mathcal{A},h}\sum^{k}_{m=1}\theta^{(k)}_{k-m}+C_{9}C_{10}K_{f}\big{(}\tau^{2}+H^{7}\big{)}\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}+2\sqrt{3}\sum^{n}_{k=1}\sum^{k}_{m=1}\theta^{(k)}_{k-m}\|R^{m}\|_{h}.$

Therefore, for $\tau\leq 1/(4\sqrt{3}K_{f})$, we use Lemma 3.8 and estimate
(3.29) to obtain

$\displaystyle\|e_{h}^{n}\|_{\mathcal{A},h}$ $\displaystyle\leq
4\sqrt{3}K_{f}\sum^{n-1}_{k=1}\tau_{k}\|e^{k}\|_{\mathcal{A},h}+(8\sqrt{3}C_{5}+4\sqrt{3}C_{6}T+2C_{9}C_{10}K_{f}T)\tau^{2}$
(4.10) $\displaystyle\qquad+12\sqrt{3}C_{7}Th^{4}+2C_{9}C_{10}K_{f}TH^{7},$

from which an application of the discrete Grönwall inequality in Lemma 3.9 and
Lemma 2.1 yield the following estimate

$\displaystyle\|e_{h}^{n}\|_{h}$
$\displaystyle\leq\sqrt{3}\|e_{h}^{n}\|_{\mathcal{A},h}$ $\displaystyle\leq
2\sqrt{3}\exp(4\sqrt{3}K_{f}T)\left((4\sqrt{3}C_{5}+2\sqrt{3}C_{6}T+C_{9}C_{10}K_{f}T)\tau^{2}+6\sqrt{3}C_{7}Th^{4}+C_{9}C_{10}K_{f}TH^{7}\right)$
$\displaystyle\leq C_{11}(\tau^{2}+h^{4}+H^{7}).$

The proof is completed.

##  5 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
$\bar{\omega}_{\kappa}$

$\mathcal{V}^{p}_{\kappa}=\left\\{v|v\in\mathcal{V}_{\kappa}\ \ \text{and}\ \
v\ \text{is periodic}\right\\}.$

Furthermore, for any grid functions $w,q\in\mathcal{V}^{p}_{\kappa}$, the
discrete inner product and corresponding norms are redefined as

$(w,q)_{\kappa}=\kappa_{x}\kappa_{y}\sum^{N^{\kappa}_{x}}_{i=1}\sum^{N^{\kappa}_{y}}_{j=1}w_{i,j}q_{i,j},~{}~{}\|w\|_{\kappa}=\sqrt{(w,w)_{\kappa}},~{}~{}\|w\|_{\mathcal{A},\kappa}=\sqrt{(\mathcal{A}_{\kappa}w,w)_{\kappa}},$

and a useful lemma is listed below.

######  Lemma 5.18 ([28])

For any $w\in\mathcal{V}^{p}_{\kappa}$, we have
$\frac{2}{3}\|w\|_{\kappa}\leq\|w\|_{\mathcal{A},\kappa}\leq\|w\|_{\kappa}.$

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^{\infty}$ norms.

###### Lemma 5.19

For any $w\in\mathcal{V}_{H}^{p}$, the following estimate holds

$\displaystyle\|\Pi_{H}w\|_{h}\leq C_{3}\|w\|_{H},$

where $C_{3}=4\left(\frac{3+27^{2}+(10+7\sqrt{7})^{2}}{27^{2}}\right)$.

Proof. Denote $\xi_{i,j}=\Pi_{H,x}w_{i,j}$, and then $\|\Pi_{H}w\|^{2}_{h}$
can be rewritten as

$\displaystyle\|\Pi_{H}w\|^{2}_{h}=h_{x}h_{y}\sum^{N^{h}_{x}}_{i=1}\sum^{N^{h}_{y}}_{j=1}\left(\Pi_{H}w\right)_{i,j}^{2}=h_{x}\sum^{N^{h}_{x}}_{i=1}\left(I_{1}+I_{2}+I_{4}\right),$
(5.1)

where $I_{1}$, $I_{2}$ are defined in Lemma 2.5 and

$\displaystyle
I_{4}:=h_{y}\sum^{N^{h}_{y}}_{j=M_{y}(N^{H}_{y}-1)+1}\left(\Pi_{H,y}\xi_{i,j}\right)^{2}.$

While, the application of Lemma 2.4 shows

$\displaystyle I_{4}\leq
4H_{y}\left(\frac{3}{27^{2}}\xi_{i,N^{H}_{y}-3}^{2}+\frac{(7\sqrt{7}-10)^{2}}{27^{2}}\xi_{i,N^{H}_{y}-2}^{2}+\frac{(7\sqrt{7}+10)^{2}}{27^{2}}\xi_{i,N^{H}_{y}-1}^{2}+\xi_{i,N^{H}_{y}}^{2}\right).$

This together with (2.6)–(2.7) implies

$\displaystyle I_{1}+I_{2}+I_{4}\leq
C_{3}H_{y}\sum^{N^{H}_{y}-1}_{j=1}\xi_{i,j}^{2}+H_{y}\frac{24+8\times
27^{2}}{27^{2}}\xi_{i,N^{H}_{y}}^{2}\leq
C_{3}H_{y}\sum^{N^{H}_{y}}_{j=1}\xi_{i,j}^{2}.$ (5.2)

Finally, the claimed result can be derived immediately by following a
completely similar process as (2.10)–(2.12).

###### Lemma 5.20

For any $w\in\mathcal{V}_{H}^{p}$, the following estimate holds

$\displaystyle\|\Pi_{H}w\|_{h,\infty}\leq C_{4}\|w\|_{H,\infty},$

where $C_{4}=\left(\frac{54+2\sqrt{3}}{27}\right)^{2}$.

Proof. Similar as the proof of Lemma 2.6, estimate (2.13) holds for $M_{y}\leq
j^{*}\leq(N^{H}_{y}-1)M_{y}$. Besides, if $1\leq j^{*}\leq M_{y}$, we have

$\displaystyle\|\Pi_{H}w\|_{h,\infty}$
$\displaystyle\leq\frac{7\sqrt{7}+10}{27}\left|\xi_{i^{*},1}\right|+\frac{7\sqrt{7}-10}{27}\left|\xi_{i^{*},2}\right|+\frac{\sqrt{3}}{27}\left|\xi_{i^{*},3}\right|+\frac{\sqrt{3}}{27}\left|\xi_{i^{*},N_{y}^{H}}\right|$
$\displaystyle\leq\frac{54+2\sqrt{3}}{27}\max_{1\leq k\leq
N^{H}_{y}}\left|\xi_{i^{*},k}\right|,$

and similarly, if $(N^{H}_{y}-1)M_{y}\leq j^{*}\leq N^{h}_{y}$, we have

$\displaystyle\|\Pi_{H}w\|_{h,\infty}\leq\frac{54+2\sqrt{3}}{27}\max_{1\leq
k\leq N^{H}_{y}}\left|\xi_{i^{*},k}\right|.$

Finally, analogous to (2.16), we see
$\left|\xi_{i^{*},k}\right|\leq\frac{54+2\sqrt{3}}{27}\max_{1\leq i\leq
N^{H}_{x}}\left|w_{i,k}\right|$ 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

$\mathcal{D}_{2}\mathcal{A}_{H}u^{n}_{H,i,j}-c\Lambda_{H}u^{n}_{H,i,j}=\mathcal{A}_{H}f(u^{n}_{H,i,j})+\mathcal{A}_{H}g^{n}_{i,j},\quad(i,j)\in\omega_{H},$
(5.3)

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

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}u^{n}_{h,i,j}-c\Lambda_{h}u^{n}_{h,i,j}=\mathcal{A}_{h}F^{n}_{i,j}+\mathcal{A}_{h}g^{n}_{i,j},\quad(i,j)\in\omega_{h},$
(5.4)

subject to the initial condition (1.2) and periodic boundary condition.

Following the proofs of Corollary 4.15, Theorems 4.14 and 4.17, together with
Lemmas 5.18–5.20, 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.

###### Theorem 5.21

Assume that the solution of (1.1)–(1.3) satisfy regularity assumption
condition (3.1) and the adjacent temporal stepsize ratios $r_{k}$ satisfy
$0<r_{k}<4.8645$. If the maximum temporal stepsize $\tau=o(H^{\frac{1}{2}})$
and $\tau\leq\frac{1}{6K_{f}}$, the two-grid compact difference scheme
(5.3)–(5.4) admits a unique solution satisfying

$\|U^{n}-u_{h}^{n}\|_{h}\leq C_{12}\left(\tau^{2}+h^{4}+H^{7}\right)\quad{\rm
for}\ 1\leq n\leq N,$

where
$C_{12}:=3\exp(6K_{f}T)\max\\{6C_{5}+3C_{6}T+C_{9}C_{10}K_{f}T,9C_{7}T\\}$.

##  6 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 Newton-type iterative procedure with tolerance
error $1.0\times 10^{-13}$ is performed to solve the nonlinear algebra systems
at each time level.

###  6.1 Accuracy and numerical stability tests on uniform temporal grids

In this subsection, we shall compare the numerical accuracy of the two-grid
compact difference scheme (4.1)–(4.2) with the standard nonlinear scheme
(3.10) as well as the following implicit-explicit scheme

$\displaystyle\mathcal{D}_{2}\mathcal{A}_{h}u^{n}_{i,j}-c\Lambda_{h}u^{n}_{i,j}=\mathcal{A}_{h}f(u^{n,*}_{i,j})+\mathcal{A}_{h}g^{n}_{i,j},$
(6.1)

where

$\begin{array}[]{l}u^{n,*}:=\left\\{\begin{split}&2u^{n-1}-u^{n-2},\qquad
n\geq 2,\\\ &u^{0},\qquad\qquad\quad\quad\ n=1.\end{split}\right.\end{array}$

For this purpose, we consider the following semilinear parabolic equation

$\displaystyle u_{t}-\Delta u=u-u^{3}+g,\quad(x,y)\in(0,1)^{2},\ t\in(0,\pi],$
(6.2)

where the linear part $g(x,y,t)$ is given such that the exact solution is one
of the following three types, i.e.,

  * 1.

Case I:
$u(x,y,t)=\left[5\sin\left(t\right)+2\sin\left(5t\right)\right]\sin\left(2\pi
x\right)\sin\left(2\pi y\right)$;

  * 2.

Case II:
$u(x,y,t)=\left[10\sin\left(t\right)+5\sin\left(2t\right)+2\sin\left(5t\right)+\sin\left(10t\right)\right]\sin\left(2\pi
x\right)\sin\left(2\pi y\right)$;

  * 3.

Case III:
$u(x,y,t)=\left[10\sin\left(t\right)+50\sin\left(2t\right)+30\sin\left(5t\right)+10\sin\left(10t\right)\right]\sin\left(2\pi
x\right)\sin\left(2\pi y\right)$.

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.

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

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.

Table 1: Spatial convergence of nonlinear scheme (3.10), two-grid scheme (4.1)–(4.2) and implicit-explicit scheme (6.1) for Case I on uniform temporal grids | Nonlinear scheme | Two-grid scheme | Implicit-explicit scheme
---|---|---|---
$(N,N^{h}_{x})$ | Error | Order | Error | Order | Error | Order
$(80,100)$ | $8.16\times 10^{-4}$ | — | $8.16\times 10^{-4}$ | — | $9.40\times 10^{-4}$ | —
$(180,150)$ | $1.64\times 10^{-4}$ | 3.96 | $1.64\times 10^{-4}$ | 3.96 | $1.77\times 10^{-4}$ | 4.12
$(320,200)$ | $5.21\times 10^{-5}$ | 3.99 | $5.21\times 10^{-5}$ | 3.99 | $5.49\times 10^{-5}$ | 4.07
$(500,250)$ | $2.14\times 10^{-5}$ | 3.99 | $2.14\times 10^{-5}$ | 3.99 | $2.23\times 10^{-5}$ | 4.04
Table 2: Spatial convergence of nonlinear scheme (3.10), two-grid scheme (4.1)–(4.2) and implicit-explicit scheme (6.1) for Case II on uniform temporal grids | Nonlinear scheme | Two-grid scheme | Implicit-explicit scheme
---|---|---|---
$(N,N^{h}_{x})$ | Error | Order | Error | Order | Error | Order
$(80,100)$ | $2.23\times 10^{-3}$ | — | $2.23\times 10^{-3}$ | — | Inf | —
$(180,150)$ | $4.81\times 10^{-4}$ | 3.78 | $4.81\times 10^{-4}$ | 3.78 | Inf | —
$(320,200)$ | $1.55\times 10^{-4}$ | 3.94 | $1.55\times 10^{-4}$ | 3.94 | Inf | —
$(500,250)$ | $6.40\times 10^{-5}$ | 3.97 | $6.40\times 10^{-5}$ | 3.97 | $6.77\times 10^{-5}$ | —
Table 3: Spatial convergence of nonlinear scheme (3.10), two-grid scheme (4.1)–(4.2) and implicit-explicit scheme (6.1) for Case III on uniform temporal grids | Nonlinear scheme | Two-grid scheme | Implicit-explicit scheme
---|---|---|---
$(N,N^{h}_{x})$ | Error | Order | Error | Order | Error | Order
$(80,100)$ | $1.84\times 10^{-2}$ | — | $1.84\times 10^{-2}$ | — | Inf | —
$(180,150)$ | $4.03\times 10^{-3}$ | 3.75 | $4.03\times 10^{-3}$ | 3.75 | Inf | —
$(320,200)$ | $1.30\times 10^{-3}$ | 3.92 | $1.30\times 10^{-3}$ | 3.92 | Inf | —
$(500,250)$ | $5.38\times 10^{-4}$ | 3.96 | $5.38\times 10^{-4}$ | 3.96 | Inf | —

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|\rightarrow\infty$), 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.

Table 4: Temporal convergence of nonlinear scheme (3.10), two-grid scheme (4.1)–(4.2) and implicit-explicit scheme (6.1) for Case I on uniform temporal grids | Nonlinear scheme | Two-grid scheme | Implicit-explicit scheme
---|---|---|---
$N$ | Error | Order | Error | Order | Error | Order
$32$ | $4.65\times 10^{-3}$ | — | $4.65\times 10^{-3}$ | — | $5.25\times 10^{-3}$ | —
$64$ | $1.26\times 10^{-3}$ | 1.88 | $1.26\times 10^{-3}$ | 1.88 | $1.48\times 10^{-3}$ | 1.82
$128$ | $3.23\times 10^{-4}$ | 1.97 | $3.23\times 10^{-4}$ | 1.97 | $3.56\times 10^{-4}$ | 2.06
$256$ | $8.13\times 10^{-5}$ | 1.99 | $8.13\times 10^{-5}$ | 1.99 | $8.63\times 10^{-5}$ | 2.05
Table 5: Temporal convergence of nonlinear scheme (3.10), two-grid scheme (4.1)–(4.2) and implicit-explicit scheme (6.1) for Case II on uniform temporal grids | Nonlinear scheme | Two-grid scheme | Implicit-explicit scheme
---|---|---|---
$N$ | Error | Order | Error | Order | Error | Order
$128$ | $9.30\times 10^{-4}$ | — | $9.30\times 10^{-4}$ | — | Inf | —
$256$ | $2.41\times 10^{-4}$ | 1.95 | $2.41\times 10^{-4}$ | 1.95 | Inf | —
$512$ | $6.10\times 10^{-5}$ | 1.98 | $6.10\times 10^{-5}$ | 1.98 | $6.45\times 10^{-5}$ | —
$1024$ | $1.53\times 10^{-5}$ | 1.99 | $1.53\times 10^{-5}$ | 1.99 | $1.61\times 10^{-5}$ | 2.01
Table 6: Temporal convergence of nonlinear scheme (3.10), two-grid scheme (4.1)–(4.2) and implicit-explicit scheme (6.1) for Case III on uniform temporal grids | Nonlinear scheme | Two-grid scheme | Implicit-explicit scheme
---|---|---|---
$N$ | Error | Order | Error | Order | Error | Order
$128$ | $7.70\times 10^{-3}$ | — | $7.70\times 10^{-3}$ | — | Inf | —
$256$ | $2.00\times 10^{-3}$ | 1.94 | $2.00\times 10^{-3}$ | 1.94 | Inf | —
$512$ | $5.13\times 10^{-4}$ | 1.96 | $5.13\times 10^{-4}$ | 1.96 | Inf | —
$1024$ | $1.29\times 10^{-4}$ | 1.99 | $1.29\times 10^{-4}$ | 1.99 | Inf | —

###  6.2 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}\times(0,1]$ with
$c=\frac{1}{8\pi^{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\pi x)\sin(2\pi y)$.
The variable-step temporal grids are generated randomly by

$\tau_{k}:=T\frac{\theta_{k}}{S},\quad{\rm with}\ S=\sum^{N}_{k=1}\theta_{k},$

where $\theta_{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
fourth-order accuracy in space and second-order accuracy in time as proved in
Theorems 3.12 and 4.17. 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!

Table 7: Spatial convergence of nonlinear scheme (3.10) and two-grid scheme (4.1)–(4.2) on variable-step temporal grids  |  | Nonlinear scheme | Two-grid scheme
---|---|---|---
$(N,N^{h}_{x})$ | $\max r_{k}$ | Error | Order | Error | Order
$(80,100)$ | 4.3125 | $1.76\times 10^{-5}$ | — | $1.76\times 10^{-5}$ | —
$(180,150)$ | 4.4777 | $3.67\times 10^{-6}$ | 3.87 | $3.67\times 10^{-6}$ | 3.87
$(320,200)$ | 4.5252 | $1.14\times 10^{-6}$ | 4.06 | $1.14\times 10^{-6}$ | 4.06
$(500,250)$ | 4.5481 | $4.65\times 10^{-7}$ | 4.02 | $4.65\times 10^{-7}$ | 4.02
Table 8: Temporal convergence of nonlinear scheme (3.10) and two-grid scheme (4.1)–(4.2) on variable-step temporal grids  |  | Nonlinear scheme | Two-grid scheme
---|---|---|---
$(N,N^{h}_{x})$ | $\max r_{k}$ | Error | Order | Error | Order
$(40,40)$ | 3.9128 | $7.36\times 10^{-5}$ | — | $6.17\times 10^{-5}$ | —
$(80,80)$ | 3.9673 | $1.80\times 10^{-5}$ | 2.03 | $1.82\times 10^{-5}$ | 1.76
$(160,160)$ | 4.4221 | $4.66\times 10^{-6}$ | 1.95 | $4.53\times 10^{-6}$ | 2.01
$(320,320)$ | 4.5826 | $1.16\times 10^{-6}$ | 2.01 | $1.18\times 10^{-6}$ | 1.94
Table 9: Errors and CPU times of nonlinear scheme (3.10) and two-grid scheme (4.1)–(4.2) on variable-step temporal grids  |  | Nonlinear scheme | Two-grid scheme
---|---|---|---
$(N,N^{h}_{x})$ | $N^{H}_{x}$ | Error | CPU times | Error | CPU times
(240, 240) | 60 | $2.05\times 10^{-6}$ | 33 m 42 s | $2.05\times 10^{-6}$ | 7 m 17 s
(320, 320) | 80 | $1.15\times 10^{-6}$ | 1 h 7 m 20 s | $1.15\times 10^{-6}$ | 22 m 35 s
(400, 400) | 100 | $7.47\times 10^{-7}$ | 1 h 52 m 28 s | $7.47\times 10^{-7}$ | 41 m 27 s
(480, 480) | 120 | $5.19\times 10^{-7}$ | 2 h 27 m 7 s | $5.19\times 10^{-7}$ | 1 h 2 m 18 s

###  6.3 Effectiveness of adaptive temporal stepsize strategy

In this test, we consider model (1.1)–(1.3) on $(0,1)^{2}\times(0,4]$ with
$c=1$ and $f(u)=\sin u$, and the exact solution is chosen as
$u(x,y,t)=\big{[}1+20e^{-40(t-1)^{2}}+30e^{-60(t-4)^{2}}\big{]}\sin(2\pi
x)\sin(2\pi 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 two-grid scheme
based on the adaptive temporal stepsize strategy [17, 25] will be adopted to
improve the temporal accuracy

$\displaystyle\tau_{n+1}=\min\left\\{\max\left\\{\tau_{\min},\frac{\tau_{\max}}{\sqrt{1+\eta\|\partial_{\tau}u^{n}\|_{h}^{2}}}\right\\},r_{\max}\tau_{n}\right\\},$
(6.3)

where $\partial_{\tau}u^{n}:=\nabla_{\tau}u^{n}/\tau_{n}$ and $r_{\max}=4.8$
which satisfies the restrictions in Theorems 3.12 and 4.17. Here $\tau_{\min}$
and $\tau_{\max}$ are the pre-determined minimum and maximum temporal stepsize
and $\eta$ is a pre-chosen parameter. In this example, we uniformly set
$\tau_{\max}=0.2$, $\eta=500$ and gradually reduce $\tau_{\min}$ to generate
grids with distinct stepsizes.

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 Table 10: Temporal convergence of two-grid scheme and adaptive two-grid scheme  |  | Two-grid scheme | Adaptive two-grid scheme
---|---|---|---
$\tau_{\min}$ | N | CPU times | Error | Order | CPU times | Error | Order
0.02 | 67 | 27.96 s | $3.06\times 10^{-1}$ | — | 28.73 s | $2.18\times 10^{-2}$ | —
0.01 | 114 | 45.63 s | $9.93\times 10^{-2}$ | 2.12 | 48.43 s | $4.63\times 10^{-3}$ | 2.91
0.005 | 203 | 1 m 21 s | $2.53\times 10^{-2}$ | 2.37 | 1 m 22 s | $1.01\times 10^{-3}$ | 2.64
0.002 | 489 | 3 m 13 s | $3.12\times 10^{-3}$ | 2.38 | 3 m 19 s | $1.41\times 10^{-4}$ | 2.24
0.001 | 808 | 5 m 30 s | $1.00\times 10^{-3}$ | 2.27 | 5 m 28 s | $3.46\times 10^{-5}$ | 2.80

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.

###  6.4  Application to phase-field Allen–Cahn equation

In this subsection, we consider the following Allen–Cahn equation with a
polynomial double-well potential, subject to periodic boundary condition

$\displaystyle u_{t}-\varepsilon^{2}\Delta u=u-u^{3},\quad(x,y)\in\Omega,\
t\in(0,T],$

where $\varepsilon$ is the interaction length that describes the thickness of
the transition boundary between materials. It is well known that the energy
dissipation law [38, 14]

$\displaystyle E[u](t)\leq E[u](s),\quad\forall t>s$

holds for the Allen–Cahn equation, where $E[u](t)$ represents the Lyapunov
energy functional, namely

$\displaystyle E[u](t)=\int_{\Omega}\frac{\varepsilon^{2}}{2}|\nabla
u|^{2}+F(u)d\mathbf{x}\ \ {\rm with}\ \ F(u)=\frac{1}{4}(1-u^{2})^{2}.$

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 $\mathcal{E}[u^{n}]$ as

$\displaystyle\mathcal{E}[u^{n}]:=-\frac{\varepsilon^{2}}{2}h_{x}h_{y}\sum^{N^{h}_{x}}_{i=1}\sum^{N^{h}_{y}}_{j=1}u^{n}_{i,j}\Delta_{h}u^{n}_{i,j}+\frac{1}{4}h_{x}h_{y}\sum^{N^{h}_{x}}_{i=1}\sum^{N^{h}_{y}}_{j=1}\left(1-(u^{n}_{i,j})^{2}\right)^{2}.$

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 $\partial_{\tau}u^{n}$
with $\partial_{\tau}\mathcal{E}[u^{n}]$.

###### Example 6.1

In this example, we set $\varepsilon=0.02$ and apply these three methods to
simulate the merging of four bubbles with an initial condition

$\displaystyle u_{0}(x,y)=$
$\displaystyle-\tanh\left(\left(\left(x-0.3\right)^{2}+y^{2}-0.2^{2}\right)/\varepsilon\right)\tanh\left(\left(\left(x+0.3\right)^{2}+y^{2}-0.2^{2}\right)/\varepsilon\right)$
$\displaystyle\times\tanh\left(\left(x^{2}+\left(y-0.3\right)^{2}-0.2^{2}\right)/\varepsilon\right)\tanh\left(\left(x^{2}+\left(y+0.3\right)^{2}-0.2^{2}\right)/\varepsilon\right).$

In this simulation, a $384\times 384$ uniform mesh is taken to discretize the
spatial domain $\Omega=(-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
$\tau=0.1$ until time $T=100$. Parameters in the adaptive temporal stepsize
strategy (6.3) are selected as $\tau_{\min}=0.1$, $\tau_{\max}=1$ and
$\eta=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.

A: Nonlinear scheme with fixed temporal stepsize $\tau=0.1$

B: Two-grid scheme with fixed temporal stepsize $\tau=0.1$

C: Adaptive two-grid scheme with $\tau_{\min}=0.1$, $\tau_{\max}=1$ and
$\eta=3200$ 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:  Evolutions of energy (left) and time steps (right) for the nonlinear scheme, two-grid scheme and adaptive two-grid scheme until time $T=30$ Table 11: CPU times and the total number of temporal steps for three schemes | Nonlinear scheme | Two-grid scheme | Adaptive two-grid scheme
---|---|---|---
T | $N$ | CPU times | $N$ | CPU times | $N$ | CPU times
10 | 100 | 57 m 42 s | 100 | 26 m 10 s | 34 | 9 m 31 s
30 | 300 | 3 h 50 m 55 s | 300 | 1 h 47 m 24 s | 71 | 20 m 2 s
50 | 500 | 5 h 58 m 39 s | 500 | 2 h 50 m 56 s | 96 | 37 m 25 s
100 | 1000 | 9 h 43 m 51 s | 1000 | 5 h 9 m 43 s | 156 | 1 h 2 m 22 s

###### Example 6.2

In this example, we consider the coarsening process governed by the Allen–Cahn
equation with the model parameter $\varepsilon=0.01$ and computational domain
$\Omega=(0,1)^{2}$. Here we choose a random initial condition
$u_{0}(x,y)=-0.05+0.1\times rand(x,y)$.

A: Nonlinear scheme with fixed temporal stepsize $\tau=1$

B: Nonlinear scheme with fixed temporal stepsize $\tau=0.01$

C: Two-grid scheme with fixed temporal stepsize $\tau=0.01$

D: Adaptive two-grid scheme with $\tau_{\min}=0.01$, $\tau_{\max}=1$ and
$\eta=8\times 10^{4}$ 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:  Evolutions of energy (left) and time steps (right) for the nonlinear scheme, two-grid scheme and adaptive two-grid scheme until time $T=100$ Table 12: CPU times and the total number of temporal steps for three schemes |  Nonlinear scheme with $\tau=0.01$ |  Two-grid scheme with $\tau=0.01$ | Adaptive two-grid scheme
---|---|---|---
T | $N$ | CPU times | $N$ | CPU times | $N$ | CPU times
20 | 2000 | 4 h 29 m 59 s | 2000 | 2 h 12 m 22 s | 384 | 19 m 53 s
50 | 5000 | 10 h 33 m 53 s | 5000 | 4 h 38 m 8 s | 588 | 27 m 40 s
80 | 8000 | 15 h 41 m 54 s | 8000 | 6 h 41 m 51 s | 722 | 31 m 35 s
100 | 10000 | 18 h 43 m 29 s | 10000 | 8 h 4 m 26 s | 795 | 32 m 57 s

This simulation is performed under $N_{x}^{h}=N_{y}^{h}=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 two-grid scheme by
implementing the nonlinear algorithm up to $T=0.5$. In Figure 5, it displays
the 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 $\tau=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 $\tau=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 $\tau=0.01$
consumes more than 8 hours, even worse the implementation of the nonlinear
scheme runs nearly 19 hours.

##  7 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^{\infty}$ 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}:=\tau_{k}/\tau_{k-1}<4.8645$ and the maximum
temporal stepsize condition $\tau=o(H^{\frac{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.

## CRediT authorship contribution statement

Bingyin Zhang: Methodology, Formal analysis, Software, Writing- Original
draft. Hongfei Fu: Conceptualization, Supervision, Writing- Reviewing and
Editing, Funding acquisition.

## Declaration of Competing Interest

The authors declare that they have no competing interests.

## References

  * [1] Akrivis, G.: Finite difference discretization of the cubic Schr${\rm\ddot{o}}$dinger equation.  IMA J. Numer. Anal. 13, 115–124 (1993)
  * [2] Akrivis, G., Crouzeix, M., Makridakis, C.: Implicit-explicit multistep methods for quasilinear parabolic equations.  Numer. Math. 82, 521–541 (1999)
  * [3] Akrivis, G., Larsson, S.: Linearly implicit finite element methods for the time-dependent Joule heating problem.  BIT Numer. Math. 45, 429–442 (2005)
  * [4] Becker, J.: A second order backward difference method with variable steps for a parabolic problem.  BIT Numer. Math. 38, 644–662 (1998)
  * [5] Chen, C., Li, K., Chen, Y., Huang, Y.: Two-grid finite element methods combined with Crank-Nicolson scheme for nonlinear Sobolev equations.  Adv. Comput. Math. 45, 611–630 (2019)
  * [6] Chen, W., Wang, X., Yan, Y., Zhang, Z.: A second order BDF numerical scheme with variable steps for the Cahn–Hilliard equation.  SIAM J. Numer. Anal. 57, 495–525 (2019)
  * [7] Chen, Y., Huang, Y., Yu, D.: A two-grid method for expanded mixed finite-element solution of semilinear reaction-diffusion equations.  Int. J. Numer. Meth. Eng. 57, 193–209 (2003)
  * [8] Dawson, C.N., Wheeler, M.F., Woodward, C.S.: A two-grid finite difference scheme for nonlinear parabolic equations.  SIAM J. Numer. Anal. 35, 435–452 (1998)
  * [9] Di, Y., Wei, Y., Zhang, J., Zhao, C.: Sharp error estimate of an implicit BDF2 scheme with variable time steps for the phase field crystal model.  J. Sci. Comput. 92, 65 (2022)
  * [10] Du, Y., Liu, Y., Li, H., Fang, Z., He, S.: Local discontinuous Galerkin method for a nonlinear time-fractional fourth-order partial differential equation.  J. Comput. Phys. 344, 108–126 (2017)
  * [11] Durango, F., Novo, J.: Two-grid mixed finite-element approximations to the Navier-Stokes equations based on a Newton-type step.  J. Sci. Comput. 74, 456–473 (2018)
  * [12] de Frutos, J., Novo, J.: A posteriori error estimation with the p-version of the finite element method for nonlinear parabolic differential equations.  Comput. Method. Appl. M. 191, 4893–4904 (2002)
  * [13] Gear, C., Tu, K.: The effect of variable mesh size on the stability of multistep methods.  SIAM J. Numer. Anal. 11, 1025–1043 (1974)
  * [14] Hou, D., Qiao, Z.: An implicit-explicit second-order BDF numerical scheme with variable steps for gradient flows.  J. Sci. Comput. 94, 39 (2023)
  * [15] Hou, T., Jiang, W., Yang, X., Leng, H.: Two-grid ${P}^{2}_{1}$-${P}_{1}$ mixed finite element methods combined with Crank-Nicolson scheme for a class of nonlinear parabolic equations.  Appl. Numer. Math. 137, 136–150 (2019)
  * [16] Hou, T., Tang, T., Yang, J.: Numerical analysis of fully discretized Crank-Nicolson scheme for fractional-in-space Allen–Cahn equations.  J. Sci. Comput. 72, 1214–1231 (2017)
  * [17] Huang, J., Yang, C., Wei, Y.: Parallel energy-stable solver for a coupled Allen–Cahn and Cahn–Hilliard system.  SIAM J. Sci. Comput. 42, C294–C312 (2020)
  * [18] Li, B., Sun, W.: Error analysis of linearized semi-implicit Galerkin finite element methods for nonlinear parabolic equations.  Int. J. Numer. Anal. Mod. 10, 622–633 (2013)
  * [19] Li, B., Ueda, Y., Zhou, G.: A second-order stabilization method for linearizing and decoupling nonlinear parabolic systems.  SIAM J. Numer. Anal. 58, 2736–2763 (2020)
  * [20] Li, D., Li, X., Mei, M., Yuan, W.: A structure-preserving and variable-step BDF2 Fourier pseudo-spectral method for the two-mode phase field crystal model.  Math. Comput. Simulat. 205, 483–506 (2023)
  * [21] Li, M., Wang, L., Wang, N.: Variable-time-step BDF2 nonconforming VEM for coupled Ginzburg-Landau equations.  Appl. Numer. Math. 186, 378–410 (2023)
  * [22] Li, X., Chen, Y., Chen, C.: An improved two-grid technique for the nonlinear time-fractional parabolic equation based on the block-centered finite difference method.  J. Comput. Math. 40, 455–473 (2021)
  * [23] Li, X., Rui, H.: A two-grid block-centered finite difference method for the nonlinear time-fractional parabolic equation.  J. Sci. Comput. 72, 863–891 (2017)
  * [24] Li, X., Shen, J., Rui, H.: Energy stability and convergence of SAV block-centered finite difference method for gradient flows.  Math. Comput. 88, 2047–2068 (2019)
  * [25] Liao, H., Ji, B., Wang, L., Zhang, Z.: Mesh-robustness of an energy stable BDF2 scheme with variable steps for the Cahn–Hilliard model.  J. Sci. Comput. 92, 52 (2022)
  * [26] Liao, H., Ji, B., Zhang, L.: An adaptive BDF2 implicit time-stepping method for the phase field crystal model.  IMA J. Numer. Anal. 42, 649–679 (2022)
  * [27] Liao, H., Song, X., Tang, T., Zhou, T.: Analysis of the second-order BDF scheme with variable steps for the molecular beam epitaxial model without slope selection.  Sci. China. Math. 63, 887–902 (2021)
  * [28] Liao, H., Sun, Z., Shi, H.: Error estimate of fourth-order compact scheme for linear Schr${\rm\ddot{o}}$dinger equations.  SIAM J. Numer. Anal. 47, 4381–4401 (2010)
  * [29] Liao, H., Tang, T., Zhou, T.: On energy stable, maximum-principle preserving, second-order BDF scheme with variable steps for the Allen–Cahn equation.  SIAM J. Numer. Anal. 58, 2294–2314 (2020)
  * [30] Liao, H., Zhang, Z.: Analysis of adaptive BDF2 scheme for diffusion equations.  Math. Comput. 90, 1207–1226 (2021)
  * [31] Liu, Q., Jing, J., Yuan, M., Chen, W.: A positivity-preserving, energy stable BDF2 scheme with variable steps for the Cahn–Hilliard equation with logarithmic potential.  J. Sci. Comput. 95, 37 (2023)
  * [32] Liu, W., Cui, J.: A two-grid block-centered finite difference algorithm for nonlinear compressible Darcy-Forchheimer model in porous media.  J. Sci. Comput. 74, 1786–1815 (2018)
  * [33] Marion, M., Xu, J.: Error estimates on a new nonlinear Galerkin method based on two-grid finite elements.  SIAM J. Numer. Anal. 32, 1170–1184 (1995)
  * [34] Qiao, Z., Zhang, Z., Tang, T.: An adaptive time-stepping strategy for the molecular beam epitaxy models.  SIAM J. Sci. Comput. 33, 1395–1414 (2011)
  * [35] Quarteroni, A., Sacco, R., Saleri, F.: Numerical Mathematics.  Springer (2007)
  * [36] Rui, H., Liu, W.: A two-grid block-centered finite difference method for Darcy-Forchheimer flow in porous media.  SIAM J. Numer. Anal. 53, 1941–1962 (2015)
  * [37] Shampine, L., Reichelt, M.: The MATLAB ODE suite.  SIAM J. Sci. Comput. 18, 1–22 (1997)
  * [38] Shen, J., Xu, J., Yang, J.: A new class of efficient and robust energy stable schemes for gradient flows.  SIAM Rev. 61, 474–506 (2019)
  * [39] Sun, H., Zhao, X., Cao, H., Yang, R., Zhang, M.: Stability and convergence analysis of adaptive BDF2 scheme for the Swift-Hohenberg equation.  Commun. Nonlinear. Sci. 111, 106412 (2022)
  * [40] Sun, Z.: Numerical methods of partial differential equations, 2th edn.  Science Press, Beijing (2012)
  * [41] Wang, H.: An optimal-order error estimate for a family of ELLAM-MFEM approximations to porous medium flow.  SIAM J. Numer. Anal. 46, 2133–2152 (2008)
  * [42] Wang, N., Li, M.: Unconditional error analysis of a linearized BDF2 virtual element method for nonlinear Ginzburg-Landau equation with variable time step.  Commun. Nonlinear. Sci. 116, 106889 (2023)
  * [43] Wang, W., Chen, Y., Fang, H.: On the variable two-step IMEX BDF method for parabolic integro-differential equations with nonsmooth initial data arising in finance.  SIAM J. Numer. Anal. 57, 1289–1317 (2019)
  * [44] Wang, Y., Chen, Y., Huang, Y., Yi, H.: A family of two-grid partially penalized immersed finite element methods for semi-linear parabolic interface problems.  J. Sci. Comput. 88, 80 (2021)
  * [45] Xu, J.: Two-grid discretization techniques for linear and nonlinear PDEs.  SIAM J. Numer. Anal. 33, 1759–1777 (1996)
  * [46] Xu, J., Xie, S., Fu, H.: A two-grid block-centered finite difference method for the nonlinear regularized long wave equation.  Appl. Numer. Math. 172, 128–148 (2022)
  * [47] Zhao, C., Liu, N., Ma, Y., Zhang, J.: Unconditionally optimal error estimate of a linearized variable-time-step BDF2 scheme for nonlinear parabolic equations.  Commun. Math. Sci. 21, 775–794 (2023)
  * [48] Zhao, X., Sun, Z., Hao, Z.: A fourth-order compact ADI scheme for two-dimensional nonlinear space fractional Schr${\rm\ddot{o}}$dinger equation.  SIAM J. Sci. Comput. 36, A2865–A2886 (2014)