# Efficient Exponential Integrator Finite Element Method for
Semilinear Parabolic Equations

Jianguo Huang111J. Huang’s research was partially supported by National
Natural Science Foundation of China under grant number 11571237. L. Ju’s
research was partially supported by US National Science Foundattion under
grant number DMS-2109633 and by US Department of Energy under grant number DE-
SC0020270<EMAIL_ADDRESS>School of Mathematical Sciences and MOE-LSC,
Shanghai Jiao Tong University, Shanghai 200240, China  Lili Ju222L. Ju’s
research was partially supported by US National Science Foundattion under
grant number DMS-2109633 and by US Department of Energy under grant number DE-
SC0020270<EMAIL_ADDRESS>Department of Mathematics, University of South
Carolina, Columbia, SC 29208, USA  Yuejin Xu<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose an efficient exponential integrator finite element
method for solving a class of semilinear parabolic equations in rectangular
domains. The proposed method first performs the spatial discretization of the
model equation using the finite element approximation with continuous
multilinear rectangular basis functions, and then takes the explicit
exponential Runge-Kutta approach for time integration of the resulting semi-
discrete system to produce fully-discrete numerical solution. Under certain
regularity assumptions, error estimates measured in $H^{1}$-norm are
successfully derived for the proposed schemes with one and two RK stages. More
remarkably, the mass and coefficient matrices of the proposed method can be
simultaneously diagonalized with an orthogonal matrix, which provides a fast
solution process based on tensor product spectral decomposition and fast
Fourier transform. Various numerical experiments in two and three dimensions
are also carried out to validate the theoretical results and demonstrate the
excellent performance of the proposed method.

###### keywords:

Semilinear parabolic equations, exponential integrator, finite element method,
fast Fourier transform, Runge-Kutta, error estimates

##  1 Introduction

In this paper, we are devoted to studying numerical solution of the semilinear
parabolic equation taking the following form:

$u_{t}=D\Delta u+f(t,u),\quad\;\;\bm{x}\in\Omega,\;t\geq 0,$ (1.1)

where $\Omega$ is an open rectangular domain in $\mathbb{R}^{d}$, $D>0$ is the
constant diffusion coefficient, $u(t,\bm{x})$ is the unknown function and
$f(t,u)$ is the nonlinear reaction term of the underlying system. It is well
known that the equation (1.1) together with appropriate initial value and
boundary condition has been widely used in mathematical models for various
scientific and engineering applications, such as the Allen-Cahn equation and
some other phase field models [11] for describing the phase transition and
separation [1], the time-dependent advection diffusion equation and the
Navier-Stokes equations for fluid dynamics [39], and the Ginzburg-Landau
equation [12] for superconductivity. Many related numerical methods,
especially for its temporal discretization, have been proposed and analyzed in
the past few decades. Some typical ones include fully-implicit scheme[16]
implicit-explicit method [2, 43], integrating factor (IF) method [28, 34, 41],
split step method [3, 37], sliders methods [10], exponential time differencing
(ETD) method [8, 32], invariant energy quadratization (IEQ) method [44],
integrating factor Runge-Kutta (IFRK) method [24, 29], scalar auxiliary
variable (SAV) method [38] and so on.

The IF, ETD and IFRK methods fall into the category of exponential integrator-
based methods and have received much attention in the past two decades due to
their effectiveness and stability in handling stiff semilinear systems, such
as the equation (1.1) with very small $D$ and highly nonlinear $f$. All the
three methods can maintain the exponential behavior of the system in some
extent but also have significant differences. The IF and ETD methods both
start with the Duhamel formula, but they differ on the way of evaluating the
time integration of the terms resulting from multiplication of the exponential
integrator and $f$. The IF method uses quadrature rules directly to the whole
integrand, while the ETD method first approximates the nonlinear part using
polynomial interpolation and then performs exact integration on the resulting
integrands. Instead, the IFRK method first applies an exponential-type
variable transformation to eliminate the linear term from the equation and
then use the standard Runge-Kutta method for time integration of the
transformed system. The $L^{2}$-stability and $L^{\infty}$-stability of ETD
methods were established in [15] by means of Fourier transformation (FFT) and
maximum bound principle of the underlying equations. Some researches have been
done to further enhance the stability of ETD method [15, 27, 33] and improve
their numerical advantages by exploring the stability factor numerically [8].
Motivated by [35], a family of exponential Runge-Kutta methods [19, 18, 20,
31] were systematically developed with corresponding consistency and order
conditions based on the theory of semigroup. Exponential integrator-based
methods have also been extended to various areas, such as solution of systems
with nondiagonal operators [27], Rosenbrock-type method [21], low-regularity
integrators [36], adaptive time stepping [42], and general order parabolic
equations [22, 23].

A significant feature of the exponential integrator-based methods is that they
all require evaluations of the products of matrix exponentials and vectors.
Such operations generally are implemented with Krylov subspace method [14],
which is still time-consuming for large scale systems. To overcome this
difficulty, some fast FFT-based algorithms were proposed for efficient
implementation of ETD methods in rectangular domains, but they only work for
the case of using finite difference method for spatial discretization [26, 45,
13]. The key idea is to utilize tensor product representation of discrete
spatial operators and spectral decomposition in each dimension. Moreover,
convergence analysis of ETD methods to the equation (1.1) highly depends on
discrete maximum bound principles, which does not always hold for the classic
finite element or pseudo-spectral discretizations, see [14, 25] and references
cited therein.

In this paper, we propose an efficient exponential integrator finite element
(EIFE) method for solving the semilinear parabolic equation (1.1) in
rectangular domains. In the proposed method, we first carry out the spatial
discretization by using continuous multilinear rectangular finite elements to
obtain a semi-discrete (in space) system, then apply the explicit exponential
Runge-Kutta approach for time integration of the resulting system to achieve
the fully-discrete scheme. Error estimates measured in $H^{1}$-norm are also
successfully derived for the EIFE schemes with one and two RK stages for the
problem with Dirichlet boundary condition. Specifically, we first rewrite the
semi-discrete system as a finite dimensional spatial operator equation and
then error estimate of the semi-discrete solution is naturally obtained by
using the variational formulation and the energy method. Next, the fully-
discrete system is regarded as a finite dimensional evolution equation, we
then adapt and verify the arguments from [18] to estimate the errors of the
fully-discrete solutions. Furthermore, it is observed that the resulting mass
and coefficient matrices in the EIFE method can be diagonalized simultaneously
with an orthogonal matrix, whose multiplication with a vector can be
efficiently computed using FFT, which is held for both Dirichlet and periodic
boundary condition cases. Consequently, following [26], fast solution
algorithms are able to be constructed based on the tensor product spectral
decomposition. To the best of our knowledge, the work presented in this paper
is the first study on fast numerical method with rigorous fully-discrete error
estimates, which combines exponential integrator in time and finite element
discretization in space.

The rest of the paper is organized as follows. The EIFE method is first
proposed in Section 2, and then its fully-discrete error analysis is given in
Section 3. The FFT-based fast implementation of EIFE method is illustrated in
Section 4. In Section 5, various numerical experiments are carried out to
validate the theoretical results and demonstrate the excellent performance of
the EIFE method. Finally, some concluding remarks are drawn in Section 6.

##  2 An exponential integrator finite element method

In this section, we will propose an exponential integrator finite element
method for solving the model equation (1.1). Let us start with some standard
notations for later representations. Given a bounded Lipschitz domain
$G\subset\mathbb{R}^{d}$ and a non-negative integer $s\geq 0$, denote by
$H^{s}(G)$ the standard Sobolev spaces on $G$ with the norm $\|\cdot\|_{s,G}$
and the semi-norm $|\cdot|_{s,G}$, by $(\cdot,\,\cdot)_{G}$ the $L^{2}$-inner
product on $G$, and by $H_{0}^{s}(G)$ the closure of $C_{0}^{\infty}(G)$ with
respect to the norm $\|\cdot\|_{s,G}$. We also define $\|v\|_{k,\infty,G}={\rm
ess\ sup}_{|\bm{\alpha}|\leq k}\|D^{\bm{\alpha}}v\|_{L^{\infty}(\Omega)}$ for
any function $v$ such that the right-hand side term makes sense, where
$\bm{\alpha}=(\alpha_{1},\cdots,\alpha_{d})$ is a multi-index and
$|\bm{\alpha}|=\alpha_{1}+\cdots+\alpha_{d}$. We may omit the subscript for
$G=\Omega$ when there is no confusion caused. For any non-negative integer
$\ell\geq 0$, $\mathbb{P}_{\ell}(G)$ stands for the set of all polynomials on
$G$ with the total degree at most $\ell$. Moreover, given two quantities $a$
and $b$, “$a\lesssim{b}$” abbreviates “$a\leq{C}b$”, where $C$ is a generic
positive constant; “$a\eqsim b$” is equivalent to “$a\lesssim b\lesssim a$”.

Let us consider the semilinear parabolic equation (1.1) with homogeneous
Dirichlet boundary condition and an initial configuration $u_{0}\in
H^{2}(\Omega)\cap H_{0}^{1}(\Omega)$, that is

$\left\\{\begin{split}&u_{t}=D\Delta
u+f(t,u),\quad\;\;\bm{x}\in\Omega,\;\;0\leq t\leq T,\\\
&u(0,\bm{x})=u_{0}(\bm{x}),\qquad\quad\,\bm{x}\in\Omega,\\\
&u(t,\bm{x})=0,\qquad\qquad\quad\;\bm{x}\in{\partial\Omega},\;\;0\leq t\leq
T,\;\;\end{split}\right.$ (2.1)

where $T\geq 0$ is the terminal time.

###  2.1 Semi-discretization in space by finite element approximation

The variational formulation of the above problem (2.1) is to find $u\in
L^{2}(0,T;H_{0}^{1}(\Omega))$ and $u_{t}\in L^{2}(0,T;L^{2}(\Omega))$ such
that

$\left\\{\begin{split}&(u_{t},v)+a(u,v)=(f(t,u),v),\quad\quad\forall\,v\in
H_{0}^{1}(\Omega),\;0\leq t\leq T,\\\ &u(0)=u_{0},\end{split}\right.$ (2.2)

where the bilinear operator $a(\cdot,\cdot)$ is symmetric and defined by

$a(w,v)=\int_{\Omega}D\nabla w\cdot\nabla v\,{\rm
d}\bm{x},\quad\forall\,w,v\in H_{0}^{1}(\Omega).$ (2.3)

We now define a finite element space $V_{h}$ for approximation of
$H_{0}^{1}(\Omega)$. Since $\Omega\in{\mathbb{R}}^{d}$ is a rectangular
domain, let us assume $\overline{\Omega}:=\prod_{i=1}^{d}[a_{i},b_{i}]$. For
each $i=1,\cdots,d$, we make a uniform partition of $[a_{i},b_{i}]$ with the
subinterval size $h_{i}=\frac{b_{i}-a_{i}}{N_{i}}$, to get the nodes
$x_{i}^{j}=a_{i}+jh_{i},j=0,\cdots,N_{i}$ as

$a_{i}=x_{i}^{0}<x_{i}^{1}<\cdots<x_{i}^{N_{i}}=b_{i}.$

With this regular partition, we obtain a one-dimensional continuous piecewise
linear finite element space for $[a_{i},b_{i}]$ as

$\displaystyle V_{h_{i}}^{i}(a_{i},b_{i})$ $\displaystyle:=$
$\displaystyle\left\\{v\in\mathcal{C}[a_{i},b_{i}]:\
v|_{\big{[}x_{i}^{j-1},x_{i}^{j}\big{]}}\in\mathbb{P}_{1}(\big{[}x_{i}^{j-1},x_{i}^{j}\big{]}),1\leq
j\leq N_{i}\right\\}\cap H_{0}^{1}(a_{i},b_{i})$ $\displaystyle=$
$\displaystyle{\rm
span}\left\\{\phi_{i}^{1}(x_{i}),\cdots,\phi_{i}^{N_{i}-1}(x_{i})\right\\},$

where $\phi_{i}^{j}(x_{i})$ is the $j$-th nodal basis function of
$V_{h}^{i}(a_{i},b_{i})$. By using the tensor product of all above finite
element spaces, we can obtain a finite element space for $\Omega$ as follows:

$\displaystyle V_{h}$ $\displaystyle:=$ $\displaystyle
V_{h_{1}}^{1}(a_{1},b_{1})\otimes\cdots\otimes V_{h_{d}}^{d}(a_{d},b_{d})$
$\displaystyle=$ $\displaystyle{\rm
span}\\{{\phi_{1}^{i_{1}}}(x_{1})\cdots\phi_{d}^{i_{d}}(x_{d}):\ 1\leq
i_{1}\leq N_{1}-1,\cdots,1\leq i_{d}\leq N_{d}-1\\}.$

It is evident that $V_{h}\subset H_{0}^{1}(\Omega)$. Define
$h=\textstyle\max_{1\leq i\leq d}h_{i}$ as the mesh size of the corresponding
uniformly rectangular partition $\mathcal{T}_{h}$ for generating $V_{h}$. For
the forthcoming error analysis, we assume the partition $\mathcal{T}_{h}$ is
quasi-uniform, i.e., $h\eqsim h_{i}$ for all $1\leq i\leq d$. The finite
element approximation in space for (2.2) is then to find $u_{h}\in
L^{2}(0,T;V_{h})$ such that

$\left\\{\begin{split}&(u_{h,t},v_{h})+a(u_{h},v_{h})=(f(t,u_{h}),v_{h}),\qquad\forall\,v_{h}\in
V_{h},\;0\leq t\leq T,\\\ &u_{h}(0)=P_{h}u_{0},\end{split}\right.$ (2.5)

where $P_{h}:L^{2}(\Omega)\rightarrow V_{h}$ is the $L^{2}$ orthogonal
projection operator. Since $\mathcal{T}_{h}$ is quasi-uniform, it is easy to
show, using the similar arguments in [4], that $P_{h}$ is stable with respect
to $L^{2}$-norm or $H^{1}$-norm, i.e., $\|P_{h}u\|_{0}\lesssim\|u\|_{0}$ and
$\|P_{h}u\|_{1}\lesssim\|u\|_{1}$ for any $u\in H_{0}^{1}(\Omega)$. Applying
the inverse inequality for finite elements [5, 7] we know

$a(w_{h},v_{h})\lesssim|w_{h}|_{1}|v_{h}|_{1}\lesssim
h^{-2}\|w_{h}\|_{0}\|v_{h}\|_{0},\quad\forall\,w_{h},v_{h}\in V_{h}.$ (2.6)

where the hidden constants are independent of $h$, which means
$a(\cdot,\cdot)$ is a bounded bilinear form over $V_{h}$ with respect to
$L^{2}$-norm. Therefore, by the Riesz representation theorem, there exists a
bounded linear operator $L_{h}:V_{h}\rightarrow V_{h}$ such that

$a(w_{h},v_{h})=(L_{h}w_{h},v_{h}),\quad\forall\,w_{h},v_{h}\in V_{h}.$ (2.7)

Making use of the projection operator $P_{h}$, we can finally reformulate the
problem (2.5) as the following equivalent semi-discrete (in space) system:

$\left\\{\begin{split}&u_{h,t}+L_{h}u_{h}=P_{h}f(t,u_{h}),\quad\;\bm{x}\in\Omega,\;\;0\leq
t\leq T,\\\
&u_{h}(0)=P_{h}u_{0},\qquad\qquad\qquad\;\bm{x}\in\Omega.\end{split}\right.$
(2.8)

###  2.2 Explicit exponential integrator in time

Let us divide the time interval $[0,T]$ into $N_{T}>0$ subintervals
$[t_{n},t_{n+1}]$, $n=0,1,\cdots,N_{T}-1,$ with $\tau_{n}=t_{n+1}-t_{n}>0$
being the time step size at $t_{n}$. Let $\\{e^{-tL_{h}}\\}_{t\geq 0}$ denote
the semigroup on $V_{h}$ with the infinitesimal generator $(-L_{h})$. For
simplicity of presentation, we define $u(t):=u(t,\cdot)$ and
$u_{h}(t):=u_{h}(t,\cdot)$. By the Duhamel principle, the semi-discrete
solution $u_{h}$ to the problem (2.8) can be equivalently expressed as

$u_{h}(t_{n+1})=e^{-\tau_{n}L_{h}}u_{h}(t_{n})+\int_{0}^{\tau_{n}}e^{-(\tau_{n}-\sigma)L_{h}}P_{h}f(t_{n}+\sigma,u_{h}(t_{n}+\sigma))\,{\rm
d}\sigma.$ (2.9)

Denote by $u_{h}^{n}$ the fully-discrete numerical solution at the time step
$t_{n}$ after temporal discretization of (2.9). We then apply the classic
explicit exponential Runge-Kutta approach [20] (see Table 1) to (2.9) and
obtain a fully-discrete numerical method for solving the semilinear parobolic
problem (2.1) as follows: for $n=0,1,\cdots,N_{T}-1$,

$\left\\{\begin{split}u_{h}^{n+1}&=\chi(-\tau_{n}L_{h})u_{h}^{n}+\tau_{n}\textstyle\sum\limits_{i=1}^{s}b_{i}(-\tau_{n}L_{h})G_{ni},\\\
U_{ni}&=\chi_{i}(-\tau_{n}L_{h})u_{h}^{n}+\tau_{n}\textstyle\sum\limits_{j=1}^{i-1}a_{ij}(-\tau_{n}L_{h})G_{nj},\quad
i=1,\cdots,s,\\\ G_{ni}&=P_{h}f(t_{n}+c_{i}\tau_{n},U_{ni}),\quad
i=1,\cdots,s,\end{split}\right.$ (2.10)

where the integer $s>0$ denotes the number of stages for exponential Runge-
Kutta method.

Table 1: Explicit exponential Runge-Kutta tableau for (2.9) with $s$ denoting the number of stages. $c_{1}$ |  |  |  |  | $\chi_{1}(-\tau_{n}L_{h})$
---|---|---|---|---|---
$c_{2}$ | $a_{21}(-\tau_{n}L_{h})$ |  |  |  | $\chi_{2}(-\tau_{n}L_{h})$
$\vdots$ | $\vdots$ | $\ddots$ |  |  | $\vdots$
$c_{s}$ | $a_{s1}(-\tau_{n}L_{h})$ | $\cdots$ | $a_{s,s-1}(-\tau_{n}L_{h})$ |  | $\chi_{s}(-\tau_{n}L_{h})$
| $b_{1}(-\tau_{n}L_{h})$ | $\cdots$ | $b_{s-1}(-\tau_{n}L_{h})$ | $b_{s}(-\tau_{n}L_{h})$ | $\chi(-\tau_{n}L_{h})$

Here the interpolation nodes $c_{1},\cdots,c_{s}$ are $s$ distinct nodes
selected in $[0,1]$, and the weights are given by

$\left\\{\begin{split}b_{i}(-\tau_{n}L_{h})&=\int_{0}^{1}e^{-\tau_{n}(1-\theta)L_{h}}l_{i}(\theta)\,{\rm
d}\theta,\quad i=1,\cdots,s,\\\
a_{ij}(-\tau_{n}L_{h})&=\frac{1}{\tau_{n}}\int_{0}^{c_{i}\tau_{n}}e^{-(c_{i}\tau_{n}-\tau)L_{h}}l_{j}(\tau)\,{\rm
d}\tau,\quad i=1,\cdots,s,\;j=1,\cdots,i-1,\end{split}\right.$

with $\\{l_{i}(\theta)\\}_{i=1}^{s}$ being the classic Lagrange interpolation
polynomials

$l_{i}(\theta):=\prod\limits_{m=1,m\neq i}^{s}\frac{\theta-
c_{m}}{c_{i}-c_{m}},\quad i=1,\cdots s.$

The coefficients $\chi,\chi_{i}$ are constructed from exponential functions or
approximations of such functions evaluated at the operator $(-\tau_{n}L_{h})$.
Let us define

$\varphi_{k}(z)=\int_{0}^{1}e^{(1-\theta)z}\frac{\theta^{k-1}}{(k-1)!}\,{\rm
d}\theta,\;\;k\geq 1.$

It is easy to check that $\varphi_{0}(z)=e^{z}$ and
$\varphi_{k+1}(z)=\frac{\varphi_{k}(z)-\varphi_{k}(0)}{z}$ for $k>0$. Then,

$\qquad\varphi_{k}(-\tau_{n}L_{h})=\frac{1}{\tau_{n}^{k}}\int_{0}^{\tau_{n}}e^{-(\tau_{n}-\tau)L_{h}}\frac{\tau^{k-1}}{(k-1)!}\,{\rm
d}\tau,\;\;k\geq 1.$

To ensure the consistency of the explicit exponential Runge-Kutta method, it
is always assumed that the following conditions hold [20]:

$\left\\{\begin{split}\chi(-\tau_{n}L_{h})=e^{-\tau_{n}L_{h}},\qquad&\\\
\qquad\chi_{1}(-\tau_{n}L_{h})=1,\qquad&\chi_{i}(-\tau_{n}L_{h})=e^{-c_{i}\tau_{n}L_{h}},\quad
1\leq i\leq s,\\\
\textstyle\sum\limits_{j=1}^{s}b_{j}(-\tau_{n}L_{h})=\varphi_{1}(-\tau_{n}L_{h}),\qquad&\textstyle\sum\limits_{j=1}^{i-1}a_{ij}(-\tau_{n}L_{h})=c_{i}\varphi_{1}(-c_{i}\tau_{n}L_{h}),\;\;1\leq
i\leq s.\end{split}\right.$ (2.11)

With the help of the above consistency condition (2.11), we can find $c_{1}=0$
and obtain the explicit forms of the functions $\chi$ and $\chi_{i}$ in
(2.10). Finally, the proposed fully-discrete scheme (2.10) can be rewritten as

$\left\\{\begin{split}u_{h}^{n+1}&=e^{-\tau_{n}L_{h}}u_{h}^{n}+\tau_{n}\textstyle\sum\limits_{i=1}^{s}b_{i}(-\tau_{n}L_{h})G_{ni},\\\
U_{ni}&=e^{-c_{i}\tau_{n}L_{h}}u_{h}^{n}+\tau_{n}\textstyle\sum\limits_{j=1}^{i-1}a_{ij}(-\tau_{n}L_{h})G_{nj},\quad
i=1,\cdots,s,\\\ G_{ni}&=P_{h}f(t_{n}+c_{i}\tau_{n},U_{ni}),\quad
i=1,\cdots,s.\end{split}\right.$ (2.12)

We refer to the above fully-discrete scheme (2.12) as the exponential
integrator finite element (EIFE) method. In particular, if the number of RK
stages $s=1$, then the so-called first-order Euler exponential scheme is
obtained for time stepping, and (2.12) can correspondingly be expressed as

$u_{h}^{n+1}=e^{-\tau_{n}L_{h}}u_{h}^{n}+\tau_{n}\varphi_{1}(-\tau_{n}L_{h})P_{h}f(t_{n},u_{h}^{n}),$
(2.13)

which is called EIFE1. If the number of RK stages $s=2$, then the two
interpolation nodes are taken as $c_{1}=0$ and $c_{2}\in(0,1]$ and the two-
stage second-order exponential Runge-Kutta scheme is obatined for time
stepping, and (2.12) correspondingly reads

$\begin{split}&u_{h}^{n+1}=\;e^{-\tau_{n}L_{h}}u_{h}^{n}+\tau_{n}\big{(}(\varphi_{1}(-\tau_{n}L_{h})-\frac{1}{c_{2}}\varphi_{2}(-\tau_{n}L_{h}))P_{h}f(t_{n},u_{h}^{n})\\\
&+\frac{1}{c_{2}}\varphi_{2}(-\tau_{n}L_{h})P_{h}f(t_{n}+c_{2}\tau_{n},e^{-c_{2}\tau_{n}L_{h}}u_{h}^{n}+c_{2}\tau_{n}\varphi_{1}(-c_{2}\tau_{n}L_{h})P_{h}f(t_{n},u_{h}^{n}))\big{)},\end{split}$
(2.14)

which is denoted as EIFE2.

Fully-discrete error analysis of (2.13) and (2.14) will be carefully studied
in Section 3. For even higher order ($\geq 3$) explicit exponential Runge-
Kutta schemes, more complicated order conditions are needed and we can refer
to [20] for details. It is worth noting that the proposed EIFE method is so
far described in an abstract setting for the convenience of further
theoretical analysis. If we express the finite element numerical solution in
terms of nodal basis functions of $V_{h}$, we can further reformulate the
equations (2.5) and (2.8) in the matrix form. More remarkably, the resulting
mass and coefficient matrices can be diagonalized simultaneously with an
orthogonal matrix and the product of the $e^{-\tau_{n}L_{h}}$ with any vector
consequently can be realized efficiently using tensor product spectral
decomposition and FFT, which render a fast implementation for the EIFE method
(2.12) as described in Section 4.

##  3 Error analysis of the EIFE method

In this section, we focus on error estimates of fully-discrete solutions
produced by the EIFE method (2.12) for solving the semilinear parabolic
problem (2.1) (i.e., the case with homogeneous Dirichlet boundary condition).
Note that the analysis also work for the problem with nonhomogeneous Dirichlet
boundary condition after minor modifications. From now on, we always assume
the dimension of the problem $d\leq 3$.

###  3.1 Some preliminary lemmas

We first present some results related to the semigroup
$\\{e^{-tL_{h}}\\}_{t\geq 0}$ and some terms used in the exponential Runge-
Kutta schemes, which are important to the forthcoming analysis of the proposed
EIFE method.

######  Lemma 3.1.

  * (i)

For any given parameter $\gamma\geq 0$, it holds

$\|e^{-\tau L_{h}}\|_{0}+\|\tau^{\gamma}L_{h}^{\gamma}e^{-\tau
L_{h}}\|_{0}\lesssim 1,\quad\forall\,\tau>0,\,\forall\,h>0.$ (3.1)

  * (ii)

For any given parameter $0\leq\gamma\leq 1$, it holds

$\Big{\|}\tau L_{h}^{\gamma}\textstyle\sum\limits_{j=1}^{n-1}e^{-j\tau
L_{h}}\Big{\|}_{0}\lesssim 1,\quad\forall\,\tau>0,\,\forall\,h>0.$ (3.2)

  * (iii)

For any given parameter $0\leq\gamma\leq 1$, it holds

$\|\phi(-\tau L_{h})\|_{0}+\|\tau^{\gamma}L_{h}^{\gamma}\phi(-\tau
L_{h})\|_{0}\lesssim 1,\quad\forall\,\tau>0,\,\forall\,h>0,$ (3.3)

where $\phi(-\tau L_{h})=b_{i}(-\tau L_{h})$ or $\phi(-\tau
L_{h})=a_{ij}(-\tau L_{h})$, $i,j=1,\cdots,s$.

###### Proof.

(i) Since $a(\cdot,\cdot)$ is a symmetric bilinear operator, we easily know
$L_{h}$ is a linear symmetric operator on $V_{h}$. On the other hand, by the
Poincaré inequality, there exists a constant $\alpha>0$ such that

$\alpha\|v_{h}\|_{0}^{2}\leq a(v_{h},v_{h}),\quad\forall\,v_{h}\in V_{h}.$

This combined with (2.6) gives $0\leq\alpha\leq\lambda\lesssim h^{-2}$ for all
$\lambda\in\lambda(L_{h})$ by means of Rayleigh representation theorem, where
$\lambda(L_{h})$ denotes the set of all eigenvalues of $L_{h}$. Therefore, we
obtain

$\|e^{-\tau L_{h}}\|_{0}\leq e^{-\alpha\tau}<1.$ (3.4)

On the other hand, $L_{h}$ is a symmetric operator on $V_{h}$, so is
$\tau^{\gamma}L_{h}^{\gamma}e^{-\tau L_{h}}$. Thus it holds

$\|\tau^{\gamma}L_{h}^{\gamma}e^{-\tau
L_{h}}\|_{0}=\textstyle\max_{\lambda\in\lambda(L_{h})}|(\tau\lambda)^{\gamma}e^{-\tau\lambda}|.$
(3.5)

Let us consider an auxiliary function $g(x)=x^{\gamma}e^{-x}$ for $x\geq 0$.
The derivative with respect to $x$ is
$g^{\prime}(x)=x^{\gamma-1}e^{-x}(\gamma-x)$, so the maximum of $g(x)$ is
taken at $x=\gamma$, which implies

$g(x)\leq\gamma^{\gamma}e^{-\gamma},\quad\forall\,x\geq 0.$ (3.6)

Combination of (3.6) with (3.5) immediately gives us
$\|\tau^{\gamma}L_{h}^{\gamma}e^{-\tau
L_{h}}\|_{0}\leq\gamma^{\gamma}e^{-\gamma},$ which together with (3.4) then
directly deduces (3.1).

(ii) Follow the similar arguments for deriving (3.5),

$\begin{split}\Big{\|}\tau
L_{h}^{\gamma}\textstyle\sum\limits_{j=1}^{n-1}e^{-j\tau
L_{h}}\Big{\|}_{0}&=\textstyle\max\limits_{\lambda\in\lambda(L_{h})}\Big{|}\tau\lambda^{\gamma}\textstyle\sum\limits_{j=1}^{n-1}e^{-j\tau\lambda}\Big{|}\leq\textstyle\max\limits_{\lambda\in\lambda(L_{h})}\lambda^{\gamma-1}\Big{|}\frac{\tau\lambda}{e^{\tau\lambda}-1}\Big{|}.\end{split}$

Since the auxiliary function $g(x)=\Big{|}\frac{x}{e^{x}-1}\Big{|}\leq 1$ for
any $x>0$, $0\leq\gamma\leq 1$, we have

$\Big{\|}\tau L_{h}^{\gamma}\textstyle\sum\limits_{j=1}^{n-1}e^{-j\tau
L_{h}}\Big{\|}_{0}\lesssim 1.$ (3.7)

(iii) If $\phi(x)=b_{i}(x)$ for $1\leq i\leq s$, then

$\begin{split}\|\phi(-\tau
L_{h})\|_{0}&=\Big{\|}\int_{0}^{1}e^{-\tau(1-\theta)L_{h}}l_{i}(\theta)\,{\rm
d}\theta\Big{\|}_{0}\leq\int_{0}^{1}\|e^{-\tau(1-\theta)L_{h}}\|_{0}|l_{i}(\theta)|\,{\rm
d}\theta\lesssim 1,\end{split}$

since the interpolation nodes $c_{1},\cdots,c_{s}$ are fixed. For all
$\gamma\in[0,1]$, follow the similar arguments, we have

$\begin{split}\|\tau^{\gamma}L_{h}^{\gamma}\phi(-\tau
L_{h})\|_{0}&=\textstyle\max\limits_{\lambda\in\lambda(L_{h})}\Big{|}\tau^{\gamma}\lambda^{\gamma}\int_{0}^{1}e^{-\tau(1-\theta)\lambda}l_{i}(\theta)\,{\rm
d}\theta\Big{|}\\\
&\lesssim\textstyle\max\limits_{\lambda\in\lambda(L_{h})}\Big{|}\tau^{\gamma}\lambda^{\gamma}\int_{0}^{1}e^{-\tau(1-\theta)\lambda}\,{\rm
d}\theta\Big{|}=\textstyle\max\limits_{\lambda\in\lambda(L_{h})}(\tau\lambda)^{\gamma-1}(1-e^{-\tau\lambda}).\end{split}$
(3.8)

Let us consider the auxiliary function $g(x)=x^{\gamma-1}(1-e^{-x})$ for $x>0$
where $0\leq\gamma\leq 1$. Since $\lim\limits_{x\to 0^{+}}f(x)=0$ and
$\lim\limits_{x\to+\infty}f(x)=0$, the function $f(x)$ must take its maximum
at a point $x_{0}\in(0,+\infty)$, and $f^{\prime}(x_{0})=0$ by the Fermat
lemma. Note that

$f^{\prime}(x)=x^{\gamma-2}e^{-x}((\gamma-1)(e^{x}-1)+x),$

which implies that $(\gamma-1)(e^{x_{0}}-1)+x_{0}=0$, and $x_{0}$ depends only
on $\gamma$. Hence, $f(x)$ is bounded by a constant $C(\gamma)>0$. This
combined with the above estimate (3.8) gives

$\|\tau^{\gamma}L_{h}^{\gamma}\phi(-\tau L_{h})\|_{0}\lesssim 1.$ (3.9)

Similarly,the above results are also held for $\phi(x)=a_{ij}(x)$,
$i,j=1,\cdots,s$. ∎

In addition, recalling the definitions (2.3) and (2.7) and noting that $L_{h}$
is a symmetric positive definite operator, we have the following important
property:

$\|v\|_{1}\eqsim\|L_{h}^{\frac{1}{2}}v\|_{0},\quad\forall\ v\in V_{h},$ (3.10)

which establishes an important relation between the $H^{1}$ and $L^{2}$ norms
over $V_{h}$, and will be frequently used later on. Next let us introduce the
mild growth condition for the function $f$ as given in [40] and some
regularity conditions required for the exact solution $u(t)$ in order to carry
out convergence and error analysis of EIFE method.

######  Assumption 3.1.

The function $f(t,\zeta)$ grows mildly with respect to $\zeta$, i.e., there
exists a number $p>0$ for $d=1,2$ or $p\in(0,2]$ for $d=3$ such that

$\Big{|}\frac{\partial f}{\partial\zeta}(t,\zeta)\Big{|}\lesssim
1+|\zeta|^{p},\quad\forall\,t,\,\zeta\in\mathbb{R}.$ (3.11)

######  Assumption 3.2.

The function $f(t,\zeta)$ is sufficiently smooth with respect to $t$, i.e.,
for any given constant $K>0$, it holds

$\textstyle\sum\limits_{|\bm{\alpha}|\leq
2}\Big{|}D^{\bm{\alpha}}f(t,\zeta)\Big{|}\lesssim
1,\quad\forall\,t\in[0,T],\zeta\in[-K,K].$ (3.12)

######  Assumption 3.3.

The exact solution $u(t)$ satisfies some of the following regularity
conditions:

$\displaystyle\textstyle\sup\limits_{0\leq t\leq T}\|u(t)\|_{2,\Omega}$
$\displaystyle\lesssim 1,$ (3.13a) $\displaystyle\textstyle\sup\limits_{0\leq
t\leq T}\|u_{t}(t)\|_{0,\infty,\Omega}$ $\displaystyle\lesssim 1,$ (3.13b)
$\displaystyle\textstyle\sup\limits_{0\leq t\leq
T}\|u_{tt}(t)\|_{0,\infty,\Omega}$ $\displaystyle\lesssim 1,$ (3.13c) where
the hidden constants may depend on $T$.

We then have the following result on the locally-Lipschitz continuity of $f$.

######  Lemma 3.2.

Suppose that the function $f$ satisfies Assumption 3.1, and the exact solution
$u(t)$ fulfills (3.13a) in Assumptions 3.3. Then $f$ is locally-Lipschitz
continuous in a strip along the exact solution $u(t)$, i.e., for any given
constant $R>0$,

$\|f(t,v)-f(t,w)\|_{0}\lesssim\|v-w\|_{1},$ (3.14)

for any $t\in[0,T]$ and $v,w\in V_{h}$ satisfying

$\textstyle\max\\{\|v-u(t)\|_{1},\|w-u(t)\|_{1}\\}\leq R.$

###### Proof.

It follows from (3.13a) and (3.2) that

$\textstyle\max\\{\|v\|_{1},\|w\|_{1}\\}\leq R+\|u(t)\|_{1}\lesssim 1+R.$
(3.15)

On the other hand, by the Lagrange mean value theorem and (3.11),

$\begin{split}\|f(t,v)-f(t,w)\|_{0}^{2}&=\Big{\|}\frac{\partial f}{\partial
u}(t,\xi)(v-w)\Big{\|}_{0}^{2}\leq\int_{\Omega}(1+|\xi|^{p})^{2}|v-w|^{2}\,{\rm
d}\bm{x}\\\ &\leq\int_{\Omega}(1+|v|^{p})^{2}|v-w|^{2}\,{\rm
d}\bm{x}+\int_{\Omega}(1+|w|^{p})^{2}|v-w|^{2}\,{\rm d}\bm{x},\end{split}$
(3.16)

where $\xi(\bm{x})=\theta(\bm{x})v(\bm{x})+(1-\theta(\bm{x}))w(\bm{x})$ for
some $\theta(\bm{x})\in[0,1]$. It is clear that

$\int_{\Omega}(1+|v|^{p})^{2}|v-w|^{2}\,{\rm
d}\bm{x}=\int_{\Omega}|v|^{2p}|v-w|^{2}\,{\rm
d}\bm{x}+2\int_{\Omega}|v|^{p}|v-w|^{2}\,{\rm
d}\bm{x}+\int_{\Omega}|v-w|^{2}\,{\rm d}\bm{x}.$ (3.17)

It suffices to show below the bound of the first term in the right-hand side
of (3.17) since the other two terms can be treated similarly.

Case I: $d=1$ or $2$. In this case, $0<p<\infty$. We choose $q_{1},q_{2}$
satisfying

$\frac{1}{q_{1}}+\frac{1}{q_{2}}=1,\quad 1\leq q_{1}<\infty,\quad 1\leq
q_{2}<\infty,\quad 2pq_{1}\geq 1.$

Then by using Holder’s inequality, Sobolev embedding theorem and (3.15),

$\begin{split}\int_{\Omega}|v|^{2p}|v-w|^{2}\,{\rm
d}x&\leq\Big{(}\int_{\Omega}|v|^{2pq_{1}}\,{\rm
d}\bm{x}\Big{)}^{\frac{1}{q_{1}}}\Big{(}\int_{\Omega}|v-w|^{2q_{2}}\,{\rm
d}\bm{x}\Big{)}^{\frac{1}{q_{2}}}\\\
&=\|v\|_{L^{2pq_{1}}}^{2p}\|v-w\|_{L^{2q_{2}}}^{2}\lesssim\|v\|_{1}^{2p}\|v-w\|_{1}^{2}\lesssim\|v-w\|_{1}^{2}.\end{split}$
(3.18)

Case II: $d=3$. In this case, $0<p\leq 2$. In view of the Sobolev embedding
theorem, $H^{1}(\Omega)\hookrightarrow L^{q}(\Omega)$, where $1\leq q\leq 6$.
We select appropriate $q_{1},q_{2}$ satisfying

$\frac{1}{q_{1}}+\frac{1}{q_{2}}=1,\quad 1\leq q_{1}<\infty,\quad 1\leq
q_{2}\leq 3,\quad 1\leq 2pq_{1}\leq 6.$

In fact, these conditions hold if $(q_{1},q_{2})$ satisfies the conditions

$\textstyle\max\Big{\\{}\frac{3}{2},\frac{1}{2p}\Big{\\}}\leq
q_{1}\leq\frac{3}{p},\quad 1\leq q_{2}\leq 3.$

Noting $0<p\leq 2$, we can find the existence of such a pair $(q_{1},q_{2})$
in terms of the above conditions. Follow the similar arguments for deriving
(3.18), we have

$\begin{split}\int_{\Omega}|v|^{2p}|v-w|^{2}\,{\rm
d}\bm{x}&\leq\Big{(}\int_{\Omega}|v-w|^{2q_{2}}\,{\rm
d}\bm{x}\Big{)}^{\frac{1}{q_{2}}}\Big{(}\int_{\Omega}|v|^{2pq_{1}}\,{\rm
d}\bm{x}\Big{)}^{\frac{1}{q_{1}}}\\\
&=\|v-w\|_{L^{2q_{2}}}^{2}\|v\|_{L^{2pq_{1}}}^{2p}\lesssim\|v-w\|_{1}^{2}\|v\|_{1}^{2p}\lesssim\|v-w\|_{1}^{2}.\end{split}$
(3.19)

Based on (3.18), (3.19) and (3.17), we obtain that

$\int_{\Omega}(1+|v|^{p})^{2}|v-w|^{2}\,{\rm d}\bm{x}\lesssim\|v-w\|_{1}^{2},$
(3.20)

and similarly

$\int_{\Omega}(1+|w|^{p})^{2}|v-w|^{2}\,{\rm d}\bm{x}\lesssim\|v-w\|_{1}^{2},$
(3.21)

Finally, the combination of (3.16), (3.20) and (3.21) leads to (3.14). ∎

###  3.2 Fully-discrete error estimates

In the rest of this section, for simplicity we will assume the time partition
is uniform, i.e., $\Delta t=\tau_{0}=\cdots=\tau_{N_{T}-1}$ and $t_{n}=n\Delta
t$. Let $u_{h}(t)$ be the solution of the semi-discrete (in space) problem
(2.8) (or (2.9)), and $\\{u_{h}^{n}\\}$ the fully-discrete solution produced
by the EIFE method (2.12). For the error between the exact solution $u(t)$ and
the fully-discrete solution $\\{u_{h}^{n}\\}$ measured in the $H^{1}$-norm, we
have by the triangle inequality

$\|u(t_{n})-u_{h}^{n}\|_{1}\leq\|u(t_{n})-u_{h}(t_{n})\|_{1}+\|u_{h}(t_{n})-u_{h}^{n}\|_{1}.$
(3.22)

Thus we will estimate $\|u(t_{n})-u_{h}(t_{n})\|_{1}$ and
$\|u_{h}(t_{n})-u_{h}^{n}\|_{1}$ respectively.

First of all, the following result readily comes from Theorem 14.2 in [40].

######  Lemma 3.3.

Suppose that the function $f$ satisfies Assumptions 3.1 and 3.2, and the exact
solution $u(t)$ fulfills (3.13a) in Assumptions 3.3. There exists a constant
$h_{0}>0$ such that if the spatial mesh size $h\leq h_{0}$, then

$\|u_{h}(t)-u(t)\|_{1}\lesssim h,\quad\forall\,t\in[0,T],$ (3.23)

where the hidden constant is independent of $h$.

A direct consequence of Lemma 3.3 is

$\|u(t_{n})-u_{h}(t_{n})\|_{1}\lesssim h,\quad\forall\,n=0,\cdots,N_{T},$
(3.24)

where the hidden constant is independent of $h$. Then the remaining part for
us is to bound $\|u_{h}(t_{n})-u_{h}^{n}\|_{1}$. Note that the semi-discrete
solution $u_{h}(t)$ obtained from (2.8) itself depends on the spatial mesh
size $h$, on the other hand we need remove the dependence of the hidden
constants on $h$ in estimating $\|u_{h}(t_{n})-u_{h}^{n}\|_{1}$. In order to
achieve this goal, we convert the semi-discrete solution $u_{h}(t_{n+1})$
($n=0,1,\cdots,N_{T}-1$) into the sum of the following two parts for further
analysis:

$\begin{split}u_{h}(t_{n+1})=\;&e^{-\Delta
tL_{h}}u_{h}(t_{n})+\int_{0}^{\Delta t}e^{-(\Delta
t-\sigma)L_{h}}P_{h}f(t_{n}+\sigma,u_{h}(t_{n}+\sigma))\,{\rm d}\sigma\\\
=\;&e^{-\Delta tL_{h}}u_{h}(t_{n})+\int_{0}^{\Delta t}e^{-(\Delta
t-\sigma)L_{h}}P_{h}f(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm d}\sigma\\\
&+\int_{0}^{\Delta t}e^{-(\Delta
t-\sigma)L_{h}}\big{(}P_{h}f(t_{n}+\sigma,u_{h}(t_{n}+\sigma))-P_{h}f(t_{n}+\sigma,u(t_{n}+\sigma))\big{)}\,{\rm
d}\sigma.\end{split}$ (3.25)

Define the following functions:

$\left\\{\begin{split}\psi_{i}(-\Delta tL_{h})&=\varphi_{i}(-\Delta
tL_{h})-\textstyle\sum\limits_{k=1}^{s}b_{k}(-\Delta
tL_{h})\frac{c_{k}^{i-1}}{(i-1)!},\;i=1,\cdots,s,\\\ \psi_{j,i}(-\Delta
tL_{h})&=\varphi_{j}(-c_{i}\Delta
tL_{h})c_{i}^{j}-\textstyle\sum\limits_{k=1}^{i-1}a_{ik}(-\Delta
tL_{h})\frac{c_{k}^{j-1}}{(j-1)!},\;i,j=1,\cdots,s.\end{split}\right.$ (3.26)

We also denote $f^{(k)}(t,u(t))=\frac{d^{k}}{dt^{k}}f(t,u(t))$ as the $k$-th
full differentiation of $f$ with respect to $t$. By comparing (3.25) with the
fully-discrete scheme (2.12), we then obtain

$\displaystyle u_{h}(t_{n}+c_{i}\Delta t)$ $\displaystyle=$ $\displaystyle
e^{-c_{i}\Delta tL_{h}}u_{h}(t_{n})+\Delta
t\textstyle\sum\limits_{j=1}^{i-1}a_{ij}(-\Delta tL_{h})$ (3.27)
$\displaystyle\qquad P_{h}f(t_{n}+c_{j}\Delta t,u(t_{n}+c_{j}\Delta
t))+\delta_{ni},$ $\displaystyle u_{h}(t_{n+1})$ $\displaystyle=$
$\displaystyle e^{-\Delta tL_{h}}u_{h}(t_{n})+\Delta
t\textstyle\sum\limits_{i=1}^{s}b_{i}(-\Delta tL_{h})$ (3.28)
$\displaystyle\qquad P_{h}f(t_{n}+c_{i}\Delta t,u(t_{n}+c_{i}\Delta
t))+\delta_{n+1},$

where the defect terms $\\{\delta_{ni}\\}_{i=1}^{s}$ and $\delta_{n+1}$ are
respectively given by

$\begin{split}\delta_{ni}&=\textstyle\sum\limits_{j=1}^{r}\Delta
t^{j}\psi_{j,i}(-\Delta
tL_{h})P_{h}f^{(j-1)}(t_{n},u(t_{n}))+\delta_{ni}^{[r]},\\\
\delta_{n+1}&=\textstyle\sum\limits_{i=1}^{r}\Delta t^{i}\psi_{i}(-\Delta
tL_{h})P_{h}f^{(i-1)}(t_{n},u(t_{n}))+\delta_{n+1}^{[r]},\\\ \end{split}$

with the remainders $\delta_{ni}^{[r]}$ and $\delta_{n+1}^{[r]}$ defined
respectively by

$\begin{split}\delta_{ni}^{[r]}=&\int_{0}^{c_{i}\Delta t}e^{-(c_{i}\Delta
t-\tau)L_{h}}\int_{0}^{\tau}\frac{(\tau-\sigma)^{r-1}}{(r-1)!}P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm
d}\sigma\,{\rm d}\tau\\\ &-\Delta
t\textstyle\sum\limits_{k=1}^{i-1}a_{ik}(-\Delta
tL_{h})\displaystyle\int_{0}^{c_{k}\Delta t}\frac{(c_{k}\Delta
t-\sigma)^{r-1}}{(r-1)!}P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm
d}\sigma\\\ &+\int_{0}^{c_{i}\Delta t}e^{-(c_{i}\Delta
t-\sigma)L_{h}}\big{(}P_{h}f(t_{n}+\sigma,u_{h}(t_{n}+\sigma))-P_{h}f(t_{n}+\sigma,u(t_{n}+\sigma))\big{)}\,{\rm
d}\sigma,\end{split}$ $\begin{split}\delta_{n+1}^{[r]}=&\int_{0}^{\Delta
t}e^{-(\Delta
t-\tau)L_{h}}\int_{0}^{\tau}\frac{(\tau-\sigma)^{r-1}}{(r-1)!}P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm
d}\sigma\,{\rm d}\tau\\\ &-\Delta
t\textstyle\sum\limits_{i=1}^{s}b_{i}(-\Delta
tL_{h})\displaystyle\int_{0}^{c_{i}\Delta t}\frac{(c_{i}\Delta
t-\sigma)^{r-1}}{(r-1)!}P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm
d}\sigma\\\ &+\int_{0}^{\Delta t}e^{-(\Delta
t-\sigma)L_{h}}\big{(}P_{h}f(t_{n}+\sigma,u_{h}(t_{n}+\sigma))-P_{h}f(t_{n}+\sigma,u(t_{n}+\sigma))\big{)}\,{\rm
d}\sigma.\end{split}$

Here $r$ can be any nonnegative integers such that $f^{(r)}(t,u(t))$ exists
and is continuous.

In what follows, we will adopt the arguments proposed in [18] to bound
$\|u_{h}^{n}-u_{h}(t_{n})\|_{1}$. For brevity, let us define
$e_{n}=u_{h}^{n}-u_{h}(t_{n})$ and $E_{ni}=U_{ni}-u_{h}(t_{n}+c_{i}\Delta t)$
for $i=1,\cdots,s$. Then we arrive at the following recurrence relations:

$\displaystyle E_{ni}=e^{-c_{i}\Delta tL_{h}}e_{n}+\Delta
t\textstyle\sum\limits_{j=1}^{i-1}a_{ij}(-\Delta
tL_{h})\big{(}P_{h}f(t_{n}+c_{j}\Delta t,U_{nj})$ $\displaystyle\hskip
142.26378pt-P_{h}f(t_{n}+c_{j}\Delta t,u_{h}(t_{n}+c_{j}\Delta
t))\big{)}-\delta_{ni}.$ (3.29) $\displaystyle e_{n+1}=e^{-\Delta
tL_{h}}e_{n}+\Delta t\textstyle\sum\limits_{i=1}^{s}b_{i}(-\Delta
tL_{h})\big{(}P_{h}f(t_{n}+c_{i}\Delta t,U_{ni})$ $\displaystyle\hskip
142.26378pt-P_{h}f(t_{n}+c_{i}\Delta t,u_{h}(t_{n}+c_{i}\Delta
t))\big{)}-\delta_{n+1},$ (3.30)

We first have the following result on the defect terms in (3.27) and (3.28).

######  Lemma 3.4.

Given an integer $r=1$ or $2$. Suppose that the function $f$ satisfies
Assumptions 3.1 and 3.2, and the exact solution $u(t)$ fulfills (3.13a) and
(3.13b) in Assumptions 3.3. Suppose that $u(t)$ additionally fulfills (3.13c)
if $r=2$. Then for $n=0,\cdots,N_{T}$, $i=1,\cdots,s$, it holds that

$\displaystyle\textstyle\|\delta_{ni}^{[r]}\|_{1}\lesssim(\Delta
t)^{r+1}\textstyle\sup\limits_{0\leq\eta\leq 1}\|f^{(r)}(t_{n}+\eta\Delta
t,u(t_{n}+\eta\Delta t))\|_{1}+h,$ (3.31a)
$\displaystyle\textstyle\Big{\|}\textstyle\sum\limits_{j=0}^{n-1}e^{-j\Delta
tL_{h}}\delta_{n-j}^{[r]}\Big{\|}_{1}\lesssim(\Delta
t)^{r}\textstyle\sup\limits_{0\leq t\leq T}\|f^{(r)}(t,u(t))\|_{1}+h.$ (3.31b)
Note that the above hidden constants are independent of $h$ and $\Delta t$.

###### Proof.

Recalling (3.1) in Lemma 3.1 and the relation (3.10), we have after some
direct manipulations

$\begin{split}&\Big{\|}\int_{0}^{c_{i}\Delta t}e^{-(c_{i}\Delta
t-\tau)L_{h}}\int_{0}^{\tau}\frac{(\tau-\sigma)^{r-1}}{(r-1)!}P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm
d}\sigma\,{\rm d}\tau\Big{\|}_{1}\\\ &\;\lesssim\Big{\|}\int_{0}^{c_{i}\Delta
t}e^{-(c_{i}\Delta
t-\tau)L_{h}}\int_{0}^{\tau}\frac{(\tau-\sigma)^{r-1}}{(r-1)!}L_{h}^{\frac{1}{2}}P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm
d}\sigma\,{\rm d}\tau\Big{\|}_{0}\\\ &\;\lesssim\textstyle(\Delta
t)^{r+1}\sup\limits_{0\leq\tau\leq c_{i}\Delta t}\|e^{-(c_{i}\Delta
t-\tau)L_{h}}\|_{0}\textstyle\sup\limits_{0\leq\eta\leq
1}\|L_{h}^{\frac{1}{2}}P_{h}f^{(r)}(t_{n}+\eta\Delta t,u(t_{n}+\eta\Delta
t))\|_{0}\\\ &\;\lesssim(\Delta t)^{r+1}\textstyle\sup\limits_{0\leq\eta\leq
1}\|P_{h}f^{(r)}(t_{n}+\eta\Delta t,u(t_{n}+\eta\Delta t))\|_{1}\\\
&\;\lesssim(\Delta t)^{r+1}\textstyle\sup\limits_{0\leq\eta\leq
1}\|f^{(r)}(t_{n}+\eta\Delta t,u(t_{n}+\eta\Delta t))\|_{1},\end{split}$
(3.32)

where the last inequality is due to that $P_{h}$ is $H^{1}$-stable. Similarly,
we also have

$\begin{split}&\Big{\|}\int_{0}^{\Delta t}e^{-(\Delta
t-\sigma)L_{h}}\Big{(}P_{h}f(t_{n}+\sigma,u_{h}(t_{n}+\sigma))-P_{h}f(t_{n}+\sigma,u(t_{n}+\sigma))\Big{)}\,{\rm
d}\sigma\Big{\|}_{1}\\\ &\lesssim\Big{\|}\int_{0}^{\Delta
t}L_{h}^{\frac{1}{2}}e^{-(\Delta t-\sigma)L_{h}}\,{\rm
d}\sigma\Big{\|}_{0}\textstyle\sup\limits_{0\leq\eta\leq
1}\|P_{h}f(t_{n}+\eta\Delta t,u_{h}(t_{n}+\eta\Delta t))\\\ &\quad-
P_{h}f(t_{n}+\eta\Delta t,u(t_{n}+\eta\Delta t))\Big{\|}_{0}.\end{split}$
(3.33)

In view of the similar arguments for proving Lemma 3.1, we have

$\begin{split}\Big{\|}\int_{0}^{\Delta t}L_{h}^{\frac{1}{2}}e^{-(\Delta
t-\sigma)L_{h}}\,{\rm
d}\sigma\Big{\|}_{0}&=\textstyle\max\limits_{\lambda\in\lambda(L_{h})}\Big{|}\int_{0}^{\Delta
t}\lambda^{\frac{1}{2}}e^{-(\Delta t-\sigma)\lambda}\,{\rm d}\sigma\Big{|}\\\
&\leq\textstyle\max\limits_{\lambda\in\lambda(L_{h})}|\lambda^{-\frac{1}{2}}|\lesssim
1.\end{split}$ (3.34)

Since $P_{h}$ is $L^{2}$ stable and the function $f$ is locally-Lipschitz
continuous (Lemma 3.2), we further obtain from (3.33)-(3.34) and Lemma 3.3
that

$\begin{split}&\Big{\|}\int_{0}^{\Delta t}e^{-(\Delta
t-\sigma)L_{h}}\big{(}P_{h}f(t_{n}+\sigma,u_{h}(t_{n}+\sigma))-P_{h}f(t_{n}+\sigma,u(t_{n}+\sigma))\big{)}\,{\rm
d}\sigma\Big{\|}_{1}\\\ &\lesssim\textstyle\sup\limits_{0\leq\eta\leq
1}\|P_{h}f(t_{n}+\eta\Delta t,u_{h}(t_{n}+\eta\Delta
t))-P_{h}f(t_{n}+\eta\Delta t,u(t_{n}+\eta\Delta t))\|_{0}\\\
&\lesssim\textstyle\sup\limits_{0\leq\eta\leq 1}\|u_{h}(t_{n}+\eta\Delta
t)-u(t_{n}+\eta\Delta t)\|_{1}\lesssim h.\end{split}$ (3.35)

According to (3.3) in Lemma 3.1 and the similar arguments for deriving (3.32),

$\begin{split}\quad&\Big{\|}\Delta
t\textstyle\sum\limits_{k=1}^{i-1}a_{ik}(-\Delta
tL_{h})\displaystyle\int_{0}^{c_{k}\Delta t}\frac{(c_{k}\Delta
t-\sigma)^{r-1}}{(r-1)!}P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm
d}\sigma\Big{\|}_{1}\\\ &\lesssim\Big{\|}\Delta
t\textstyle\sum\limits_{k=1}^{i-1}a_{ik}(-\Delta
tL_{h})\displaystyle\int_{0}^{c_{k}\Delta t}\frac{(c_{k}\Delta
t-\sigma)^{r-1}}{(r-1)!}L_{h}^{\frac{1}{2}}\\\ &\qquad\quad\cdot
P_{h}f^{(r)}(t_{n}+\sigma,u(t_{n}+\sigma))\,{\rm d}\sigma\Big{\|}_{0}\\\
&\leq(\Delta t)^{r+1}\textstyle\sum\limits_{k=1}^{i-1}\|a_{ik}(-\Delta
tL_{h})\|_{0}\sup\limits_{0\leq\eta\leq 1}\|P_{h}f^{(r)}(t_{n}+\eta\Delta
t,u(t_{n}+\eta\Delta t))\|_{1}\\\ &\lesssim(\Delta
t)^{r+1}\textstyle\sup\limits_{0\leq\eta\leq 1}\|f^{(r)}(t_{n}+\eta\Delta
t,u(t_{n}+\eta\Delta t))\|_{1}.\end{split}$ (3.36)

Now, using the triangle inequality, the regularity assumptions for $u(t)$ and
$f$, and the estimates (3.32), (3.35) and (3.36), we get

$\|\delta_{ni}^{[r]}\|_{1}\lesssim(\Delta t)^{r+1}\|f^{(r)}(t_{n}+\eta\Delta
t,u(t_{n}+\eta\Delta t))\|_{1}+h,\quad\forall\,i=1,\cdots,s,$ (3.37)

which leads to (3.31a). Also (3.31b) can be derived in the similar manner. ∎

For the EIFE1 scheme (2.13) (i.e., $s=1$), the following error equation holds:

$e_{n+1}=e^{-\Delta tL_{h}}e_{n}+\Delta t\varphi_{1}(-\Delta
tL_{h})\big{(}P_{h}f(t_{n},u_{h}^{n})-P_{h}f(t_{n},u(t_{n}))\big{)}-\delta_{n+1}.$
(3.38)

######  Theorem 3.1 (Error estimate for the EIFE1 scheme).

Suppose the function $f$ satisfies Assumptions 3.1 and 3.2, and the exact
solution $u(t)$ fulfills (3.13a) and (3.13b). There exists a constant
$h_{0}>0$ such that if the spatial mesh size $h\leq h_{0}$, then the numerical
solution $\\{u_{h}^{n}\\}$ produced by the EIFE1 scheme (2.13) satisfies

$\|u(t_{n})-u_{h}^{n}\|_{1}\lesssim\Delta t+h,\quad\forall\,n=1,\cdots,N_{T},$
(3.39)

where the hidden constant is independent of $h$ and $\Delta t$.

###### Proof.

By repeatedly applying (3.38), we get

$\begin{split}e_{n}=\;&\Delta
t\textstyle\sum\limits_{j=1}^{n-1}e^{-(n-j-1)\Delta tL_{h}}\varphi_{1}(-\Delta
tL_{h})\big{(}P_{h}f(t_{j},u_{h}^{j})-P_{h}f(t_{j},u(t_{j}))\big{)}\\\
&-\textstyle\sum\limits_{j=0}^{n-1}e^{-j\Delta
tL_{h}}\delta_{n-j}.\end{split}$ (3.40)

Since $s=1$, based on the definition of $\psi_{i}(-\Delta tL_{h})$ in (3.26)
and Lagrangian interpolation theorem, we can derive that $\psi_{1}(-\Delta
tL_{h})=0$, which implies $\delta_{j}=\delta_{j}^{[1]}$ for $j=1,\cdots,n$.
Recalling (3.31b) in Lemma 3.4 (with $r=1$), we can immediately obtain the
$H^{1}$-norm estimation of the second term in (3.40) as

$\displaystyle\Big{\|}\textstyle\sum\limits_{j=0}^{n-1}e^{-j\Delta
tL_{h}}\delta_{n-j}\Big{\|}_{1}$
$\displaystyle=\Big{\|}\textstyle\sum\limits_{j=0}^{n-1}e^{-j\Delta
tL_{h}}\delta_{n-j}^{[1]}\Big{\|}_{1}\lesssim\Delta
t\textstyle\sup\limits_{0\leq t\leq T}\|f^{(1)}(t,u(t))\|_{1}+h.$ (3.41)

As for the $H^{1}$-norm evaluation of the first term in (3.40), it follows
from (3.10) and the triangle inequality that

$\begin{split}&\Big{\|}\Delta
t\textstyle\sum\limits_{j=1}^{n-1}e^{-(n-j-1)\Delta tL_{h}}\varphi_{1}(-\Delta
tL_{h})\big{(}P_{h}f(t_{j},u_{h}^{j})-P_{h}f(t_{j},u(t_{j}))\big{)}\Big{\|}_{1}\\\
&\lesssim\|\Delta tL_{h}^{\frac{1}{2}}\varphi_{1}(-\Delta
tL_{h})\|_{0}\|P_{h}f(t_{n-1},u_{h}^{n-1})-P_{h}f(t_{n-1},u(t_{n-1}))\|_{0}\\\
&\quad+\Big{\|}\textstyle\sum\limits_{j=0}^{n-2}\Delta
tL_{h}^{\frac{1}{2}}e^{-(n-j-1)\Delta tL_{h}}\varphi_{1}(-\Delta
tL_{h})\big{(}P_{h}f(t_{j},u_{h}^{j})-P_{h}f(t_{j},u(t_{j}))\big{)}\Big{\|}\\\
&:={\rm I}_{1}+{\rm I}_{2}.\end{split}$ (3.42)

By applying (3.1) in Lemma 3.1, we get

$\begin{split}{\rm I}_{1}=\,&\Big{\|}L_{h}^{\frac{1}{2}}\int_{0}^{\Delta
t}e^{-(\Delta t-\sigma)L_{h}}\,{\rm
d}\sigma\Big{\|}_{0}\|P_{h}f(t_{n-1},u_{h}^{n-1})-P_{h}f(t_{n-1},u(t_{n-1}))\|_{0}\\\
\lesssim\,&\Delta t\textstyle\sup\limits_{0\leq\sigma\leq\Delta
t}\|L_{h}^{\frac{1}{2}}e^{-(\Delta
t-\sigma)L_{h}}\|_{0}\|P_{h}f(t_{n-1},u_{h}^{n-1})-P_{h}f(t_{n-1},u(t_{n-1}))\|_{0}\\\
\lesssim\,&(\Delta t)^{\frac{1}{2}}\|u_{h}^{n-1}-u(t_{n-1})\|_{1}.\end{split}$
(3.43)

For the second term in (3.42), it follows from (3.1) in Lemma 3.1 that

$\|\varphi_{1}(-\Delta tL_{h})\|_{0}=\Big{\|}\frac{1}{\Delta
t}\int_{0}^{\Delta t}e^{-(\Delta t-\tau)L_{h}}\,{\rm
d}\tau\Big{\|}_{0}\leq\textstyle\sup\limits_{0\leq\tau\leq\Delta
t}\|e^{-(\Delta t-\tau)}\|_{0}\lesssim 1.$

By (3.1) and (3.2) in Lemma 3.1 and similar arguments for deriving (3.35), we
have

$\begin{split}{\rm I}_{2}\lesssim\,&\Big{\|}\Delta
tL_{h}^{\frac{1}{2}}\textstyle\sum\limits_{j=0}^{n-2}e^{-(n-1-j)\Delta
tL_{h}}\Big{\|}_{0}\textstyle\sup\limits_{0\leq t\leq
T}\|P_{h}f(t,u_{h}(t))-P_{h}f(t,u(t))\|_{0}\\\ &+\Delta
t\textstyle\sum\limits_{j=0}^{n-2}\|L_{h}^{\frac{1}{2}}e^{-(n-1-j)\Delta
tL_{h}}\|_{0}\|P_{h}f(t_{j},u_{h}^{j})-P_{h}f(t_{j},u_{h}(t_{j}))\|_{0}\\\
\lesssim\,&\textstyle\sup\limits_{0\leq t\leq
T}\|P_{h}f(t,u_{h}(t))-P_{h}f(t,u(t))\|_{0}+\Delta
t\textstyle\sum\limits_{j=0}^{n-2}t_{n-j-1}^{-\frac{1}{2}}\|P_{h}f(t_{j},u_{h}^{j})\\\
&\qquad-P_{h}f(t_{j},u_{h}(t_{j}))\|_{0}\\\ \lesssim\,&\Delta
t\textstyle\sum\limits_{j=0}^{n-2}t_{n-j-1}^{-\frac{1}{2}}\|e_{j}\|_{1}+h.\end{split}$
(3.44)

In terms of (3.40) and the estimates (3.41)-(3.44), we arrive at

$\begin{split}\|e_{n}\|_{1}&\lesssim(\Delta
t)^{\frac{1}{2}}\|e_{n-1}\|_{1}+\Delta
t\textstyle\sum\limits_{j=0}^{n-2}t_{n-j-1}^{-\frac{1}{2}}\|e_{j}\|_{1}+\Delta
t+h\\\ &\lesssim\Delta
t\textstyle\sum\limits_{j=1}^{n-1}t_{n-j}^{-\frac{1}{2}}\|e_{j}\|_{1}+\Delta
t+h.\end{split}$

Then we have by the discrete Gronwall inequality (Theorem 6.1 with
$\alpha=\frac{1}{2}$ in [9])

$\|u_{h}(t_{n})-u_{h}^{n}\|_{1}\lesssim\Delta t+h.$ (3.45)

The combination of (3.22), (3.24) and (3.45) immediately leads to (3.39). ∎

######  Lemma 3.5.

Suppose the function $f$ satisfies Assumptions 3.1 and 3.2, and the exact
solution $u(t)$ fulfills (3.13a)-(3.13b). If $s\geq 2$, then it holds for any
$0\leq n<N_{T}$,

$\|E_{ni}\|_{1}\lesssim\|e_{n}\|_{1}+(\Delta
t)^{2}\textstyle\sup\limits_{0\leq\eta\leq 1}\|f^{\prime}(t_{n}+\eta\Delta
t,u(t_{n}+\eta\Delta t))\|_{1}+h,\quad\forall\,i=1,\cdots,s,$ (3.46)

where the hidden constant is independent of $h$ and $\Delta t$.

###### Proof.

According to the definition of $\psi_{j,i}$ in (3.26), we have by some
manipulations that $\psi_{1,j}=0,j=1,\cdots,s$ when the consistency conditions
(2.11) is fulfilled. Therefore, the estimation of $\|\delta_{ni}\|_{1}$ can be
converted to that of $\|\delta_{ni}^{[1]}\|_{1}$. Using the similar arguments
for deriving the estimate (3.35) and (3.3) in Lemma 3.1, we have

$\displaystyle\begin{split}&\Big{\|}\Delta
t\textstyle\sum\limits_{j=1}^{i-1}a_{ij}(-\Delta
tL_{h})\big{(}P_{h}f(t_{n}+c_{j}\Delta t,U_{nj})-P_{h}f(t_{n}+c_{j}\Delta
t,u(t_{n}+c_{j}\Delta t))\big{)}\Big{\|}_{1}\\\
&\lesssim\textstyle\sum\limits_{j=1}^{i-1}(\Delta t)^{\frac{1}{2}}\|\Delta
t^{\frac{1}{2}}L_{h}^{\frac{1}{2}}a_{ij}(-\Delta tL_{h})\|_{0}\\\
&\qquad\cdot\textstyle\max\limits_{2\leq j\leq i-1}\|P_{h}f(t_{n}+c_{j}\Delta
t,U_{nj})-P_{h}f(t_{n}+c_{j}\Delta t,u(t_{n}+c_{j}\Delta t))\|_{0}\\\
&\lesssim(\Delta t)^{\frac{1}{2}}\textstyle\max\limits_{2\leq j\leq
i-1}\|P_{h}f(t_{n}+c_{j}\Delta t,U_{nj})-P_{h}f(t_{n}+c_{j}\Delta
t,u(t_{n}+c_{j}\Delta t))\|_{0}\\\ &\lesssim(\Delta
t)^{\frac{1}{2}}\textstyle\max\limits_{2\leq j\leq
i-1}\|E_{nj}\|_{1}+h.\end{split}$

Note that $\|\delta_{ni}^{[1]}\|_{1}$ is uniformly bounded for $i=1,\cdots,s$
(see (3.31a) in Lemma 3.4 with $r=1$). Recalling the relation (3.29), we have
by the triangle inequality that

$\begin{split}\|E_{ni}\|_{1}\lesssim&\;\|e^{-c_{i}\Delta
tL_{h}}e_{n}\|_{1}+(\Delta t)^{2}\textstyle\sup\limits_{0\leq\eta\leq
1}\|P_{h}f^{\prime}(t_{n}+\eta\Delta t,u(t_{n}+\eta\Delta t))\|_{1}+h\\\
&+\Big{\|}\Delta t\textstyle\sum\limits_{j=1}^{i-1}a_{ij}(-\Delta
tL_{h})\big{(}P_{h}f(t_{n}+c_{j}\Delta t,U_{nj})\\\ &\qquad-
P_{h}f(t_{n}+c_{j}\Delta t,u(t_{n}+c_{j}\Delta t))\big{)}\Big{\|}_{1}\\\
\lesssim&\;\|e_{n}\|_{1}+(\Delta t)^{2}\textstyle\sup\limits_{0\leq\eta\leq
1}\|P_{h}f^{\prime}(t_{n}+\eta\Delta t,u(t_{n}+\eta\Delta t))\|_{1}+h\\\
&+(\Delta t)^{\frac{1}{2}}\textstyle\max\limits_{2\leq j\leq
i-1}\|E_{nj}\|_{0}\\\ \lesssim&\;\|e_{n}\|_{1}+(\Delta
t)^{2}\textstyle\sup\limits_{0\leq\eta\leq 1}\|f^{\prime}(t_{n}+\eta\Delta
t,u(t_{n}+\eta\Delta t))\|_{1}+h\\\ &+(\Delta
t)^{\frac{1}{2}}\textstyle\max\limits_{2\leq j\leq
i-1}\|E_{nj}\|_{1}.\end{split}$

Finally (3.46) is obtained by recursively using the above inequality. ∎

######  Theorem 3.2 (Error estimate for the EIFE2 scheme).

Suppose that the function $f$ satisfies Assumptions 3.1 and 3.2, and the exact
solution $u(t)$ fulfills (3.13a)-(3.13c) in Assumptions 3.3. There exists a
constant $h_{0}>0$ such that if the spatial mesh size $h\leq h_{0}$, then the
numerical solution $\\{u_{h}^{n}\\}$ produced by EIFE2 scheme (2.14) satisfies

$\|u(t_{n})-u_{h}^{n}\|_{1}\lesssim(\Delta
t)^{2}+h,\quad\forall\,n=1,\cdots,N_{T},$ (3.47)

where the hidden constant is independent of $h$ and $\Delta t$.

###### Proof.

Recalling the definition of $\psi_{i}$ in (3.26), we can check that
$\psi_{1}(-\Delta tL_{h})=\psi_{2}(-\Delta tL_{h})=0$, which implies that
$\delta_{n+1}=\delta_{n+1}^{[2]}$ by Lagrangian interpolation theorem for
$s=2$. By (3.3) in Lemma 3.1, we have

$\begin{split}&\Big{\|}\Delta
t\textstyle\sum\limits_{j=0}^{n-1}e^{-(n-1-j)\Delta
tL_{h}}\textstyle\sum\limits_{i=1}^{s}b_{i}(-\Delta
tL_{h})\big{(}P_{h}f(t_{j}+c_{i}\Delta t,U_{ji})\\\ &\qquad-
P_{h}f(t_{j}+c_{i}\Delta t,u(t_{j}+c_{i}\Delta t))\big{)}\Big{\|}_{1}\\\
&\lesssim\Big{\|}\Delta t\textstyle\sum\limits_{i=1}^{s}b_{i}(-\Delta
tL_{h})\big{(}P_{h}f(t_{n-1}+c_{i}\Delta t,U_{n-1,i})\\\ &\qquad-
P_{h}f(t_{n-1}+c_{i}\Delta t,u(t_{n-1}+c_{i}\Delta t)\big{)}\Big{\|}_{1}\\\
&\qquad+\Big{\|}\Delta t\textstyle\sum\limits_{j=0}^{n-2}e^{-(n-1-j)\Delta
tL_{h}}\textstyle\sum\limits_{i=1}^{s}\big{(}P_{h}f(t_{j}+c_{i}\Delta
t,U_{ji})\\\ &\qquad-P_{h}f(t_{j}+c_{i}\Delta t,u(t_{j}+c_{i}\Delta
t))\big{)}\Big{\|}_{1}\\\ &:={\rm II}_{1}+{\rm II}_{2}.\end{split}$ (3.48)

Using the similar arguments for deriving (3.35) and (3.2) and (3.3) in Lemma
3.1, we can obtain

$\begin{split}{\rm II}_{1}\lesssim\,&\textstyle\sum\limits_{i=1}^{s}(\Delta
t)^{\frac{1}{2}}\|(\Delta t)^{\frac{1}{2}}L_{h}^{\frac{1}{2}}b_{i}(-\Delta
tL_{h})\|_{0}\textstyle\max\limits_{1\leq i\leq s}\|P_{h}f(t_{n-1}+c_{i}\Delta
t,U_{n-1,i})\\\ &\qquad\qquad-P_{h}f(t_{n-1}+c_{i}\Delta
t,u(t_{n-1}+c_{i}\Delta t))\|_{0}\\\ \lesssim\,&(\Delta
t)^{\frac{1}{2}}\textstyle\max\limits_{1\leq i\leq
s}\|P_{h}f(t_{n-1}+c_{i}\Delta t,U_{n-1,i})\\\ &\qquad\qquad-
P_{h}f(t_{n-1}+c_{i}\Delta t,u(t_{n-1}+c_{i}\Delta t))\|_{0}\\\
\lesssim\,&(\Delta t)^{\frac{1}{2}}\textstyle\max\limits_{1\leq i\leq
s}\|E_{n-1,i}\|_{1}+h,\end{split}$ (3.49)

and

$\begin{split}{\rm II}_{2}\lesssim\,&\Big{\|}\Delta
tL_{h}^{\frac{1}{2}}\textstyle\sum\limits_{j=0}^{n-2}e^{-(n-1-j)\Delta
tL_{h}}\Big{\|}_{0}\textstyle\sup\limits_{0\leq t\leq
T}\|P_{h}f(t,u(t))-P_{h}f(t,u_{h}(t))\|_{0}\\\
&+\textstyle\sum\limits_{j=0}^{n-2}\Delta
t\|L_{h}^{\frac{1}{2}}e^{-(n-j-1)\Delta
tL_{h}}\|_{0}\textstyle\max\limits_{1\leq i\leq s}\|P_{h}f(t_{j}+c_{i}\Delta
t,U_{ji})\\\ &\qquad-P_{h}f(t_{j}+c_{i}\Delta t,u_{h}(t_{j}+c_{i}\Delta
t))\|_{0}\\\ \lesssim\,&\Delta
t\textstyle\sum\limits_{j=0}^{n-2}t_{n-j-1}^{-\frac{1}{2}}\textstyle\max\limits_{1\leq
i\leq s}\|P_{h}f(t_{j}+c_{i}\Delta t,U_{ji})\\\ &\qquad-
P_{h}f(t_{j}+c_{i}\Delta t,u_{h}(t_{j}+c_{i}\Delta t))\|_{0}+h\\\
\lesssim\,&\Delta
t\textstyle\sum\limits_{j=0}^{n-2}t_{n-j-1}^{-\frac{1}{2}}\textstyle\max\limits_{1\leq
i\leq s}\|E_{ji}\|_{1}+h.\end{split}$ (3.50)

With the help of the estimates (3.48), (3.49) and (3.50), it follows from the
relation (3.30) and (3.31b) in Lemma 3.4 (with $r=2$) that

$\begin{split}\|e_{n}\|_{1}\leq\,&\Big{\|}\Delta
t\textstyle\sum\limits_{j=0}^{n-1}e^{-(n-1-j)\Delta
tL_{h}}\textstyle\sum\limits_{i=1}^{s}b_{i}(-\Delta
tL_{h})\big{(}P_{h}f(t_{j}+c_{i}\Delta t,U_{ji})\\\
&\;-P_{h}f(t_{j}+c_{i}\Delta t,u(t_{j}+c_{i}\Delta
t))\big{)}\Big{\|}_{1}+\Big{\|}\textstyle\sum\limits_{j=0}^{n-1}e^{-j\Delta
tL_{h}}\delta_{n-j}^{[2]}\Big{\|}_{1}\\\ \lesssim\,&(\Delta
t)^{\frac{1}{2}}\textstyle\max\limits_{1\leq i\leq s}\|E_{n-1,i}\|_{1}+\Delta
t\textstyle\sum\limits_{j=0}^{n-2}t_{n-j-1}^{-\frac{1}{2}}\textstyle\max\limits_{1\leq
i\leq s}\|E_{ji}\|_{1}\\\ &+\Delta t^{2}\textstyle\sup\limits_{0\leq t\leq
T}\|f^{(2)}(t,u(t))\|_{1}+h\\\ \lesssim\,&\Delta
t\textstyle\sum\limits_{j=0}^{n-1}t_{n-j-1}^{-\frac{1}{2}}\textstyle\max\limits_{1\leq
i\leq s}\|E_{ji}\|_{1}+(\Delta t)^{2}+h.\end{split}$

This combined with the estimation of $\|E_{j,i}\|_{1}$ in Lemma 3.5 and the
discrete Gronwall inequality leads to

$\|u_{h}^{n}-u_{h}(t_{n})\|_{1}\lesssim(\Delta t)^{2}+h.$ (3.51)

Finally, the combination of (3.22), (3.24) and (3.51) immediately gives
(3.47). ∎

######  Remark 3.1.

It is worth noting that no restriction on the time step size $\Delta t$ is
imposed in Theorems 3.1 and 3.2, which implies that the proposed EIFE method
is stable with large time stepping.

######  Remark 3.2.

When $s\geq 3$, the order conditions for the explicit exponential Runge-Kutta
method do not hold automatically like the case of $s=1$ or $2$, and
convergence analysis of the proposed EIFE method (2.12) becomes much more
complicated since the estimates of $\|e_{n}\|_{1}$ and $\|E_{ni}\|_{1}$ will
be coupled together. We refer the reader to [18] for some details along this
line, and rigorous error estimates of higher-order EIFE schemes would be an
interesting open question.

##  4 FFT-based fast solution of the EIFE method

We now present fast solution of the EIFE method (2.12) for solving the
semilinear parabolic equation (1.1). The main idea is to take the advantage of
simultaneously diagonalizing the mass and coefficient matrices of the finite
element system (2.5) with an orthogonal matrix whose multiplication with a
vector can be efficiently implemented by FFT and tensor product spectral
decomposition.

We first focus on the case of homogeneous Dirichlet boundary condition (2.1),
and the extension to that of nonhomogeneous Dirichlet boundary condition is
straightforward by incorporating the given values of the exact solution on the
boundary nodes into the discrete method as part of $f$. Let us take the three
dimensional problem for illustration, which also naturally works in any other
$d$ dimensions. Suppose that $\Omega:=\prod_{i=1}^{3}[a_{i},b_{i}]$ is divided
into $N_{x},N_{y},N_{z}$ parts uniformly along the $x,y,z$ directions with the
meshsizes $h_{x},h_{y},h_{z}$, respectively. Thus, as given in Section 2, the
finite element space $V_{h}$ is spanned by
$\\{\phi_{i}(x)\phi_{j}(y)\phi_{k}(z)\\}$ where
$i=1,\cdots,N_{x}-1,j=1,\cdots,N_{y}-1,k=1,\cdots,N_{z}-1$. Therefore, the
finite element solution $u_{h}(t)$ of the semi-discrete (in space) problem
(2.5) can be expressed as

$u_{h}(t,x,y,z)=\textstyle\sum\limits_{i=1}^{N_{x}-1}\textstyle\sum\limits_{j=1}^{N_{y}-1}\textstyle\sum\limits_{k=1}^{N_{z}-1}u_{i,j,k}(t)\phi_{i}(x)\phi_{j}(y)\phi_{k}(z).$

Define $U(t)=[u_{i,i,k}(t)]_{(N_{x}-1)\times(N_{y}-1)\times(N_{z}-1)}$ and let
$F(t,U)$ be the tensor with entries
$(f(t,u_{i,j,k}(t)),\phi_{i}(x)\phi_{j}(y)\phi_{k}(z)).$ As in [26], we
introduce some tensor operations for later uses. Given three matrices
$M_{x}\in{\mathbb{R}}^{(N_{x}-1)\times(N_{x}-1)}$,
$M_{y}\in{\mathbb{R}}^{(N_{y}-1)\times(N_{y}-1)}$,
$M_{z}\in{\mathbb{R}}^{(N_{z}-1)\times(N_{z}-1)}$, let us define the following
operations with $U$:
$\big{(}M_{x}ⓧU\big{)}_{ijk}:=\textstyle\sum\limits_{r=1}^{N_{x}-1}(M_{x})_{ir}U_{rjk}$,
$\big{(}M_{y}ⓨU\big{)}_{ijk}:=\textstyle\sum\limits_{r=1}^{N_{y}-1}(M_{y})_{jr}U_{irk}$,
$\big{(}M_{z}ⓩU\big{)}_{ijk}:=\textstyle\sum\limits_{r=1}^{N_{z}-1}(M_{z})_{kr}U_{ijr}$.

With the help of the above operations, we have by some direct manipulations
that the variational scheme (2.5) can be reformulated to the following form

$\left\\{\begin{split}&A_{x}ⓧA_{y}ⓨA_{z}ⓩ\frac{\,{\rm d}{U}}{\,{\rm
d}t}+D(B_{x}ⓧA_{y}ⓨA_{z}ⓩU\\\
&\qquad\qquad\qquad+A_{x}ⓧB_{y}ⓨA_{z}ⓩU+A_{x}ⓧA_{y}ⓨB_{z}ⓩU)=F(t,U),\\\
&U(0)=\\{P_{h}u_{0}\\},\end{split}\right.$ (4.1)

where

$A_{\Diamond}=\textstyle\frac{h_{\Diamond}}{6}R_{N_{\Diamond}-1},\quad
B_{\Diamond}=\textstyle\frac{1}{h_{\Diamond}}G^{D}_{N_{\Diamond}-1},\quad\mbox{for}\;\;{\Diamond}=x,y,z,$

and

$\begin{split}R_{p}=\left[\begin{matrix}4&1&&\\\ 1&4&\ddots&\\\
&\ddots&\ddots&1\\\ &&1&4\end{matrix}\right]_{p\times p},\quad
G^{D}_{p}=\left[\begin{matrix}2&-1&&\\\ -1&2&\ddots&\\\ &\ddots&\ddots&-1\\\
&&-1&2\end{matrix}\right]_{p\times p}.\end{split}$

An important observation is that $A_{x}$ and $B_{x}$ can be diagonalized
simultaneously by an orthogonal matrix ([22]), and so do the other two pairs
$A_{y}$ and $B_{y}$, $A_{z}$ and $B_{z}$. Specifically, we have the following
decompositions:

$A_{\Diamond}=P_{\Diamond}\Lambda_{A_{\Diamond}}P_{\Diamond}^{T},\quad
B_{\Diamond}=P_{\Diamond}\Lambda_{B_{\Diamond}}P_{\Diamond}^{T},\quad\mbox{for}\;\;{\Diamond}=x,y,z$

with

$(P_{\Diamond})_{i,j}=\sin\Big{(}\textstyle\frac{ij\pi}{N_{\Diamond}}\Big{)},\quad
i,j=1,\cdots,N_{\Diamond}-1,$
$\Lambda_{A_{\Diamond}}=\text{diag}(\lambda^{A_{\Diamond}}_{1},\cdots,\lambda^{A_{\Diamond}}_{N_{\Diamond}-1}),\quad\lambda^{A_{\Diamond}}_{i}=\textstyle\frac{h_{\Diamond}}{6}\Big{(}6-4\sin\Big{(}\frac{i\pi}{2N_{\Diamond}}\Big{)}\Big{)},\quad
i=1,\cdots,N_{\Diamond}-1,$
$\Lambda_{B_{\Diamond}}=\text{diag}(\lambda^{B_{\Diamond}}_{1},\cdots,\lambda^{B_{\Diamond}}_{N_{\Diamond}-1}),\quad\lambda^{B_{\Diamond}}_{i}=\textstyle\frac{4}{h_{\Diamond}}\Big{(}\sin^{2}\Big{(}\frac{i\pi}{2N_{\Diamond}}\Big{)}\Big{)},\quad
i=1,\cdots,N_{\Diamond}-1.$

Let us define $\widetilde{U}=P_{x}^{T}ⓧP_{y}^{T}ⓨP_{z}^{T}ⓩU$ (i.e.,
$U=P_{x}ⓧP_{y}ⓨP_{z}ⓩ\widetilde{U}$). Then the equation (4.1) can be
transformed to

$\begin{split}&\Lambda_{A_{x}}ⓧ\Lambda_{A_{y}}ⓨ\Lambda_{A_{z}}ⓩ\frac{\,{\rm
d}\widetilde{U}}{\,{\rm
d}t}+D\Big{(}\Lambda_{B_{x}}ⓧ\Lambda_{A_{y}}ⓨ\Lambda_{A_{z}}ⓩ+\Lambda_{A_{x}}ⓧ\Lambda_{B_{y}}ⓨ\Lambda_{A_{z}}ⓩ\\\
&\qquad\quad+\Lambda_{A_{x}}ⓧ\Lambda_{A_{y}}ⓨ\Lambda_{B_{z}}ⓩ\Big{)}\widetilde{U}=P_{x}^{T}ⓧP_{y}^{T}ⓨP_{z}^{T}ⓩF(t,P_{x}ⓧP_{y}ⓨP_{z}ⓩ\widetilde{U}).\end{split}$
(4.2)

Similar to [45], we define the tensor
$H=(h_{ijk})_{(N_{x}-1)\times(N_{y}-1)\times(N_{z}-1)}$ with

$h_{ijk}=D\Big{(}\textstyle\frac{(\Lambda_{B_{x}})_{ii}}{(\Lambda_{A_{x}})_{ii}}+\frac{(\Lambda_{B_{y}})_{jj}}{(\Lambda_{A_{y}})_{jj}}+\frac{(\Lambda_{B_{z}})_{kk}}{(\Lambda_{A_{z}})_{kk}}\Big{)}$

and another tensor
$\widehat{H}=(\widehat{h}_{ijk})_{(N_{x}-1)\times(N_{y}-1)\times(N_{z}-1)}$
with

$\hat{h}_{ijk}=\textstyle\frac{1}{(\Lambda_{A_{x}})_{ii}(\Lambda_{A_{y}})_{jj}(\Lambda_{A_{z}})_{kk}}.$

Define the operator $(e^{*})$ as taking the exponential of each entry of a
tensor as $((e^{*})^{H})_{ijk}=e^{h_{ijk}},$ and another operator $\odot$ for
element by element multiplication between two arrays of same sizes as $(A\odot
B)_{ijk}=A_{ijk}B_{ijk}.$ Thus we can rewrite the system (4.2) to the
following ODE system:

$\frac{\,{\rm d}\widetilde{U}}{\,{\rm
d}t}+H\odot\widetilde{U}=\widehat{H}\odot\big{(}P_{x}^{T}ⓧP_{y}^{T}ⓨP_{z}^{T}ⓩF(t,P_{x}ⓧP_{y}ⓨP_{z}ⓩ\widetilde{U})\big{)}.$
(4.3)

Hence, by the Duhamel principle, we can get the following expression

$\begin{split}&\widetilde{U}(t_{n+1})=(e^{*})^{-H\Delta
t}\odot\widetilde{U}(t_{n})+\int_{t_{n}}^{t_{n+1}}(e^{*})^{-H(t_{n+1}-\tau)}\odot\widehat{H}\\\
&\qquad\qquad\qquad\qquad\odot\big{(}P_{x}^{T}ⓧP_{y}^{T}ⓨP_{z}^{T}ⓩF(\tau,P_{x}ⓧP_{y}ⓨP_{z}ⓩ\widetilde{U})\big{)}\,{\rm
d}\tau,\end{split}$ (4.4)

which is equivalently the tensor representation of (2.9).

We remark that all the tensor product operations can be realized by the FFT.
More precisely, the above transformation from $U$ to $\widetilde{U}$ (resp.
from $\widetilde{U}$ to $U$) can be implemented by applying a Discrete Sine
Transform (DST) (resp. an inverse DST, abbreviated as iDST) in each dimension
separately. That means, for a $d$-dimensional problem, the computational cost
of the EIFE method (2.12) is of $O(\log_{2}(N)\prod_{i=1}^{d}N_{i})$ per time
step, where $N=\textstyle\max\\{N_{1},\cdots,N_{d}\\}$.

Next, we also would like to briefly discuss the equation (1.1) with periodic
boundary condition (although not theoretically analyzed in this paper). In
this case, the finite element space $V_{h}$ is spanned by
$\\{\phi_{i}(x)\phi_{j}(y)\phi_{k}(z)\\}$ where
$i=1,\cdots,N_{x},j=1,\cdots,N_{y},k=1,\cdots,N_{z}$ by periodic extension.
The same EIFE method and efficient implementation algorithm can be directly
applied with slight modifications. In fact, the unknowns in this case is a
three-order tensor given by $U(t)=[u_{i,j,k}(t)]_{N_{x}\times N_{y}\times
N_{z}}$, and the corresponding mass and coefficient matrices become
respectively

$\begin{split}A_{\Diamond}=\textstyle\frac{h_{\Diamond}}{6}R_{N_{\Diamond}},\quad
B_{\Diamond}=\textstyle\frac{1}{h_{\Diamond}}G^{P}_{N_{\Diamond}},\quad\mbox{for}\;\;{\Diamond}=x,y,z,\end{split}$

where

$G^{P}_{p}=\left[\begin{matrix}2&-1&&-1\\\ -1&2&\ddots&\\\
&\ddots&\ddots&-1\\\ -1&&-1&2\end{matrix}\right]_{p\times p}.$

Then we have

$A_{\Diamond}=P_{\Diamond}\Lambda_{A_{\Diamond}}P_{\Diamond}^{T},\quad
B_{\Diamond}=P_{\Diamond}\Lambda_{B_{\Diamond}}P_{\Diamond}^{T},\quad{\Diamond}=x,y,z$

with

$(P_{\Diamond})_{i,j}=\sin\Big{(}\textstyle\frac{ij\pi}{N_{\Diamond}}\Big{)},\quad
i,j=1,\cdots,N_{\Diamond},$
$\Lambda_{A_{\Diamond}}=\text{diag}(\lambda^{A_{\Diamond}}_{1},\cdots,\lambda^{A_{\Diamond}}_{N_{\Diamond}}),\quad\lambda^{A_{\Diamond}}_{i}=\textstyle\frac{h_{\Diamond}}{6}\Big{(}6-4\sin\Big{(}\frac{(i-1)\pi}{2N_{\Diamond}}\Big{)}\Big{)},\quad
i=1,\cdots,N_{\Diamond},$
$\Lambda_{B_{\Diamond}}=\text{diag}(\lambda^{B_{\Diamond}}_{1},\cdots,\lambda^{B_{\Diamond}}_{N_{\Diamond}}),\quad\lambda^{B_{\Diamond}}_{i}=\textstyle\frac{4}{h_{\Diamond}}\Big{(}\sin^{2}\Big{(}\frac{(i-1)\pi}{2N_{\Diamond}}\Big{)}\Big{)},\quad
i=1,\cdots,N_{\Diamond}.$

The other steps of implementation are similar to the ones as in the Dirichlet
case, but note that the transformation from $U$ to $\widetilde{U}$ (resp. from
$\widetilde{U}$ to $U$) in the case can be implemented by applying a Discrete
Fourier Transform (DFT) (resp. an inverse DFT, abbreviated as iDFT) in each
dimension separately.

##  5 Numerical experiments

In this section we will present some numerical experiments to verify the error
estimates obtained in Section 3 and demonstrate the performance of the EIFE
method. All tests are done using Matlab on a Laptop with Intel i5-8250U,
1.80GHz CPU and 8GB memory. Specifically, we choose the EIFE1 scheme (2.13)
and the EIFE2 scheme (2.14) with $c_{2}=\frac{1}{2}$.

###  5.1 Convergence tests

We verify the error estimates obtained in Theorem 3.1 for EIFE1 scheme and
Theorem 3.2 for EIFE2 scheme, and numerical errors
$\|u(t_{n})-u_{h}^{n}\|_{0}$ and $\|u(t_{n})-u_{h}^{n}\|_{1}$ are all
evaluated at the terminal time $T$.

######  Example 5.1.

In this example, we consider the following two-dimensional linear reaction-
diffusion problem with homogeneous Dirichlet boundary condition:

$\left\\{\begin{split}&u_{t}=\frac{1}{2}\Delta
u-\frac{1}{2}\pi^{2}u+\frac{1}{2}\pi^{2}e^{-\pi^{2}t}\sin(\pi x)\sin(\pi
y),\quad(x,y)\in\Omega,\,0\leq t\leq T\\\ &u(0,x,y)=(\sin(\pi x)-1)\sin(\pi
y),\quad(x,y)\in\Omega,\\\ \end{split}\right.$

where $\Omega=(\frac{1}{2},\frac{5}{2})\times(0,1)$ and the terminal time
$T=1$. The exact solution is given by $u(t,x,y,z)=e^{-\pi^{2}t}(\sin(\pi
x)-1)\sin(\pi y)$.

For the spatial accuracy tests, we run the EIFE2 scheme with fixed
$N_{T}=1024$ (i.e., $\Delta t=T/N_{T}=1/1024$) and uniformly refined spatial
meshes with $N_{x}\times N_{y}=8\times 4$, $16\times 8$, $32\times 16$ and
$64\times 32$, respectively, so that the spatial mesh sizes are much coarser
compared to the time step size. For the temporal accuracy tests, we run the
EIFE1 and EIFE2 schemes with fixed $N_{x}\times N_{y}=2048\times 1024$ and
uniform time step with $N_{T}=16,32,64,128$. Also, the overall cost per time
step of EIFE method is tested with fixed $N_{T}=50$ and uniform spatial meshes
with $N_{x}\times N_{y}=512\times 256$, $1024\times 512$, $2048\times 1024$
and $4096\times 2048$. Since the analysis of computing cost of EIFE1 scheme is
same as EIFE2 scheme, we only test the running time with EIFE2 scheme. All
numerical results are reported in Table 2, including the solution errors
measured in the $L^{2}$ and $H^{1}$ norms and corresponding convergence rates.
We observe the roughly second-order spatial convergence with respect to both
the $L^{2}$ and $H^{1}$ norms, so the convergence order is one-order higher
than the expected value one in terms of the $H^{1}$-norm in this case. It is
also easy to find the first-order temporal convergence for the EIFE1 scheme
and the second-order temporal convergence for the EIFE2 scheme in both the
$L^{2}$ and $H^{1}$ norms, which coincide very well with the error estimates
derived in Theorems 3.1 and 3.2. Table 3 reports the average CPU time costs
(seconds) per step for the EIFE2 scheme and corresponding growth factors along
the refinement of the spatial mesh. The results clearly show that the
computational cost grows almost linearly along with the number of mesh nodes,
which matches well with the property of FFT and demonstrates the high
efficiency of our EIFE method.

Table 2: Numerical results on the solution errors measured in the $L^{2}$ and $H^{1}$ norms and corresponding convergence rates for the EIFE1 and EIFE2 schemes in Example 5.1.  $N_{T}$ | $N_{x}\times N_{y}$ | $\|u_{h}^{n}-u(t_{n})\|_{0}$ | CR | $\|u_{h}^{n}-u(t_{n})\|_{1}$ | CR
---|---|---|---|---|---
Spatial accuracy tests for EIFE2
1024 | $8\times 4$ | 2.1975e-05 | - | 5.8018e-05 | -
1024 | $16\times 8$ | 6.8220e-06 | 1.69 | 2.0817e-05 | 1.48
1024 | $32\times 16$ | 1.8046e-06 | 1.92 | 6.2344e-06 | 1.74
1024 | $64\times 32$ | 4.5693e-07 | 1.98 | 1.6276e-06 | 1.94
Temporal accuracy tests for EIFE1
16 | $2048\times 1024$ | 1.6807e-05 | - | 6.1597e-05 | -
32 | $2048\times 1024$ | 9.2840e-06 | 0.86 | 3.3985e-05 | 0.86
64 | $2048\times 1024$ | 4.8687e-06 | 0.93 | 1.7814e-05 | 0.93
128 | $2048\times 1024$ | 2.4919e-06 | 0.97 | 9.1158e-06 | 0.97
Temporal accuracy tests for EIFE2
16 | $2048\times 1024$ | 8.5334e-06 | - | 3.1234e-05 | -
32 | $2048\times 1024$ | 1.6087e-06 | 2.41 | 5.8873e-06 | 2.41
64 | $2048\times 1024$ | 3.5955e-07 | 2.16 | 1.3157e-06 | 2.16
128 | $2048\times 1024$ | 8.4991e-08 | 2.08 | 3.1121e-07 | 2.08
Table 3: The average CPU time costs (seconds) per step under different spatial meshes and corresponding growth factors with respect to the number of mesh nodes for the EIFE2 scheme in Example 5.1. $N_{x}\times N_{y}$ | Average CPU time | Growth
---|---|---
| cost per step | factor
$512\times 256$ | 1.121 | -
$1024\times 512$ | 4.444 | 0.99
$2048\times 1024$ | 17.516 | 0.99
$4096\times 2048$ | 75.944 | 1.06

######  Example 5.2.

In this example, we consider the traveling wave problem governed by the
following 3D Allen-Cahn equation with double-well potential function:

$\left\\{\begin{split}&u_{t}=\Delta
u-\frac{1}{\epsilon^{2}}(u^{3}-u),\quad(x,y,z)\in\Omega,\;0\leq t\leq T,\\\
&u(0,x,y,z)=\frac{1}{2}(1-\tanh(\frac{x}{2\sqrt{2}\epsilon})),\quad(x,y,z)\in\Omega,\\\
\end{split}\right.$

where $\Omega=(0,\sqrt{2})\times(0,\frac{1}{8})\times(0,\frac{1}{8})$. The
exact solution is given by
$u(t,x,y,z)=\frac{1}{2}(1-\tanh(\frac{x-st}{2\sqrt{2}\epsilon}))$, where
$s=\frac{3}{\sqrt{2}\epsilon}$, and the Dirichlet boundary condition is
correspondingly imposed, which is clearly nonhomogeneous. The terminal time is
taken to be $T=\frac{3\sqrt{2}\epsilon}{5}$.

We set the interface thickness parameter $\epsilon=0.05$. For the spatial
accuracy tests, we run the EIFE2 scheme with fixed $N_{T}=1024$ and uniformly
refined spatial meshes with $N_{x}\times N_{y}=64\times 4\times 4$, $128\times
8\times 8$, $256\times 16\times 16$ and $512\times 32\times 32$, respectively,
so that the spatial mesh sizes are much coarser compared to the time step
size. For the temporal accuracy tests, we run the EIFE1 and EIFE2 schemes with
fixed $N_{x}\times N_{y}=1024\times 64\times 64$ and uniformly refined time
step sizes with $N_{T}=16,32,64$ and $128$. All numerical results are reported
in Table 4, including the errors measured in the $L^{2}$ and $H^{1}$ norms and
corresponding convergence rates. We still observe the roughly second-order
spatial convergence in the $L^{2}$-norm but it drops to between 1.11 to 1.28
(but still higher than one) in the $H^{1}$-norm. It is also seen that the
temporal convergence is just slightly higher than one for the EIFE1 scheme and
two for the EIFE2 scheme, which basically match the theoretical results.

Table 4: Numerical results on the solution errors measured in the $L^{2}$ and $H^{1}$ norms and corresponding convergence rates for the EIFE1 and EIFE2 schemes in Example 5.2.  $N_{T}$ | $N_{x}\times N_{y}\times N_{z}$ | $\|u_{h}^{n}-u(t_{n})\|_{0}$ | CR | $\|u_{h}^{n}-u(t_{n})\|_{1}$ | CR
---|---|---|---|---|---
Spatial accuracy tests for EIFE2
1024 | $64\times 4\times 4$ | 5.6535e-06 | - | 1.0563e-04 | -
1024 | $128\times 8\times 8$ | 1.5135e-06 | 1.90 | 4.8833e-05 | 1.11
1024 | $256\times 16\times 16$ | 4.3188e-07 | 1.81 | 2.0047e-05 | 1.28
1024 | $512\times 32\times 32$ | 1.2193e-07 | 1.82 | 8.2719e-06 | 1.28
Temporal accuracy tests for EIFE1
16 | $1024\times 64\times 64$ | 2.4559e-04 | - | 3.3100e-02 | -
32 | $1024\times 64\times 64$ | 7.5000e-05 | 1.71 | 1.0500e-02 | 1.66
64 | $1024\times 64\times 64$ | 2.9850e-05 | 1.33 | 4.2000e-03 | 1.32
128 | $1024\times 64\times 64$ | 1.2765e-05 | 1.23 | 1.9000e-03 | 1.14
Temporal accuracy tests for EIFE2
16 | $1024\times 64\times 64$ | 9.5466e-04 | - | 3.6900e-02 | -
32 | $1024\times 64\times 64$ | 3.0789e-05 | 4.95 | 2.5000e-03 | 3.88
64 | $1024\times 64\times 64$ | 4.6133e-06 | 2.74 | 4.8362e-04 | 2.37
128 | $1024\times 64\times 64$ | 9.5582e-07 | 2.27 | 1.1057e-04 | 2.13

###  5.2 3D Grain coarsening simulations

We now illustrate the performance of the proposed EIFE method through
numerical simulation of the 3D grain coarsening process. In particular, the
EIFE2 scheme is used for its higher-order accuracy in time.

######  Example 5.3.

We consider the grain coarsening process governed by the following 3D Allen-
Cahn equation with Flory-Huggins potential function:

$u_{t}=\epsilon^{2}\Delta
u+\frac{\theta}{2}\ln\frac{1-u}{1+u}+\theta_{c}u,\quad(x,y,z)\in\Omega,\;0\leq
t\leq T,$

where $\Omega=(0,1)^{3}$. The initial data is generated by random numbers on
each mesh point ranging from $-0.9$ to $0.9$, and the periodic boundary
condition is imposed. This equation can be regarded as the $L^{2}$ gradient
flow of the following energy functional

$E(u)=\int_{\Omega}\frac{\theta}{2}\Big{(}(1+u)\ln(1+u)+(1-u)\ln(1-u)\Big{)}-\frac{\theta_{c}}{2}u^{2}+\frac{\epsilon^{2}}{2}|\nabla
u|^{2}\,{\rm d}\bm{x},$

and thus the energy monotonically decays along the time.

We set the interface thickness parameter $\epsilon=0.01$ and the Flory-Huggins
potential parameters $\theta_{c}=1.6,\theta=0.8$. This problem satisfies the
maximum bound principle with the maximum bound value $\gamma\approx 0.9575$,
i.e., $|u(t,\bm{x})|\leq\gamma$ for all $\bm{x}\in\Omega$ and $t\geq 0$ [29,
6]. The terminal time is set to be $T=20$ for our simulation. We take the
spatial mesh with $N_{x}=N_{y}=N_{z}=128$ which means the mesh size $h=1/128$,
and the total time step size $N_{T}=2048$ (i.e., $\Delta t=T/N_{T}=5/512$).
Evolutions of the supremum norm and the energy of the numerical solutions
produced by the EIFE2 scheme are plotted in Fig. 1. We observe that the
maximum bound principle is well preserved and the energy also decays
monotonically along the time. Fig. 2 presents simulated phase structures of
the numerical solutions at times $t=$ 0, 2.5, 5, 10, 15, and 20, and from
which we can clearly observe the coarsening process of the 3D grains.

Figure 1: The evolutions of supremum norm (left) and energy (right) of the
numerical solution produced by the EIFE2 scheme for Example 5.3.

Figure 2: Simulated 3D phase structures at $t=0,2.5,5,10,15,20$ (from left to
right and top to bottom) produced by the EIFE2 scheme for Example 5.3.

##  6 Conclusions

In this paper, we develop an efficient EIFE method for solving a class of
semilinear parabolic equations taking the form (1.1) in regular domains, in
which the fully-discrete solution is obtained by using first the finite
element method for spatial discretization and then explicit exponential Runge-
Kutta approximation for temporal integration. The EIFE method allows for fast
implementation based on FFT and tensor product spectral decomposition. We
successfully derive optimal error estimates in the $H^{1}$-norm for the EIFE
method with one and two RK stages when the Dirichlet boundary condition s
imposed. Some numerical examples are also presented to demonstrate the
accuracy and high efficiency of the proposed method. Rigorous error analysis
of the EIFE method in the $L^{2}$-norm and for the model problem with periodic
boundary condition still remain to be explored. In addition, the numerical
method and corresponding error analysis framework developed in this paper also
naturally enable us to further investigate the localized ETD methods [30, 17]
with solid theoretical support.

## References

  * [1] S. Allen and J. Cahn.  A microscopic theory for antiphase domain boundary motion and its application to antiphase domain coarsening.  Acta Metall., 27:1085–1095, 06 1979.
  * [2] Uri M. Ascher, Steven J. Ruuth, and Brian T. R. Wetton.  Implicit-explicit methods for time-dependent partial differential equations.  SIAM J. Numer. Anal., 32(3):797–823, 1995.
  * [3] John P. Boyd.  Chebyshev and Fourier Spectral Methods.  Dover Publications, Inc., Mineola, NY, second edition, 2001.
  * [4] James H. Bramble and J. Xu.  Some estimates for a weighted $L^{2}$ projection.  Math. Comp., 56(194):463–476, 1991.
  * [5] Susanne C. Brenner and L. Ridgway Scott.  The Mathematical Theory of Finite Element Methods.  Springer, New York, third edition, 2008.
  * [6] W. Chen, J. Jing, C. Wang, X. Wang, and Steven M. Wise.  A modified Crank-Nicolson numerical scheme for the Flory-Huggins Cahn-Hilliard model.  Commun. Comput. Phys., 31(1):60–93, 2022.
  * [7] Philippe G. Ciarlet.  The Finite Element Method for Elliptic Problems.  North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978.
  * [8] S. M. Cox and P. C. Matthews.  Exponential time differencing for stiff systems.  J. Comput. Phys., 176(2):430–455, 2002.
  * [9] J. Dixon and S. McKee.  Weakly singular discrete Gronwall inequalities.  Z. Angew. Math. Mech., 66(11):535–544, 1986.
  * [10] Tobin A. Driscoll.  A composite Runge-Kutta method for the spectral solution of semilinear PDEs.  J. Comput. Phys., 182(2):357–367, 2002.
  * [11] Q. Du and X. Feng.  The phase field method for geometric moving interfaces and their numerical approximations.  In Geometric Partial Differential Equations. Part I, volume 21, pages 425–508. Elsevier/North-Holland, Amsterdam, 2020.
  * [12] Q. Du, M. D. Gunzburger, and J. S. Peterson.  Analysis and approximation of the Ginzburg–Landau model of superconductivity.  SIAM Rev., 34(1):54–81, 1992.
  * [13] Q. Du, L. Ju, X. Li, and Z. Qiao.  Maximum principle preserving exponential time differencing schemes for the nonlocal Allen–Cahn equation.  SIAM J. Numer. Anal., 57(2):875–898, 2019.
  * [14] Q. Du, L. Ju, X. Li, and Z. Qiao.  Maximum bound principles for a class of semilinear parabolic equations and exponential time differencing schemes.  SIAM Rev., 63(2):317–359, 2021.
  * [15] Q. Du and W. Zhu.  Stability analysis and application of the exponential time differencing schemes.  J. Comput. Math., 22(2):200–209, 2004.
  * [16] Xiaobin Feng and Andreas Prohl.  Numerical analysis of the allen-cahn equation and approximation for mean curvature flows.  Numerische Mathematik, 94:33–65, 2003.
  * [17] Thi-Thao-Phuong Hoang, L. Ju, and Z. Wang.  Overlapping localized exponential time differencing methods for diffusion problems.  Commun. Math. Sci., 16(6):1531–1555, 2018.
  * [18] M. Hochbruck and A. Ostermann.  Explicit exponential Runge-Kutta methods for semilinear parabolic problems.  SIAM J. Numer. Anal., 43(3):1069–1090, 2005.
  * [19] M. Hochbruck and A. Ostermann.  Exponential Runge-Kutta methods for parabolic problems.  Appl. Numer. Math., 53(2-4):323–339, 2005.
  * [20] M. Hochbruck and A. Ostermann.  Exponential integrators.  Acta Numer., 19:209–286, 2010.
  * [21] M. Hochbruck, A. Ostermann, and J. Schweitzer.  Exponential Rosenbrock-type methods.  SIAM J. Numer. Anal., 47(1):786–803, 2008/09.
  * [22] J. Huang, L. Ju, and B. Wu.  A fast compact exponential time differencing method for semilinear parabolic equations with Neumann boundary conditions.  Appl. Math. Lett., 94:257–265, 2019.
  * [23] J. Huang, L. Ju, and B. Wu.  A fast compact time integrator method for a family of general order semilinear evolution equations.  J. Comput. Phys., 393:313–336, 2019.
  * [24] L. Isherwood, Z. J. Grant, and S. Gottlieb.  Strong stability preserving integrating factor Runge-Kutta methods.  SIAM J. Numer. Anal., 56(6):3276–3307, 2018.
  * [25] L. Ju, X. Li, Z. Qiao, and J. Yang.  Maximum bound principle preserving integrating factor Runge-Kutta methods for semilinear parabolic equations.  J. Comput. Phys., 439:110405, 2021.
  * [26] L. Ju, J. Zhang, L. Zhu, and Q. Du.  Fast explicit integration factor methods for semilinear parabolic equations.  J. Sci. Comput., 62(2):431–455, 2015.
  * [27] Aly-Khan Kassam and Lloyd N. Trefethen.  Fourth-order time-stepping for stiff PDEs.  SIAM J. Sci. Comput., 26(4):1214–1233, 2005.
  * [28] J. Douglas Lawson.  Generalized Runge-Kutta processes for stable systems with large Lipschitz constants.  SIAM J. Numer. Anal., 4:372–380, 1967.
  * [29] J. Li, X. Li, L. Ju, and X. Feng.  Stabilized integrating factor Runge-Kutta method and unconditional preservation of maximum bound principle.  SIAM J. Sci. Comput., 43(3):A1780–A1802, 2021.
  * [30] Xiao Li, Lili Ju, and Thi-Thao-Phuong Hoang.  Overlapping domain decomposition based exponential time differencing methods for semilinear parabolic equations.  BIT, 61(1):1–36, 2021.
  * [31] Vu Thai Luan and Alexander Ostermann.  Explicit exponential runge–kutta methods of high order for parabolic problems.  Journal of Computational and Applied Mathematics, 256:168–179, 2014\.
  * [32] Borislav V. Minchev and Will M. Wright.  A Review of Exponential Integrators for First Order Semilinear Problems.  preprint. Norwegian University of Science and Technology Trondheim, Norway, 2005.
  * [33] A. Mohebbi and M. Dehghan.  High-order solution of one-dimensional sine-Gordon equation using compact finite difference and DIRKN methods.  Math. Comput. Modelling, 51(5-6):537–549, 2010.
  * [34] Q. Nie, Y. Zhang, and R. Zhao.  Efficient semi-implicit schemes for stiff systems.  J. Comput. Phys., 214:521–537, 2006.
  * [35] David A. Pope.  An exponential method of numerical integration of ordinary differential equations.  Comm. ACM, 6:491–493, 1963.
  * [36] Frédéric Rousset and Katharina Schratz.  A general framework of low regularity integrators.  SIAM Journal on Numerical Analysis, 59(3):1735–1768, 2021.
  * [37] J. M. Sanz-Serna and M. P. Calvo.  Numerical Hamiltonian Problems, volume 7.  Chapman & Hall, London, 1994.
  * [38] Jie Shen, Jie Xu, and Jiang Yang.  A new class of efficient and robust energy stable schemes for gradient flows.  SIAM Review, 61(3):474–506, 2019.
  * [39] Roger Temam, editor.  The Mathematical Basis of Finite Element Methods, volume 343.  AMS Chelsea Publishing, 1984.
  * [40] Vidar Thomée.  Galerkin Finite Element Methods for Parabolic Problems.  Springer-Verlag, Berlin, second edition, 2006.
  * [41] D. Wang, L. Zhang, and Q. Nie.  Array-representation integration factor method for high dimensional systems.  J. Comput. Phys., 258:585–600, 2014.
  * [42] P. Whalen, M. Brio, and J. V. Moloney.  Exponential time-differencing with embedded Runge-Kutta adaptive step control.  J. Comput. Phys., 280:579–601, 2015.
  * [43] S. M. Wise, C. Wang, and J. S. Lowengrub.  An energy-stable and convergent finite-difference scheme for the phase field crystal equation.  SIAM Journal on Numerical Analysis, 47(3):2269–2288, 2009.
  * [44] Xiaofeng Yang.  Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends.  Journal of Computational Physics, 327:294–316, 2016.
  * [45] L. Zhu, L. Ju, and W. Zhao.  Fast high-order compact exponential time differencing Runge-Kutta methods for second-order semilinear parabolic equations.  J. Sci. Comput., 67(3):1043–1065, 2016.