Maximum bound principles for a class of semilinear parabolic equations and exponential time differencing schemes

Qiang Du²²2Department of Applied Physics and Applied Mathematics and Data Science Institute, Columbia University, New York, NY 10027, USA ([email protected]). Q. Du’s work is partially supported by US National Science Foundation grant DMS-1719699, US ARO MURI grant W911NF-15-1-0562 and US AFOSR MURI Center for Material Failure Prediction Through Peridynamics. Lili Ju³³3Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA ([email protected]). L. Ju’s work is partially supported by US National Science Foundation grant DMS-1818438 and US Department of Energy grants DE-SC0016540 and DE-SC0020270. Xiao Li⁴⁴4Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA. Current address: Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong ([email protected]). X. Li’s work is partially supported by National Natural Science Foundation of China grant 11801024. Zhonghua Qiao⁵⁵5Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong ([email protected]). Z. Qiao’s work is partially supported by the Hong Kong Research Council GRF grants 15302214 and 15325816 and the Hong Kong Polytechnic University fund 1-ZE33.

Abstract

The ubiquity of semilinear parabolic equations has been illustrated in their numerous applications ranging from physics, biology, to materials and social sciences. In this paper, we consider a practically desirable property for a class of semilinear parabolic equations of the abstract form $u_{t}={\mathcal{L}}u+f[u]$ with ${\mathcal{L}}$ being a linear dissipative operator and $f$ being a nonlinear operator in space, namely a time-invariant maximum bound principle, in the sense that the time-dependent solution $u$ preserves for all time a uniform pointwise bound in absolute value imposed by its initial and boundary conditions. We first study an analytical framework for some sufficient conditions on ${\mathcal{L}}$ and $f$ that lead to such a maximum bound principle for the time-continuous dynamic system of infinite or finite dimensions. Then, we utilize a suitable exponential time differencing approach with a properly chosen generator of contraction semigroup to develop first- and second-order accurate temporal discretization schemes, that satisfy the maximum bound principle unconditionally in the time-discrete setting. Error estimates of the proposed schemes are derived along with their energy stability. Extensions to vector- and matrix-valued systems are also discussed. We demonstrate that the abstract framework and analysis techniques developed here offer an effective and unified approach to study the maximum bound principle of the abstract evolution equation that cover a wide variety of well-known models and their numerical discretization schemes. Some numerical experiments are also carried out to verify the theoretical results.

keywords:

Semilinear parabolic equation, maximum bound principle, numerical approximation, exponential time differencing, energy stability, error estimate.

AMS:

35B50, 35K55, 65M12, 65R20

1 Introduction

Semilinear parabolic equations of the form

(1)

u_{t}={\mathcal{L}}u+f[u],

with $u$ being a time-dependent quantity of interests defined over a spatial domain $\Omega$ , ${\mathcal{L}}$ being a linear classic elliptic operator or its nonlocal variant, and $f$ representing a nonlinear operator, have been used to model numerous phenomena in nature. Many model equations like (1) and their solutions often satisfy some properties such as maximum principle, comparison principle, existence of invariant regions, and energy decay. These properties represent important physical features and are also essential for mathematical analysis and numerical simulations. For instance, classic reaction-diffusion equations can be seen as special cases of (1) where ${\mathcal{L}}$ is given by a diffusion operator and $f$ is a reaction source: the diffusion process causes the concentration $u$ of some substance to spread in space and the reaction drives the dynamics based on the concentration values. One illustration is the dimensionless ignition model [62] for a supercritical high activation energy thermal explosion of a solid fuel in a bounded container, described by

(2)

u_{t}=\Delta u+\mathrm{e}^{u}.

Since the reaction term $\mathrm{e}^{u}$ is always positive, the solution must reach its minimum either at the initial time or on the boundary of the domain, which is a popular form of the well-studied maximum principle for parabolic equations. Another popular example of (1) is the so-called Allen–Cahn equation [1], which takes on the form

(3)

u_{t}=\varepsilon^{2}\Delta u+u-u^{3},

with $\varepsilon>0$ reflecting the width of the transition regions. The equation is also a special case of the Ginzburg–Landau theory that has been introduced earlier in the modeling of superconductivity [41] and other phase transition problems. It is well-known that the equation (3) and its more general form (as discussed later) hold a maximum bound principle (MBP) [19, 29]: if the initial data and/or the boundary values are pointwisely bounded by $1$ in absolute value, then the absolute value of the solution is also bounded by $1$ everywhere and for all time. This MBP is related to the conventional version of maximum principle described for (2), but also differs somewhat in scope, so we use a new name to distinguish them. For the scalar model like the equation (3), it is equivalent to the existence of special upper and lower solutions (the constant functions $\overline{u}\equiv 1$ and $\underline{u}\equiv-1$ respectively) of (3) under Dirichlet or homogeneous Neumann boundary condition [81]. Moreover, it can be also described as the equation (3) having a time-invariant region [2, 63, 91] since the region $[-1,1]$ , of the value of the solution, remains unchanged during the time evolution. In this paper, we consider an abstract form of evolution equations given by (1) and focus on the problems having the MBP or the existence of the special form of the invariant region, namely, an upper bound on the pointwise absolute value of the solution (suitable defined for vector- and matrix-valued quantities). We note that the MBP is weaker than the conventional maximum principle in the sense that a problem satisfying a maximum principle must satisfy an MBP. Thus, the study of the conventional maximum principle, particularly with respect to the linear operator ${\mathcal{L}}$ , plays relevant and important roles.

Whether the equation (1) has a maximum principle or not depends highly on the property of the linear operator ${\mathcal{L}}$ . Concerning the linear operator ${\mathcal{L}}$ , one can find from the standard textbooks on partial differential equations (see, e.g., [28]) that the equation (1) with the uniformly elliptic linear operator ${\mathcal{L}}$ satisfies a maximum principle. During the past several decades, there have also been many studies devoted to the maximum principles for numerical approximations of linear elliptic operators. While it is too numerous to detail them all here, a partial list of earlier works includes the cases for finite difference method [5, 9, 83, 97], finite element method [3, 6, 7, 10, 61], collocation method [102, 103], and finite volume method [104]. In [9], a systematic analysis was presented on some sufficient conditions for the finite difference operator to satisfy the weak and strong maximum principles in the discrete sense. There were also works devoted to general analysis on algebraic properties of the approximating operators, which consist of the discrete maximum principles as special cases [32, 33, 73]. Extensions to parabolic cases have also been made. For example, [59] provided some sufficient conditions for the discrete maximum principles of the forward and backward Euler time-stepping methods, and additional researches can be found, e.g., in [30, 31, 60, 67, 69, 101]. In addition to the maximum principles for both linear elliptic differential operators and their finite dimensional discretizations, there are also some recent studies on the nonlocal analogues for linear nonlocal integral operators [22, 78, 94].

With the linear operator ${\mathcal{L}}$ satisfying the conditions to yield the maximum principle, a suitable nonlinear term may lead to the existence of time-invariant regions of (1). The MBP, that is a special invariant region of the Allen–Cahn equation (3), was proved in [29]. An important and natural question for numerical analysis is whether such an MBP could be preserved by some time-stepping schemes for discretizing (3). This question has been studied in a variety of works recently. In [92], the discrete MBPs of a finite difference semi-discrete scheme and its fully discrete approximations with forward and backward Euler time-stepping methods were obtained for (3) in one-dimensional space. Later, the first-order stabilized implicit-explicit schemes with finite difference spatial discretization were proved to preserve the MBP [93], which was then generalized by [88] to the case with more general nonlinear terms. The discrete MBP was also obtained in [100] by proving the uniform boundedness of the average $L^{p}$ norm of the solution for any even integer $p$ , followed by passing the limit as $p$ goes to infinity. A first-order exponential time differencing (ETD) scheme (or say, exponential integrator scheme) in the space-continuous setting was analyzed in [24], where some properties of the heat kernel were used. In addition, the MBP-preserving numerical schemes have been also studied for the analogues of (3), such as the fractional Allen–Cahn equation using the Crank–Nicolson time-stepping [51], the nonlocal Allen–Cahn equation by using first- and second-order ETD schemes [20], and the complex-valued Ginzburg–Landau model of superconductivity [19] by considering the finite volume method [13] and finite element method with the mass-lumping technique [14] in space with backward Euler time-stepping.

In this paper, we present a unified framework on the MBPs for a class of semilinear parabolic equations (1) with general boundary conditions and their spatially discretized systems, as well as the time-discrete analogues of numerical approximations based on the exponential time differencing method. The ETD method [4, 11, 48] comes from the variation-of-constants formula with the nonlinear terms approximated by polynomial interpolations, followed by exact integration of the resulting integrals. We mainly address the following question: Under what conditions does the equation (1) have the MBP and do numerical approximations of (1) preserve the MBP? We provide an abstract framework on mathematical and numerical analysis for the MBP of (1), where the ETD method is used for the temporal discretization. Our main contribution includes several aspects. First, we systematically formulate an abstract mathematical framework to illustrate the essential characteristics of the linear and nonlinear operators in (1) so that the model equation satisfies the MBP and the corresponding first- and second-order ETD schemes preserve the discrete MBP unconditionally. Second, we are able to present results valid for a large class of semilinear parabolic problems subject to different boundary conditions and constraints which significantly generalize those obtained in [20] (which is a specialized study focused only on the nonlocal Allen–Cahn equation with periodic boundary condition). We demonstrate that the theory also works for problems subject to either Dirichlet boundary condition or homogeneous Neumann boundary condition in the classic (local continuum) and nonlocal sense. These generalizations offer our theory a much wider range of applicability. Third, in order to establish a broad and abstract framework presented here, new analysis techniques are developed. They differ from the ones used in earlier works. Indeed, the abstraction allows us to discover the essential ingredients of the MBPs, which were absent from earlier analysis. In addition, we also derive error estimates for the approximate solutions of MBP-preserving ETD schemes, as well as their energy stability when applied to gradient flow models.

The necessity to take nonhomogeneous Dirichlet boundary condition into account comes from various practical considerations. For instance, phase separations occurring in hydrocarbon systems could be described by a diffuse-interface model with the Peng–Robinson equation of state [84], where a nonhomogeneous Dirichlet boundary condition is usually needed to keep the microstructures with certain phases on the boundary. In addition, nonhomogeneous Dirichlet boundary conditions are also necessary for many scalable and multiscale algorithms based on domain and subspace decompositions.

One of the distinctive features of the ETD schemes is the exact evaluation of the contribution of the linear operator, which provides good stability and accuracy even though the linear part has strong stiffness. Besides, the ETD schemes usually perform as efficient as an explicit scheme since the operator exponentials could be often implemented by some fast algorithms in regular domains. Such advantages lead to successful applications of ETD schemes on a large class of phase-field models which usually yield highly stiff ODE systems under spatial discretizations. We refer the readers to the literature, e.g., [20, 55, 57, 58, 105, 107].

The rest of this paper is organized as follows. In Section 2, we recall the semilinear equation (1) in a Banach space consisting of real scalar-valued continuous functions, declare the basic assumptions on the linear and nonlinear operators, and show that the model equation has a unique solution and satisfies the MBP under these assumptions. We also present a variety of concrete examples of the linear and nonlinear operators in the space-continuous and space-discrete settings covered by the theoretical framework. In Section 3, first- and second-order ETD time-discrete schemes are constructed and proved to preserve the discrete MBPs along with their convergence analysis. In addition, energy stability is also analyzed when the framework is applied to gradient flow models. In Section 4, we discuss the extensions to vector-valued and matrix-valued problems in the space-continuous setting with the complex-valued one being a special case. In Section 5, practical implementations of the ETD schemes are discussed and some numerical experiments are also performed to verify the theoretical results. Finally, concluding remarks are given in Section 6.

2 Maximum bound principle of the model equation

2.1 Abstract framework

Let us assume $\Omega$ is either a connected spatial region or a collection of isolated points in $\mathbb{R}^{d}$ . More precisely, we consider the following two situations:

(D1)

$\Omega$ is an open, connected and bounded set with a Lipschitz boundary denoted by $\partial\Omega$ , and $\Omega_{c}$ is a closed connected set disjoint with $\Omega$ but $\partial\Omega\subset\Omega_{c}$ ; denote $\overline{\Omega}=\Omega\cup\partial\Omega$ and $\widehat{\Omega}=\Omega\cup\Omega_{c}$ ;
(D2)

for a given pair of sets $\widetilde{\Omega}$ and $\widetilde{\Omega}_{c}$ that are described in (D1), and with $\Sigma$ consisting of all the nodes in a mesh partitioning $\widehat{\widetilde{\Omega}}$ , we let $\Omega=\widetilde{\Omega}\cap\Sigma$ , $\partial\Omega=\partial\widetilde{\Omega}\cap\Sigma$ , $\Omega_{c}=\widetilde{\Omega}_{c}\cap\Sigma$ , $\overline{\Omega}=\Omega\cup\partial\Omega$ , and $\widehat{\Omega}=\Omega\cup\Omega_{c}$ .

Note that in Case (D1), we have $\Omega_{c}=\partial\Omega$ for classic differential operators and $\Omega_{c}$ is usually a nonempty volume for nonlocal integral operators. When $\Omega_{c}=\partial\Omega$ , we define $\Omega^{*}=\Omega$ and $\Omega_{c}^{*}=\partial\Omega$ ; otherwise, we define $\Omega^{*}=\overline{\Omega}$ and $\Omega_{c}^{*}=\Omega_{c}\setminus\partial\Omega$ . For a problem with periodic boundary condition, we assume that its period cell is a rectangle in $\mathbb{R}^{d}$ , i.e., $\Omega=\prod_{i=1}^{d}(a_{i},b_{i})$ , and we use a special notation $\overline{\Omega}_{+}=\prod_{i=1}^{d}(a_{i},b_{i}]$ . Case (D2) corresponds to a discrete version of (D1). This allows us to unify the discussions presented later with the same set of notations for both space-continuous and space-discrete cases. For any set $D\subset\mathbb{R}^{d}$ , we denote by $C(D)$ the space consisting of real scalar-valued continuous functions defined on $D$ and by $C_{b}(D)$ the set of all bounded functions in $C(D)$ . Here, the continuity of functions is defined as follows [86]:

(4)

\text{$w$ is continuous at $\bm{x}^{*}\in D$}\iff\text{$\forall\,\bm{x}_{k}\to\bm{x}^{*}$ in $D$ implies $w(\bm{x}_{k})\to w(\bm{x}^{*})$}.

If $D$ is bounded and closed, then obviously $C(D)=C_{b}(D)$ .

Let $\mathcal{X}=C(\overline{\Omega})$ if $\Omega_{c}=\partial\Omega$ ; otherwise $\mathcal{X}=\{w:\widehat{\Omega}\to\mathbb{R}\,|\,w|_{\overline{\Omega}}\in C(\overline{\Omega})\text{ and }w|_{\Omega_{c}^{*}}\in C_{b}(\Omega_{c}^{*})\}$ and we write $\mathcal{X}=C(\overline{\Omega})\cap C_{b}(\Omega_{c}^{*})$ below for simplicity. Define the supremum norm on $\mathcal{X}$ by

\|w\|=\sup_{\bm{x}\in\widehat{\Omega}}|w(\bm{x})|,\quad\forall\,w\in\mathcal{X},

then $(\mathcal{X},\|\cdot\|)$ becomes a Banach space. According to the definition (4), $\mathcal{X}$ is well-defined under Cases (D1) and (D2). Let $f:C_{b}(\Omega^{*})\to C_{b}(\Omega^{*})$ be a nonlinear operator and ${\mathcal{L}}:D({\mathcal{L}})\to C_{b}(\Omega^{*})$ be a linear operator with the domain $D({\mathcal{L}})\subset\mathcal{X}$ . Then we consider the following two cases:

(C1)

$X=\{w|_{\Omega^{*}}\,|\,w\in\widehat{X}\}$ and $\partial X=C_{b}(\Omega_{c}^{*})$ , where $\widehat{X}=\{w\in\mathcal{X}\,|\,w|_{\Omega_{c}^{*}}=0\}$ ;
(C2)

$X=\{w|_{\overline{\Omega}_{+}}\,|\,w\in\widehat{X}\}$ , where $\widehat{X}=C_{\text{\rm per}}{(\overline{\Omega})}:=\{w\in C(\mathbb{R}^{d})$ is periodic with respect to the rectangle $\Omega\}$ .

It is easy to see that $X$ defined in any of Cases (C1) and (C2) can be regarded as a linear subspace of $\mathcal{X}$ in the sense of isometric isomorphism ( $X\simeq\widehat{X}$ ) through the zero or periodic extension to $\Omega_{c}^{*}$ or $\mathbb{R}^{d}$ correspondingly. We always omit the extension mapping between $X$ and $\widehat{X}$ for simplicity when there is no ambiguity. We also remark that $X$ equipped with the supremum norm $\|\cdot\|$ is a Banach space in both Cases (C1) and (C2), so is $\partial X$ in Case (C1).

For Case (D1), define $D({\mathcal{L}}_{0})=\{w\in D({\mathcal{L}})\cap X\,|\,{\mathcal{L}}w\in X\}$ which is a subspace of $X$ ; for Case (D2), define $D({\mathcal{L}}_{0})=X$ . Then we define the linear operator ${\mathcal{L}}_{0}={\mathcal{L}}|_{D({\mathcal{L}}_{0})}:D({\mathcal{L}}_{0})\to X$ . Note that for any $w\in D({\mathcal{L}}_{0})$ , it holds ${\mathcal{L}}_{0}w={\mathcal{L}}w$ in $\Omega^{*}$ , and ${\mathcal{L}}_{0}w$ is also well-defined in $\Omega^{*}_{c}$ in the sense of isomorphism.

The model problem we consider in this paper is a class of semilinear parabolic equations taking the following form

(5)

u_{t}={\mathcal{L}}u+f[u],\quad t>0,\ \bm{x}\in\Omega^{*},

where $u:[0,\infty)\times\widehat{\Omega}\to\mathbb{R}$ is the unknown function subject to the initial condition

(6)

u(0,\bm{x})=u_{0}(\bm{x}),\quad\bm{x}\in\widehat{\Omega}

with $u_{0}\in\mathcal{X}$ and either the periodic boundary condition for Case (C2) or the Dirichlet boundary condition

(7)

u(t,\bm{x})=g(t,\bm{x}),\quad t\geq 0,\ \bm{x}\in\Omega_{c}^{*}

for Case (C1) with $g\in C([0,\infty);\partial X)$ . The compatibility condition is also assumed to hold, that is, $g(0,\cdot)=u_{0}$ on $\Omega^{*}_{c}$ for Case (C1) and $u_{0}$ is $\Omega$ -periodic for Case (C2).

In order to establish the maximum bound principle (MBP) for the model problem (5), as well as its time discretizations proposed later, we make the following specific assumptions regarding the operators ${\mathcal{L}}$ , ${\mathcal{L}}_{0}$ and $f$ defined above.

Assumption 1 (Requirements on ${\mathcal{L}}$ and ${\mathcal{L}}_{0}$ ).

The linear operators ${\mathcal{L}}$ and ${\mathcal{L}}_{0}$ satisfy the following conditions:

(a) for any $w\in D({\mathcal{L}})$ and $\bm{x}_{0}\in\Omega^{*}$ , if

(8)

w(\bm{x}_{0})=\sup_{\bm{x}\in\widehat{\Omega}}w(\bm{x}),

then ${\mathcal{L}}w(\bm{x}_{0})\leq 0$ ;

(b) the domain $D({\mathcal{L}}_{0})$ is dense in $X$ ;

(c) there exists $\lambda_{0}>0$ such that the operator $\lambda_{0}\mathcal{I}-{\mathcal{L}}_{0}:D({\mathcal{L}}_{0})\to X$ is surjective, where $\mathcal{I}$ is the identity operator.

Note that it can be derived from Assumption 1-(a) that ${\mathcal{L}}$ maps any constant function in $D({\mathcal{L}})$ to the zero element in $\mathcal{X}$ .

Assumption 2 (Requirements on $f$ ).

The nonlinear operator $f$ acts as a composite function induced by a given one-variable continuously differentiable function $f_{0}:\mathbb{R}\to\mathbb{R}$ , that is,

(9)

f[w](\bm{x})=f_{0}(w(\bm{x})),\quad\forall\,w\in C_{b}(\Omega^{*}),\ \forall\,\bm{x}\in{\Omega^{*}},

and there exists a constant $\beta>0$ such that

(10)

f_{0}(\beta)\leq 0\leq f_{0}(-\beta).

Remark 1.

A more general case is that the function $f_{0}$ satisfies $f_{0}(M)\leq 0\leq f_{0}(m)$ for some $M>m$ instead of (10) in Assumption 2. In this case, one can carry out an affine transform to the unknown function $u$ . More precisely, take the affine transform $\eta:\mathbb{R}\to\mathbb{R}$ defined by

\eta(\theta)=\frac{M-m}{2\beta}\theta+\frac{M+m}{2},\quad\theta\in\mathbb{R},

and define $\tilde{f}_{0}(\theta)=\frac{2\beta}{M-m}f_{0}(\eta(\theta))$ for any $\theta\in\mathbb{R}$ , then it holds $\tilde{f}_{0}(\beta)\leq 0\leq\tilde{f}_{0}(-\beta)$ . By letting $\tilde{u}=\eta^{-1}(u)$ , we can obtain from (5) that

\tilde{u}_{t}={\mathcal{L}}\tilde{u}+\tilde{f}[\tilde{u}],\quad t>0,\ \bm{x}\in\Omega^{*},

where the nonlinear mapping $\tilde{f}$ is determined by the composite function

\tilde{f}[w](\bm{x})=\tilde{f}_{0}(w(\bm{x})),\quad\forall\,w\in C_{b}(\Omega^{*}),\ \forall\,\bm{x}\in{\Omega^{*}},

and satisfies Assumption 2.

Lemma 1.

Under Assumption 1, the following properties on ${\mathcal{L}}_{0}$ hold:

(i) ${\mathcal{L}}_{0}$ is dissipative, i.e., for any $\lambda>0$ and any $w\in D({\mathcal{L}}_{0})$ ,

\|(\lambda\mathcal{I}-{\mathcal{L}}_{0})w\|\geq\lambda\|w\|;

(ii) ${\mathcal{L}}_{0}$ is the generator of a contraction semigroup $\{S_{{\mathcal{L}}_{0}}(t)\}_{t\geq 0}$ on $X$ , i.e., $|\!|\!|S_{{\mathcal{L}}_{0}}(t)|\!|\!|\leq 1$ , where $|\!|\!|\cdot|\!|\!|$ denotes the operator norm defined by

|\!|\!|\mathcal{T}|\!|\!|=\sup_{\begin{subarray}{c}w\in X,\|w\|=1\end{subarray}}\|\mathcal{T}w\|.

Proof.

We first prove the result for Case (C1). For any $w\in D({\mathcal{L}}_{0})$ , since $w|_{\overline{\Omega}}\in C(\overline{\Omega})$ and $w|_{\Omega_{c}^{*}}=0$ , it is clear that $|w(\bm{x})|$ reaches its maximum at some point $\bm{x}_{0}\in\widehat{\Omega}$ , namely, $|w(\bm{x}_{0})|=\|w\|$ . Without loss of generality, let us assume $w(\bm{x}_{0})\geq 0$ ; otherwise, we consider $-w$ instead. If $\bm{x}_{0}\in\Omega^{*}$ , we know that ${\mathcal{L}}_{0}w(\bm{x}_{0})={\mathcal{L}}w(\bm{x}_{0})\leq 0$ by Assumption 1-(a). If $\bm{x}_{0}\in\Omega_{c}^{*}$ , then ${\mathcal{L}}_{0}w(\bm{x}_{0})=0$ by the definition of $X$ . For any $\lambda>0$ , we then have

\|(\lambda\mathcal{I}-{\mathcal{L}}_{0})w\|\geq|\lambda w(\bm{x}_{0})-{\mathcal{L}}_{0}w(\bm{x}_{0})|=\lambda w(\bm{x}_{0})-{\mathcal{L}}_{0}w(\bm{x}_{0})\geq\lambda w(\bm{x}_{0})=\lambda\|w\|,

which completes the proof of (i). Then, according to (i), Assumptions 1-(b) and 1-(c), the property (ii) follows from the Lumer–Phillips Theorem [27, Theorem II.3.15].

Now we consider Case (C2) and $\bm{x}_{0}\in\widehat{\Omega}$ such that $|w(\bm{x}_{0})|=\|w\|$ . If $\bm{x}_{0}\in\Omega_{c}^{*}$ , by the $\Omega$ -periodicity, we can regard $\bm{x}_{0}$ as a point in $\Omega^{*}$ . Then, the above analysis could be similarly done to obtain (i) and (ii). ∎

Remark 2.

The proof of Lemma 1-(i) uses the Assumption 1-(a) (which is only related to ${\mathcal{L}}$ ) and the definitions of $X$ and $D({\mathcal{L}}_{0})$ to deduce the fact that ${\mathcal{L}}_{0}w(\bm{x}_{0})\leq 0$ if (8) holds with $\bm{x}_{0}\in\widehat{\Omega}$ . Lemma 1-(ii) is the consequence of Lemma 1-(i) and Assumption 1-(b) and (c). We note that Lemma 1-(ii) is the key result to be used in later discussions. It can be established under different assumptions by using other approaches, see discussions in Example 2.6.

Next, let us introduce a stabilizing constant $\kappa>0$ and rewrite the equation (5) in the following equivalent form:

(11)

u_{t}+\kappa u={\mathcal{L}}u+\mathcal{N}[u],\quad t>0,\ \bm{x}\in\Omega^{*},

where $\mathcal{N}=\kappa\mathcal{I}+f$ . According to (9) in Assumption 2, we know

\mathcal{N}[w](\bm{x})=N_{0}(w(\bm{x})),\quad\forall\,w\in C_{b}(\Omega^{*}),\ \forall\,\bm{x}\in{\Omega^{*}},

where

(12)

N_{0}(\xi)=\kappa\xi+f_{0}(\xi),\quad\xi\in\mathbb{R}.

We impose a requirement on the selection of the stabilizing constant $\kappa$ as

(13)

\kappa\geq\max_{|\xi|\leq\beta}|f_{0}^{\prime}(\xi)|.

Note that (13) always can be reached since $f_{0}$ is continuously differentiable.

Lemma 2.

Under Assumption 2 and the requirement (13), it holds that

(i) $|N_{0}(\xi)|\leq\kappa\beta$ for any $\xi\in[-\beta,\beta]$ ;

(ii) $|N_{0}(\xi_{1})-N_{0}(\xi_{2})|\leq 2\kappa|\xi_{1}-\xi_{2}|$ for any $\xi_{1},\xi_{2}\in[-\beta,\beta]$ .

Proof.

From (12), we have $N_{0}^{\prime}(\xi)=\kappa+f_{0}^{\prime}(\xi)$ and

0\leq N_{0}^{\prime}(\xi)\leq 2\kappa,\quad\forall\,\xi\in[-\beta,\beta],

which leads to (ii) directly. From Assumption 2, for any $\xi\in[-\beta,\beta]$ , we know that

-\kappa\beta\leq-\kappa\beta+f_{0}(-\beta)=N_{0}(-\beta)\leq N_{0}(\xi)\leq N_{0}(\beta)=\kappa\beta+f_{0}(\beta)\leq\kappa\beta,

which gives (i). ∎

Now we can show that the aforementioned model problem admits a unique solution and possesses the MBP. Our main theorem can be stated as follows.

Theorem 3.

Given any constant $T>0$ . Under Assumptions 1 and 2, if

(14a)		$\|u_{0}(\bm{x})\|\leq\beta,\quad\forall\,\bm{x}\in\widehat{\Omega},$
then the equation (5) subject to the initial condition (6) and either the periodic boundary condition or the Dirichlet boundary condition (7) with
(14b)		$\|g(t,\bm{x})\|\leq\beta,\quad\forall\,t\in[0,T],\ \forall\,\bm{x}\in\Omega_{c}^{*}$

has a unique solution $u\in C([0,T];\mathcal{X})$ and it satisfies $\|u(t)\|\leq\beta$ for any $t\in[0,T]$ .

Proof.

Let us consider Case (C1) first. Denote $\mathcal{X}_{\beta}=\{w\in\mathcal{X}\,|\,\|w\|\leq\beta\}$ and $C_{g}([0,t];\mathcal{X}_{\beta})=\{w\in C([0,t];\mathcal{X}_{\beta})\,|\,w|_{[0,t]\times\Omega_{c}^{*}}=g\}$ for any $t\in(0,T]$ . For a fixed $t_{1}>0$ and a given $v\in C_{g}([0,t_{1}];\mathcal{X}_{\beta})$ , let us define $w:[0,t_{1}]\to\mathcal{X}$ to be the solution of the following linear problem:

(15)

\begin{dcases}w_{t}+\kappa w={\mathcal{L}}w+\mathcal{N}[v],&t\in(0,t_{1}],\ \bm{x}\in\Omega^{*},\\ w(t,\bm{x})=g(t,\bm{x}),&t\in[0,t_{1}],\ \bm{x}\in\Omega_{c}^{*},\\ w(0,\bm{x})=u_{0}(\bm{x}),&\bm{x}\in\widehat{\Omega},\end{dcases}

where the constant $\kappa$ satisfies (13). It is easy to know that $w$ is uniquely defined on $[0,t_{1}]\times\widehat{\Omega}$ by noticing that setting $\mathcal{N}[v]=0$ , $g=0$ , and $u_{0}=0$ in (15) leads to $w(t)\in D({\mathcal{L}}_{0})$ , and thus $w(t)=\mathrm{e}^{-\kappa t}S_{{\mathcal{L}}_{0}}(t)u_{0}=0$ for each $t\in[0,t_{1}]$ . Suppose there exists $(t^{*},\bm{x}^{*})\in(0,t_{1}]\times\widehat{\Omega}$ such that $w(t^{*},\bm{x}^{*})=\max_{0\leq t\leq t_{1}}\|w(t)\|$ . By (14b), we only need to consider the case $\bm{x}^{*}\in\Omega^{*}$ . According to Assumption 1-(a), we have $w_{t}\geq 0$ and ${\mathcal{L}}w\leq 0$ at $(t^{*},\bm{x}^{*})$ , which implies $\kappa w(t^{*},\bm{x}^{*})\leq N_{0}(v(t^{*},\bm{x}^{*}))$ . Since $|v(t^{*},\bm{x}^{*})|\leq\beta$ , according to Lemma 2-(i), we obtain $w(t^{*},\bm{x}^{*})\leq\beta$ . Similarly, if there exists $(t^{**},\bm{x}^{**})\in(0,t_{1}]\times\widehat{\Omega}$ such that $w(t^{**},\bm{x}^{**})=-\max_{0\leq t\leq t_{1}}\|w(t)\|$ , one can show $w(t^{**},\bm{x}^{**})\geq-\beta$ . Thus we obtain $\|w(t)\|\leq\beta$ for any $t\in[0,t_{1}]$ , which means $w\in C_{g}([0,t_{1}];\mathcal{X}_{\beta})$ .

Next we define a mapping $\mathcal{A}:C_{g}([0,t_{1}];\mathcal{X}_{\beta})\to C_{g}([0,t_{1}];\mathcal{X}_{\beta})$ by $\mathcal{A}[v]=w$ through (15). For $v,\tilde{v}\in C_{g}([0,t_{1}];\mathcal{X}_{\beta})$ , we see that $w=\mathcal{A}[v]$ and $\tilde{w}=\mathcal{A}[\tilde{v}]$ satisfy

w(t)-\tilde{w}(t)=\int_{0}^{t}\mathrm{e}^{-\kappa(t-s)}S_{{\mathcal{L}}_{0}}(t-s)\big{\{}\mathcal{N}[v(s)]-\mathcal{N}[\tilde{v}(s)]\big{\}}\,\mathrm{d}s,\quad t\in[0,t_{1}].

Using Lemma 2-(ii), we can estimate the difference between $w(t)$ and $\tilde{w}(t)$ by

	$\displaystyle\\|w(t)-\tilde{w}(t)\\|$	$\displaystyle\leq\int_{0}^{t}\mathrm{e}^{-\kappa(t-s)}\|\!\|\!\|S_{{\mathcal{L}}_{0}}(t-s)\|\!\|\!\|\\|\mathcal{N}[v(s)]-\mathcal{N}[\tilde{v}(s)]\\|\,\mathrm{d}s$
		$\displaystyle\leq 2\kappa\int_{0}^{t}\mathrm{e}^{-\kappa(t-s)}\\|v(s)-\tilde{v}(s)\\|\,\mathrm{d}s$
		$\displaystyle\leq 2(1-\mathrm{e}^{-\kappa t_{1}})\\|v-\tilde{v}\\|_{C([0,t_{1}];\mathcal{X})},\quad t\in[0,t_{1}],$

then,

\|\mathcal{A}[v]-\mathcal{A}[\tilde{v}]\|_{C([0,t_{1}];\mathcal{X})}=\|w-\tilde{w}\|_{C([0,t_{1}];\mathcal{X})}\leq 2(1-\mathrm{e}^{-\kappa t_{1}})\|v-\tilde{v}\|_{C([0,t_{1}];\mathcal{X})}.

When $t_{1}<\kappa^{-1}\ln 2$ , we see that $\mathcal{A}$ is a contraction. Since $\mathcal{X}_{\beta}$ is closed in $\mathcal{X}$ , we know that $C_{g}([0,t_{1}];\mathcal{X}_{\beta})$ is complete with respect to the metric reduced by the norm $\|\cdot\|_{C([0,t_{1}];\mathcal{X})}$ . Then we can apply Banach’s fixed-point theorem to get the existence of a unique solution $u\in C_{g}([0,t_{1}];\mathcal{X}_{\beta})$ to the model equation (5) on $[0,t_{1}]$ . Using standard continuation argument, we then have the global existence of the unique solution $u\in C_{g}([0,T];\mathcal{X}_{\beta})$ for any $T>0$ .

For Case (C2), the system (15) still holds after removing the second equation in it. For the case $\bm{x}^{*}\in\Omega_{c}^{*}$ , by the periodicity, we could regard $\bm{x}^{*}$ as a point in $\Omega^{*}$ . Then, the analysis above is still valid. Putting all the different cases together, we get a complete proof. ∎

Remark 3.

We note that Theorem 3 also could be proved using other classical methods, such as the method of upper and lower solutions [81] in the case of scalar second-order elliptic operator. We emphasize that, in the proof given here, the MBP of the model equation (5) has no explicit dependence on the constant $\kappa$ , even though we have assumed that $\kappa$ satisfies (13) in order to use the Lemma 2. Meanwhile, we will see later that the choice of the constant $\kappa$ plays an important role in designing MBP-preserving ETD schemes.

In the following subsections, we present some examples of the nonlinear and linear operators which are applicable to the mathematical framework established above.

2.2 Examples of the nonlinear function $f_{0}$

One usually chooses $f_{0}$ as $f_{0}=-F^{\prime}$ with $F$ being a primitive function, when considering the phase-field models, or says, the gradient flows of some free energy functionals.

Example 2.1.

Consider the function

(16)

f_{0}(s)=\lambda s(1-s^{p}),

where $\lambda>0$ and $p\in\mathbb{N}_{+}$ . According to Remark 1, $f_{0}$ satisfies that $f_{0}(M)\leq 0\leq f_{0}(m)$ with any $M\geq 1$ and $1\geq m\geq 0$ . Especially, for an even integer $p$ , one can choose $\beta\geq 1$ such that $f_{0}$ satisfies Assumption 2, then the requirement (13) becomes $\kappa\geq\lambda[(p+1)\beta^{p}-1]$ . The special case $p=1$ gives the well-known logistic function and the case $p=2$ with $\lambda=1$ gives

(17)

f_{0}(s)=s-s^{3},

the derivative of $-F$ with $F(s)=\frac{1}{4}(s^{2}-1)^{2}$ (the quartic double-well potential).

Example 2.2.

Consider the Flory–Huggins free energy

F(s)=\frac{\theta}{2}[(1+s)\ln(1+s)+(1-s)\ln(1-s)]-\frac{\theta_{c}}{2}s^{2},

where $\theta$ and $\theta_{c}$ are two constants satisfying $0<\theta<\theta_{c}$ , and

(18)

f_{0}(s)=-F^{\prime}(s)=\frac{\theta}{2}\ln\frac{1-s}{1+s}+\theta_{c}s.

Denoting by $\rho$ the positive root of $f_{0}(\rho)=0$ , i.e.,

(19)

\frac{1}{2\rho}\ln\frac{1+\rho}{1-\rho}=\frac{\theta_{c}}{\theta},

we can find that $\rho\in\Big{(}\sqrt{1-\frac{\theta}{2\theta_{c}-\theta}},1\Big{)}$ . Here, $f_{0}$ satisfies the condition in Assumption 2 with $\beta\in[\rho,1)$ , and then the requirement (13) becomes $\kappa\geq\frac{\theta}{1-\beta^{2}}-\theta_{c}$ .

Example 2.3.

The Peng–Robinson equation of state [82] is widely used in the oil industries and petroleum engineering. The Helmholtz free-energy density considered in such a model can be expressed as

F(s)=RTs(\ln s-1)-RTs\ln(1-bs)+\frac{as}{2\sqrt{2}b}\ln\frac{1+(1-\sqrt{2})bs}{1+(1+\sqrt{2})bs},

where $R$ is the universal gas constant, $T$ the temperature, $a=a(T)$ the energy parameter, and $b$ the covolume parameter. The values of these parameters could be found in [84]. Then

(20)

f_{0}(s)=-RT\ln\frac{s}{1-bs}-\frac{RTbs}{1-bs}-\frac{a}{2\sqrt{2}b}\ln\frac{1+(1-\sqrt{2})bs}{1+(1+\sqrt{2})bs}+\frac{as}{1+2bs-b^{2}s^{2}},

which has two zero points $m$ and $M$ satisfying $0<m<M<1/b$ . This example falls into the situation discussed in Remark 1.

2.3 Examples of the linear operator ${\mathcal{L}}$

2.3.1 Infinite dimensional examples

Here we present some examples of the linear operator ${\mathcal{L}}$ in the form of classic differential or nonlocal integral operators corresponding to Case (D1). We will verify that ${\mathcal{L}}$ and the induced operator ${\mathcal{L}}_{0}$ satisfy Assumption 1.

Example 2.4 (Second-order elliptic differential operator).

Consider the linear operator

(21)

{\mathcal{L}}w(\bm{x})=A(\bm{x}):\nabla^{2}w(\bm{x})+\bm{q}(\bm{x})\cdot\nabla w(\bm{x}),\quad\bm{x}\in\Omega,

where $\bm{q}\in C(\overline{\Omega};\mathbb{R}^{d})$ and $A\in C(\overline{\Omega};\mathbb{R}^{d\times d})$ is symmetric and positive definite uniformly (i.e., there exists $\theta>0$ such that $\bm{\xi}^{T}A(\bm{x})\bm{\xi}\geq\theta\bm{\xi}^{T}\bm{\xi}$ for all $\bm{x}\in\overline{\Omega}$ and $\bm{\xi}\in\mathbb{R}^{d}$ ). In this case, $\Omega^{*}=\Omega$ , $\Omega_{c}^{*}=\partial\Omega$ , $\widehat{\Omega}=\overline{\Omega}$ , the boundary condition (7) is the classic Dirichlet one, $\mathcal{X}=C(\overline{\Omega})$ , and $D({\mathcal{L}})=C^{2}(\overline{\Omega})$ . For any $w\in D({\mathcal{L}})$ and $\bm{x}_{0}\in\Omega$ such that (8) holds, we always have $\nabla w(\bm{x}_{0})=\bm{0}$ and $\nabla^{2}w(\bm{x}_{0})$ is negative semi-definite. Since $A(\bm{x}_{0})$ is symmetric, there exists an orthogonal matrix $O\in\mathbb{R}^{d\times d}$ such that

OA(\bm{x}_{0})O^{T}=\mathop{\operator@font diag}\nolimits\{\lambda_{1},\lambda_{2},\dots,\lambda_{d}\}=:\Lambda

with $\lambda_{i}\geq\theta$ for all $i=1,2,\dots,d$ . Thus, we have

A(\bm{x}_{0}):\nabla^{2}w(\bm{x}_{0})=(O^{T}\Lambda O):\nabla^{2}w(\bm{x}_{0})=\Lambda:(O^{T}\nabla^{2}w(\bm{x}_{0})O).

Since $O^{T}\nabla^{2}w(\bm{x}_{0})O$ is symmetric and negative semi-definite, its diagonal entries are all non-positive. Thus, we obtain $A(\bm{x}_{0}):\nabla^{2}w(\bm{x}_{0})\leq 0$ and then ${\mathcal{L}}w(\bm{x}_{0})\leq 0$ , which gives Assumption 1-(a). Assumption 1-(b) is guaranteed by the fact that $X=C_{0}(\overline{\Omega})$ and $D({\mathcal{L}}_{0})=\{w\in C^{2}(\overline{\Omega})\cap C_{0}(\overline{\Omega})\,|\,({\mathcal{L}}w)|_{\partial\Omega}=0\}$ for Case (C1), and $X=C_{\text{\rm per}}(\overline{\Omega})$ and $D({\mathcal{L}}_{0})=C^{2}(\overline{\Omega})\cap C_{\text{\rm per}}(\overline{\Omega})$ for Case (C2) in the sense of isomorphism. Assumption 1-(c) can be obtained from the standard analysis of the existence and regularity of the solution to the partial differential equation $\lambda w-{\mathcal{L}}_{0}w=f$ in $\Omega$ for any $\lambda>0$ and $f\in X$ .

Remark 4.

The second-order elliptic differential operator could be also defined in the divergence form:

(22)

\widetilde{{\mathcal{L}}}w(\bm{x})=\nabla\cdot(A(\bm{x})\nabla w(\bm{x}))+\widetilde{\bm{q}}(\bm{x})\cdot\nabla w(\bm{x}),\quad\bm{x}\in\Omega

with $A\in C(\overline{\Omega};\mathbb{R}^{d\times d})\cap C^{1}(\Omega;\mathbb{R}^{d\times d})$ and $\widetilde{\bm{q}}\in C(\overline{\Omega};\mathbb{R}^{d})$ . This form could be written in a non-divergence form (21) by setting $\bm{q}(\bm{x})=\nabla\cdot A(\bm{x})+\widetilde{\bm{q}}(\bm{x})$ . Similarly we can show the operator $\widetilde{{\mathcal{L}}}$ and $\widetilde{{\mathcal{L}}}_{0}$ also satisfies Assumption 1.

Example 2.5 (Nonlocal diffusion operator).

Consider the integral operator [17]

(23)

{\mathcal{L}}w(\bm{x})=\int_{\widehat{\Omega}}\gamma(\bm{x},\bm{y})(w(\bm{y})-w(\bm{x}))\,\mathrm{d}\bm{y},\quad\bm{x}\in\overline{\Omega},

where $\gamma:\widehat{\Omega}\times\widehat{\Omega}\to\mathbb{R}$ is a symmetric nonnegative kernel function, i.e., $\gamma(\bm{x},\bm{y})=\gamma(\bm{y},\bm{x})\geq 0$ . We consider the widely studied case of the nonlocal operator (23), in which the kernel is radial and parameterized by a horizon parameter $\delta>0$ much less than the size of $\Omega$ , i.e., $\gamma(\bm{x},\bm{y})=\gamma_{\delta}(|\bm{x}-\bm{y}|)$ for some nonnegative function $\gamma_{\delta}:\mathbb{R}\to\mathbb{R}$ with $\gamma_{\delta}|_{\mathbb{R}\setminus(0,\delta]}=0$ and $\gamma_{\delta}|_{(0,\delta/2]}\geq\gamma^{*}$ for some constant $\gamma^{*}>0$ . In this case, $\Omega^{*}=\overline{\Omega}$ , $\Omega_{c}^{*}=\Omega_{\delta}:=\{\bm{y}\in\mathbb{R}^{d}\setminus\overline{\Omega}\,|\,\exists\,\bm{x}\in\overline{\Omega}\text{ such that }|\bm{x}-\bm{y}|\leq\delta\}$ , $\widehat{\Omega}=\overline{\Omega}\cup\Omega_{{\delta}}$ , the boundary condition (7) becomes a volume constraint, and $\mathcal{X}=C(\overline{\Omega})\cap C_{b}(\Omega_{\delta})$ . Then, the operator (23) can be re-expressed as

(24)		$\displaystyle{\mathcal{L}}w(\bm{x})$	$\displaystyle=$	$\displaystyle\int_{B_{\delta}(\bm{\bm{x}})}\gamma_{\delta}(\|\bm{x}-\bm{y}\|)\big{(}w(\bm{y})-w(\bm{x})\big{)}\,\mathrm{d}\bm{y}$
(24)			$\displaystyle=$	$\displaystyle\frac{1}{2}\int_{B_{\delta}(\bm{0})}\gamma_{\delta}(\|\bm{y}\|)\big{(}w(\bm{x}+\bm{y})+w(\bm{x}-\bm{y})-2w(\bm{x})\big{)}\,\mathrm{d}\bm{y},\quad\bm{x}\in\overline{\Omega},$

For suitable $\gamma_{\delta}$ subject to the finite second-order moment condition:

\int_{B_{\delta}(\bm{0})}|\bm{y}|^{2}\gamma_{\delta}(|\bm{y}|)\,\mathrm{d}\bm{y}=2d,

the operator ${\mathcal{L}}$ is consistent with the standard Laplacian $\Delta$ as $\delta\to 0$ (see, e.g., [15, 17]). Assumption 1-(a) results from the nonnegativity of the kernel. Whether Assumptions 1-(b) and 1-(c) are satisfied depends on the choice of the kernel. Let us consider the situation of integrable kernels only (i.e., $\gamma_{\delta}(|\bm{y}|)\in L^{1}(\mathbb{R}^{d})$ ), for which $D({\mathcal{L}})=\mathcal{X}$ and ${\mathcal{L}}w(\bm{x})=(\gamma_{\delta}*w)(\bm{x})-\alpha_{\delta}w(\bm{x})$ for any $\bm{x}\in\overline{\Omega}$ . Here, $\gamma_{\delta}*w$ denotes the convolution and $\alpha_{\delta}=\|\gamma_{\delta}(|\cdot|)\|_{L^{1}(\mathbb{R}^{d})}$ . For Case (C1), we have $X=C(\overline{\Omega})$ and $D({\mathcal{L}}_{0})=C(\overline{\Omega})$ in the sense of isomorphism using zero extension to $\Omega_{\delta}$ thus Assumption 1-(b) follows automatically. For any $\lambda>0$ and $f\in C(\overline{\Omega})$ , there exists a unique solution $w\in L^{\infty}({\widehat{\Omega}})$ to the equation $\lambda w-{\mathcal{L}}_{0}w=f$ in $\overline{\Omega}$ , and it further can be shown $w\in C({\overline{\Omega}})$ [8, 23] by the property that the convolution of $L^{\infty}$ and $L^{1}$ functions is uniformly continuous, which verifies Assumption 1-(c). For Case (C2), we have $X=D({\mathcal{L}}_{0})=C_{\text{\rm per}}(\overline{\Omega})$ and all the assumptions are still satisfied. Regarding non-integrable kernels with strong singularity at the origin like the fractional power, we refer to discussion in the next example.

Example 2.6 (Fractional Laplace operator).

For a fixed $s\in(0,1)$ , the fractional Laplace operator $(-\Delta)^{s}$ is the pseudo-differential operator with symbol $|\bm{\xi}|^{2s}$ , that is,

\mathcal{F}[(-\Delta)^{s}w](\bm{\xi})=|\bm{\xi}|^{2s}\mathcal{F}[w](\bm{\xi}),\quad\bm{\xi}\in\mathbb{R}^{d},

where $w\in\mathscr{S}(\mathbb{R}^{d})$ , the class of Schwarz functions, and $\mathcal{F}$ denotes the Fourier transform. Denoting ${\mathcal{L}}=-(-\Delta)^{s}$ , an equivalent form of ${\mathcal{L}}$ on a bounded spatial domain $\overline{\Omega}\subset\mathbb{R}^{d}$ is given by [87]

(25)

{\mathcal{L}}w(\bm{x})={c_{d,s}}\int_{\mathbb{R}^{d}}\frac{w(\bm{y})-w(\bm{x})}{|\bm{y}|^{d+2s}}\,\mathrm{d}\bm{y},\quad\bm{x}\in\overline{\Omega}

with $c_{d,s}=\frac{2^{2s}\Gamma(d/2+s)}{\pi^{d/2}\Gamma(-s)}$ . In this case, $\Omega^{*}=\overline{\Omega}$ , $\Omega_{c}^{*}=\mathbb{R}^{d}\setminus\overline{\Omega}$ , $\widehat{\Omega}=\mathbb{R}^{d}$ and $\mathcal{X}=C(\overline{\Omega})\cap C_{b}(\mathbb{R}^{d}\setminus\overline{\Omega})$ . As shown in [96], the operator ${\mathcal{L}}$ could be regarded as the limit, when $\delta\to\infty$ , of the nonlocal operator ${\mathcal{L}}$ defined by (24) equipped with a special rescaled and truncated fractional kernel $\gamma_{\delta}(r)=c_{d,s}r^{-d-2s}$ for $r\in(0,\delta]$ . Let us consider Case (C1) only, for which $X=C(\overline{\Omega})$ in the sense of isomorphism using zero extension to $\mathbb{R}^{d}\setminus\overline{\Omega}$ . Define ${\mathcal{L}}_{0}$ as the restriction of ${\mathcal{L}}$ on $X$ . Due to the lack of full elliptic regularity, Lemma 1 is not readily applicable. However, we can still get Lemma 1-(ii) using the Beurling–Deny criteria [80]. Indeed, as proved in [36] if the boundary $\partial\Omega$ is $C^{1,1}$ and still valid if $\partial\Omega$ is Lipschitz, the equation $u_{t}-{\mathcal{L}}u=0$ in $(0,\infty)\times\overline{\Omega}$ subject to $u(t,\cdot)|_{\mathbb{R}^{d}\setminus\overline{\Omega}}=0$ with the initial data $u_{0}\in X$ has a unique weak solution $u(t,\cdot)\in C(\mathbb{R}^{d})\cap H^{s}(\mathbb{R}^{d})$ for all $t>0$ , i.e., $u(t,\cdot)|_{\mathbb{R}^{d}\setminus\overline{\Omega}}=0$ and

\int_{\Omega}u_{t}(t,\bm{x})v(t,\bm{x})\,\mathrm{d}\bm{x}+\frac{c_{d,s}}{2}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\frac{(u(t,\bm{x})-u(t,\bm{y}))(v(t,\bm{x})-v(t,\bm{y}))}{|\bm{x}-\bm{y}|^{d+2s}}\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}=0,

for any $v(t,\cdot)\in H^{s}(\mathbb{R}^{d})$ with $v(t,\cdot)|_{{\mathbb{R}^{d}\setminus\overline{\Omega}}}=0$ . To obtain Lemma 1-(ii), i.e., the solution $u(t)=S_{{\mathcal{L}}_{0}}(t)u_{0}$ induces a contraction semigroup $\{S_{{\mathcal{L}}_{0}}(t)\}_{t\geq 0}$ on $X$ with respect to the supremum norm $\|\cdot\|$ , it suffices to show that for any $u_{0}\in X$ with $u_{0}(\bm{x})\geq 0$ on $\overline{\Omega}$ , one has $u(t,\bm{x})\geq 0$ in $(0,\infty)\times\overline{\Omega}$ . The latter can be checked by taking $v=-u^{-}:=\min\{u,0\}$ in the above weak form to get that

	$\displaystyle\int_{\Omega}u_{t}(t,\bm{x})u^{-}(t,\bm{x})\,\mathrm{d}\bm{x}$	$\displaystyle=\frac{c_{d,s}}{2}\int_{u(t,\bm{x})<0}\int_{u(t,\bm{y})<0}\frac{(u(t,\bm{x})-u(t,\bm{y}))^{2}}{\|\bm{x}-\bm{y}\|^{d+2s}}\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}$
		$\displaystyle\quad+c_{d,s}\int_{u(t,\bm{x})<0}\int_{u(t,\bm{y})\geq 0}\frac{(u(t,\bm{x})-u(t,\bm{y}))u(t,\bm{x})}{\|\bm{x}-\bm{y}\|^{d+2s}}\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}\geq 0,$

which deduces

0\leq\int_{u(t,\bm{x})<0}u^{2}(t,\bm{x})\,\mathrm{d}\bm{x}\leq\int_{u_{0}(\bm{x})<0}u^{2}_{0}(\bm{x})\,\mathrm{d}\bm{x}=0,

thus $u(t,\bm{x})\geq 0$ in $(0,\infty)\times\overline{\Omega}$ . With Lemma 1-(ii) established, the first part of the proof of Theorem 3 can also be modified accordingly so that the same theorem remains valid.

To sum up, for the above examples, we can obtain the following result on the MBP according to Theorem 3.

Corollary 4.

Let the nonlinear function $f_{0}$ be given by either (16) with even $p$ or (18) and the operator ${\mathcal{L}}$ defined by any of (21), (24) and (25). If (14) holds, then the solution of the model problem (5) exists and satisfies $|u(t,\bm{x})|\leq\beta$ for any $t\geq 0$ and $\bm{x}\in\widehat{\Omega}$ , where $\beta=1$ for the case of (16) and $\beta=\rho$ for (18), respectively.

Remark 5.

When $f_{0}$ is given by either (16) with general $p$ or (20), based on Remark 1 and Corollary 4, we can derive a similar result on the MBP for the model equation (5) in the form of $m\leq u(t,\bm{x})\leq M$ for any $t\geq 0$ and $\bm{x}\in\widehat{\Omega}$ .

2.3.2 Finite dimensional examples

We consider some concrete linear operators for the situation of $\dim\mathcal{X}<\infty$ , in particular, the case of $(\Omega,\Omega_{c})$ is given by Case (D2). Let us denote by ${\mathcal{L}}_{h}$ the linear operator, with the subscript used to distinguish it from the notation in the infinite dimensional cases. Without loss of generality, we assume that $\Omega$ is the node set of a uniform Cartesian mesh on the hypercube $\widetilde{\Omega}=(0,1)^{d}$ . The choice of uniform mesh simplifies the notation but we note that the analysis can be extended to non-uniform and non-Cartesian grids as discussed later with respect to finite element (or finite volume) discretizations. With the uniform Cartesian mesh, for a given $M_{0}\in\mathbb{N}_{+}$ and the uniform mesh size $h=1/M_{0}$ , let $\Omega=\{\bm{x}_{\bm{i}}=h{\bm{i}}\,|\,{\bm{i}}\in\{1,2,\dots,M_{0}-1\}^{d}\}$ , $\overline{\Omega}=\{\bm{x}_{\bm{i}}=h{\bm{i}}\,|\,{\bm{i}}\in\{0,1,2,\dots,M_{0}\}^{d}\}$ , and $\overline{\Omega}_{+}=\{\bm{x}_{\bm{i}}=h{\bm{i}}\,|\,{\bm{i}}\in\{1,2,\dots,M_{0}\}^{d}\}$ . Thus, $X$ is isomorphic to $\mathbb{R}^{M^{d}}$ with the norm $\|\cdot\|$ equivalent to the vector infinity-norm. Here, for Case (C1), we let $M=M_{0}-1$ if $\widetilde{\Omega}_{c}=\partial\widetilde{\Omega}$ and $M=M_{0}+1$ otherwise, while $M=M_{0}$ for Case (C2). The mesh can be extended as needed to represent ${\widetilde{\Omega}}_{c}^{*}$ , and the solution on ${\Omega}_{c}^{*}$ (mesh nodes in $\widetilde{\Omega}_{c}^{*}$ ) is explicitly given for Case (C1) or can be obtained from that on $\overline{\Omega}_{+}$ by the periodicity of the solution in ${\mathbb{R}^{d}}$ for Case (C2).

We may view the grid function $w\in X$ as a vector $\bm{w}=(w_{1},w_{2},\dots,w_{M^{d}})^{T}\in\mathbb{R}^{M^{d}}$ with $w_{i}$ denoting the values of $w$ at nodes $\bm{x}_{\bm{i}}\in\Omega^{*}$ (for Case (C1)) or $\bm{x}_{\bm{i}}\in\overline{\Omega}_{+}$ (for Case (C2)) ordered lexicographically. Then, the nonlinear mapping $f$ plays a role as a vector-valued function $\mathbb{R}^{M^{d}}\to\mathbb{R}^{M^{d}}$ satisfying

(f[\bm{w}])_{i}=f_{0}(w_{i}),\quad\forall\,\bm{w}\in\mathbb{R}^{M^{d}},\ 1\leq i\leq M^{d}.

In addition, due to the discrete nature, we now have that $D({\mathcal{L}}_{h0})=X$ and the linear operator ${\mathcal{L}}_{h0}$ is equivalent to ${\mathcal{L}}_{h}|_{X}$ , which is an $M^{d}$ -by- $M^{d}$ matrix. Thus, the semigroup generated by ${\mathcal{L}}_{h0}$ is actually the matrix exponential $S_{{\mathcal{L}}_{h0}}(t)=\mathrm{e}^{t{\mathcal{L}}_{h0}}$ . In these settings, Assumption 1-(b) is trivial and Assumption 1-(c) could be verified more directly by the following result.

Proposition 5.

If a matrix ${\mathcal{L}}_{h0}=(a_{ij})_{M^{d}\times M^{d}}$ satisfies

(26)

|a_{ii}|\geq\sum_{\begin{subarray}{c}j=1\\ j\not=i\end{subarray}}^{M^{d}}|a_{ij}|,\quad a_{ii}<0,\quad a_{ij}\geq 0\ (j\not=i)

for all $\,i,j=1,2,\dots,M^{d}$ , then ${\mathcal{L}}_{h0}$ satisfies Assumption 1-(c).

Example 2.7 (Central difference operator for the Laplacian).

To make the discussion more clearly, let us begin with the one-dimensional case. Now, $\mathcal{X}\simeq\mathbb{R}^{M_{0}+1}$ . For Case (C1), we have $\partial X\simeq\mathbb{R}^{2}$ . For any $w\in\mathcal{X}$ , the second-order central difference approximation of the Laplacian is defined by

(27)

{\mathcal{L}}_{h}w(x_{i})=\frac{1}{h^{2}}\big{(}w(x_{i-1})-2w(x_{i})+w(x_{i+1})\big{)}

for $x_{i}\in\Omega$ with $w(x_{0})$ and $w(x_{M_{0}})$ determined by the boundary condition. For Case (C2), for a grid function $w\in\mathcal{X}$ , we still define ${\mathcal{L}}_{h}$ as (27) for $x_{i}\in\overline{\Omega}_{+}$ but with $w(x_{0})=w(x_{M_{0}})$ and $w(x_{M_{0}+1})=w(x_{1})$ due to the periodicity. Assumption 1-(a) holds according to (27). Define

D_{h}=\frac{1}{h^{2}}\begin{pmatrix}-2&1&~{}&~{}&c\\ 1&-2&1\\ ~{}&\ddots&\ddots&\ddots\\ ~{}&~{}&1&-2&1\\ c&~{}&~{}&1&-2\end{pmatrix}\in\mathbb{R}^{M\times M},\quad\text{with }c=\begin{cases}0,&\text{for Case (C1)},\\ 1,&\text{for Case (C2)},\end{cases}

and $D_{h}$ satisfies (26) obviously. We have ${\mathcal{L}}_{h0}=D_{h}$ , so Assumption 1-(c) is satisfied.

For the two- and three-dimensional cases, we define ${\mathcal{L}}_{h}$ as a summation of the central difference approximations in every coordinate direction given as (27). Thus, ${\mathcal{L}}_{h}$ satisfies Assumption 1-(a). Also, we can present the matrix ${\mathcal{L}}_{h0}$ by using the Kronecker products as

(28)		$\displaystyle{\mathcal{L}}_{h0}$	$\displaystyle=I_{M}\otimes D_{h}+D_{h}\otimes I_{M},\quad\;\qquad\qquad\qquad\qquad\qquad\qquad\quad\;\;\text{(2-D)}$
(29)		$\displaystyle{\mathcal{L}}_{h0}$	$\displaystyle=I_{M}\otimes I_{M}\otimes D_{h}+I_{M}\otimes D_{h}\otimes I_{M}+D_{h}\otimes I_{M}\otimes I_{M},\qquad\text{(3-D)}$

where $I_{M}$ is the $M\times M$ identity matrix, and the vectors determined by the boundary conditions could be given in the similar way. It is easy to check that ${\mathcal{L}}_{h0}$ corresponding to (28) or (29) satisfies Assumption 1-(c).

Remark 6.

It is well-known that (28) provides the central difference approximation of the Laplacian in the two-dimensional case. However, if one considers more general elliptic operators, such as (21), including mixed partial derivatives and a convection term, the corresponding discretized operator cannot be written in the form of Kronecker products. Instead, the upwind formulas are adequate to discretize the derivatives in the convection term and those in the mixed partial derivatives. The detailed analysis of the resulted linear operator is quite similar to Example 2.7 and we omit it. A general formulation for the finite difference approximation of (21) could be found in [9], where a sufficient and necessary condition (more general than Proposition 5) is given for the discrete MBP of the approximating operator.

Example 2.8 (Quadrature-based difference operator for the nonlocal diffusion).

Recalling the nonlocal diffusion operator (24), we consider its quadrature-based finite difference discretization, which is well known to be asymptotically compatible [22, 95]. The dimension of $\Omega_{c}$ clearly depends on the value of $\delta$ in this case. For a given $w\in\mathcal{X}$ , we define at $\bm{x}_{\bm{i}}\in\overline{\Omega}$ for Case (C1) or $\bm{x}_{\bm{i}}\in\overline{\Omega}_{+}$ for Case (C2) via

(30)

{\mathcal{L}}_{h}w(\bm{x}_{\bm{i}})=\sum_{\bm{0}\not=\bm{s}_{\bm{j}}\in B_{\delta}(\bm{0})}\frac{w(\bm{x}_{\bm{i}}+\bm{s}_{\bm{j}})+w(\bm{x}_{\bm{i}}-\bm{s}_{\bm{j}})-2w(\bm{x}_{\bm{i}})}{|\bm{s}_{\bm{j}}|^{2}}|\bm{s}_{\bm{j}}|_{1}\beta_{\delta}(\bm{s}_{\bm{j}}),

where the values of $w(\bm{y})$ with $\bm{y}\in\Omega_{c}^{*}$ could be supplementally defined using the boundary conditions so that the operator (30) is well-defined. Here, $|\cdot|_{1}$ stands for the vector $1$ -norm, and

\beta_{\delta}(\bm{s}_{\bm{j}})=\frac{1}{2}\int_{B_{\delta}(\bm{0})}\psi_{\bm{j}}(\bm{s})\frac{|\bm{s}|^{2}}{|\bm{s}|_{1}}\gamma_{\delta}(|\bm{s}|)\,\mathrm{d}\bm{s},

where $\psi_{\bm{j}}$ is the piecewise $d$ -multilinear basis function satisfying $\psi_{\bm{j}}(\bm{s}_{\bm{i}})=0$ when ${\bm{i}}\not={\bm{j}}$ and $\psi_{\bm{j}}(\bm{s}_{\bm{j}})=1$ . It is obvious that ${\mathcal{L}}_{h}$ given by (30) satisfies Assumption 1-(a). By ordering the nodes lexicographically, we regard the function $w\in X$ as a vector in $\mathbb{R}^{M^{d}}$ and it is easy to check that the matrix ${\mathcal{L}}_{h0}$ , the restriction of ${\mathcal{L}}_{h}$ on $X$ , is diagonally dominant in the sense of (26). Thus, Assumption 1-(c) is satisfied. Moreover, ${\mathcal{L}}_{h0}$ is a band matrix with the bandwidth depending on the value of $\delta$ .

Example 2.9 (Finite difference operator for the fractional Laplacian).

For the fractional Laplace operator (25), we consider its finite difference discretization in Case (C1), and for simplicity the homogeneous Dirichlet boundary condition is enforced so that the operators ${\mathcal{L}}_{h}$ and ${\mathcal{L}}_{h0}$ are essentially the same in this case. For the one-dimensional case, a second-order finite difference discretization [25, 52] (a limiting case of a similar scheme given in [95] for nonlocal diffusion models) is defined by:

(31)

{\mathcal{L}}_{h}w(x_{i})=\sum_{j=0}^{M_{0}}({\mathcal{L}}_{h0})_{ij}w(x_{j}),\quad x_{i}\in\overline{\Omega},

for any $w\in X$ , where the entries of the matrix ${\mathcal{L}}_{h0}$ are given by

\small({\mathcal{L}}_{h0})_{ij}=\frac{c_{1,s}}{\nu h^{2s}}\begin{dcases}-(2^{\nu}+\lfloor\gamma\rfloor-1)-\sum_{k=2}^{M_{0}-1}\frac{(k+1)^{\nu}-(k-1)^{\nu}}{k^{\gamma}}\\ \qquad\qquad\qquad\qquad-\frac{M_{0}^{\nu}-(M_{0}-1)^{\nu}}{M_{0}^{\gamma}}-\frac{\nu}{sM_{0}^{2s}},&j=i,\\ \frac{1}{2}(2^{\nu}+\lfloor\gamma\rfloor-1),&|j-i|=1,\\ \frac{(|j-i|+1)^{\nu}-(|j-i|-1)^{\nu}}{2|j-i|^{\gamma}},&2\leq|j-i|\leq M_{0}-1,\\ \frac{M_{0}^{\nu}-(M_{0}-1)^{\nu}}{2M_{0}^{\gamma}},&|j-i|=M_{0}\end{dcases}

for $0\leq i,j\leq M_{0}$ , where $\nu=(1-s)\lfloor\gamma\rfloor$ with $\gamma=2$ or $\gamma=1+s$ . Obviously, the matrix ${\mathcal{L}}_{h0}$ is diagonally dominant in the sense of (26), and thus Assumptions 1-(a) and 1-(c) are satisfied. For the two- and three-dimensional cases, the operators ${\mathcal{L}}_{h}$ and ${\mathcal{L}}_{h0}$ could be given in the similar way, we refer to [26, 52] for details, and it is easy to verify that Assumption 1 is still satisfied.

Example 2.10 (Finite element operator for the Laplacian).

Let us consider Case (C1) again. Let $\mathcal{T}_{h}=\{K\}$ be a uniform partition of $\widetilde{\Omega}$ into isosceles right triangles $K$ non-overlapping with each other based on the set of nodes $\overline{\Omega}$ . Let $V_{h}$ be the space of continuous and piecewise linear functions with respect to $\mathcal{T}_{h}$ and $V_{h}^{0}$ be its subspace with zero-trace:

V_{h}=\{w\in C(\overline{\widetilde{\Omega}})\,|\,\text{$w|_{K}$ is linear for any $K\in\mathcal{T}_{h}$}\},\quad V_{h}^{0}=\{w\in V_{h}\,|\,w|_{\partial\widetilde{\Omega}}=0\}.

We denote by $N_{I}=(M_{0}-1)^{d}$ and $N_{B}=(M_{0}+1)^{d}-N_{I}$ the numbers of the nodes in $\widetilde{\Omega}$ and on $\partial\widetilde{\Omega}$ respectively, and order the nodes such that $\Omega=\{\bm{x}_{j}\,|\,1\leq j\leq N_{I}\}$ consisting of all interior nodes and $\partial\Omega=\{\bm{x}_{j+N_{I}}\,|\,1\leq j\leq N_{B}\}$ all nodes on the boundary. Denote by $\{\phi_{j}\,|\,1\leq j\leq N_{I}+N_{B}\}$ the basis functions of $V_{h}$ satisfying $\phi_{i}(\bm{x}_{j})=\delta_{ij}$ for any $1\leq i,j\leq N_{I}+N_{B}$ . Define the operator ${\mathcal{L}}_{h}$ for $w\in\mathcal{X}\simeq V_{h}$ by

(32)

{\mathcal{L}}_{h}w(\bm{x}_{i})=\sum_{j=1}^{N_{I}+N_{B}}l_{ij}w(\bm{x}_{j}),\quad\bm{x}_{i}\in\Omega,

where

l_{ij}=-\frac{\int_{\widetilde{\Omega}}\nabla\phi_{i}\cdot\nabla\phi_{j}\,\mathrm{d}\bm{x}}{\int_{\widetilde{\Omega}}\phi_{i}\,\mathrm{d}\bm{x}},\quad 1\leq i\leq N_{I},\ 1\leq j\leq N_{I}+N_{B}.

For each $i=1,2,\dots,N_{I}$ , it is known that

(33)

\int_{\widetilde{\Omega}}\phi_{i}\,\mathrm{d}\bm{x}>0,\quad\int_{\widetilde{\Omega}}|\nabla\phi_{i}|^{2}\,\mathrm{d}\bm{x}>0,\quad\text{and}\quad\int_{\widetilde{\Omega}}\nabla\phi_{i}\cdot\nabla\phi_{j}\,\mathrm{d}\bm{x}\leq 0\ \text{for $j\not=i$}.

Noticing that $\sum\limits_{j=1}^{N_{I}+N_{B}}\phi_{j}(\bm{x})\equiv 1$ on $\overline{\widetilde{\Omega}}$ , we have that

(34)

\sum_{j=1}^{N_{I}+N_{B}}\int_{\widetilde{\Omega}}\nabla\phi_{i}\cdot\nabla\phi_{j}\,\mathrm{d}\bm{x}=\int_{\widetilde{\Omega}}\nabla\phi_{i}\cdot\nabla\bigg{(}\sum_{j=1}^{N_{I}+N_{B}}\phi_{j}\bigg{)}\,\mathrm{d}\bm{x}=0,\quad 1\leq i\leq N_{I}.

Combination of (33) and (34) yields

(35)

|l_{ii}|=\sum_{\begin{subarray}{c}j=1\\ j\not=i\end{subarray}}^{N_{I}+N_{B}}|l_{ij}|,\quad l_{ii}<0,\quad l_{ij}\geq 0\ (j\not=i)

for $1\leq i\leq N_{I}$ and $1\leq j\leq N_{I}+N_{B}$ , which implies Assumption 1-(a). The limitation of ${\mathcal{L}}_{h}$ on $X\simeq V_{h}^{0}$ is given by the matrix ${\mathcal{L}}_{h0}\in\mathbb{R}^{N_{I}\times N_{I}}$ with $({\mathcal{L}}_{h0})_{ij}=l_{ij}$ , $i,j=1,2,\dots,N_{I}$ . Then, by (35), ${\mathcal{L}}_{h0}$ satisfies (26) and Assumption 1-(c) holds.

Remark 7.

For a general partition of $\widetilde{\Omega}$ , the FEM-based linear operator could be also defined by (32) and still satisfies (26), see, e.g., [78] for discussions. In fact, it is well-known in the literature that the discrete maximum principle holds for piecewise linear finite element approximations of the Laplacian on general two-dimensional Delaunay triangulations for which the circumscribing circle of any Delaunay triangle does not contain any other vertices in its interior, see [21]. The brief note [7] analyzes a stabilized Galerkin approximation of the Laplacian guaranteeing the discrete maximum principle on arbitrary meshes. In addition, one can also apply the finite volume method (FVM) [66] to obtain the discretized operator of the Laplacian. The FVM-based linear operator possesses similar properties as those we have considered in the above.

Applying Theorem 3 to the cases with $f_{0}$ given by either Example 2.1 or 2.2 (or Example 2.3 when a transformation is first applied as explained in Remark 1) and ${\mathcal{L}}={\mathcal{L}}_{h}$ determined by any one of Examples 2.7–2.10, we obtain the following result on the MBP of (5) in the space-discrete version.

Corollary 6.

Let the nonlinear function $f_{0}$ be given by either (16) with even $p$ or (18) and the linear operator ${\mathcal{L}}={\mathcal{L}}_{h}$ defined by any of (27), (30), (31) and (32). If (14) holds, then the solution of the model problem (5) exists and satisfies $|u(t,\bm{x})|\leq\beta$ for any $t\geq 0$ and $\bm{x}\in\widehat{\Omega}$ , where $\beta=1$ for the case of (16) and $\beta=\rho$ for (18), respectively.

2.4 Discussion on homogeneous Neumann boundary conditions

In Section 2.1, the MBPs for the cases of Dirichlet and periodic boundary conditions are studied, let us now discuss the case of homogeneous Neumann boundary condition. For simplicity, we focus the analysis on the infinite dimensional operators. The key ingredient is to find a Neumann operator $\mathcal{N}_{c}$ such that the following Gaussian formula holds:

\int_{\Omega^{*}}{\mathcal{L}}u(\bm{x})\,\mathrm{d}\bm{x}=\int_{\Omega_{c}^{*}}\mathcal{N}_{c}u(\bm{x})\,\mathrm{d}\bm{x}.

The homogeneous Neumann boundary condition then could be given by

(36)

\mathcal{N}_{c}u(\bm{x})=0,\quad\bm{x}\in\Omega_{c}^{*}.

For the second-order elliptic operator in the divergence form (22), using the classic Gaussian formula and assuming $\widetilde{\bm{q}}$ is divergence-free, we obtain

\int_{\Omega}{\mathcal{L}}u(\bm{x})\,\mathrm{d}\bm{x}=\int_{\partial\Omega}(A(\bm{x})\nabla u(\bm{x})+\widetilde{\bm{q}}(\bm{x})u(\bm{x}))\cdot\bm{n}(\bm{x})\,\mathrm{d}\bm{x},

where $\bm{n}$ denotes the unit outward normal vector to $\partial\Omega$ . Thus, the Neumann operator $\mathcal{N}_{c}$ could be defined as

\mathcal{N}_{c}u(\bm{x})=A(\bm{x})\nabla u(\bm{x})\cdot\bm{n}(\bm{x})+u(\bm{x})\widetilde{\bm{q}}(\bm{x})\cdot\bm{n}(\bm{x}),\quad\bm{x}\in\partial\Omega,

and the corresponding constraint becomes the classic Robin boundary condition. Alternatively, if we apply the classic Gaussian formula only to the diffusion part ${\mathcal{L}}_{\text{diff}}u(\bm{x})=\nabla\cdot(A(\bm{x})\nabla u(\bm{x}))$ , then $\mathcal{N}_{c}$ could be also defined by

(37)

\mathcal{N}_{c}u(\bm{x})=A(\bm{x})\nabla u(\bm{x})\cdot\bm{n}(\bm{x}),\quad\bm{x}\in\partial\Omega,

and the corresponding constraint becomes the classic Neumann boundary condition.

For the nonlocal diffusion operator ${\mathcal{L}}$ defined in (23), by the symmetry of the kernel, we have

	$\displaystyle\int_{\widehat{\Omega}}\int_{\widehat{\Omega}}\gamma(\bm{x},\bm{y})(u(\bm{y})-u(\bm{x}))\,\mathrm{d}\bm{y}\mathrm{d}\bm{x}=\int_{\widehat{\Omega}}\int_{\widehat{\Omega}}\gamma(\bm{x},\bm{y})(u(\bm{x})-u(\bm{y}))\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}=0,$
	$\displaystyle\int_{\Omega_{c}^{}}\int_{\Omega_{c}^{}}\gamma(\bm{x},\bm{y})(u(\bm{y})-u(\bm{x}))\,\mathrm{d}\bm{y}\mathrm{d}\bm{x}=\int_{\Omega_{c}^{}}\int_{\Omega_{c}^{}}\gamma(\bm{x},\bm{y})(u(\bm{x})-u(\bm{y}))\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}=0.$

Then we have

\int_{\Omega^{*}}{\mathcal{L}}u(\bm{x})\,\mathrm{d}\bm{x}=\int_{\Omega_{c}^{*}}\mathcal{N}_{c}u(\bm{x})\,\mathrm{d}\bm{x},

where $\mathcal{N}_{c}$ is defined by

\mathcal{N}_{c}u(\bm{x})=\int_{\Omega}\gamma(\bm{x},\bm{y})(u(\bm{x})-u(\bm{y}))\,\mathrm{d}\bm{y},\quad\bm{x}\in\Omega_{c}^{*}.

This definition, different from those in [17, 18], can be found in [15] as a nonlocal extension of (37).

We next consider the model problem (5) subject to the homogeneous Neumann boundary condition (36). In order to show that the MBP is still satisfied for the above classic and nonlocal cases, we need to generate a contraction semigroup as before so that all the analysis conducted above could be similarly used. To this end, we still let $\mathcal{X}=C(\overline{\Omega})$ if $\Omega_{c}=\partial\Omega$ and $\mathcal{X}=C(\overline{\Omega})\cap C_{b}(\Omega_{c}^{*})$ otherwise, and Assumption 1-(a) is clearly true for the above two operators ${\mathcal{L}}$ . Then we need to define $X$ and ${\mathcal{L}}_{0}$ appropriately and verify whether all conditions related to them in Remark 2 hold.

For the classic elliptic differential operator case, without loss of generality, let us assume $A(\bm{x})\equiv I_{d}$ and $\bm{q}(\bm{x})\equiv\bm{0}$ in (21). For any $w\in C^{2}(\overline{\Omega})$ , we can directly extend the domain of $\Delta w$ from $\Omega$ to $\overline{\Omega}$ by continuity so that ${\mathcal{L}}=\Delta:C^{2}(\overline{\Omega})\to C(\overline{\Omega})$ . Let us define $X=C(\overline{\Omega})$ , $D({\mathcal{L}}_{0})=\{w\in C^{2}(\overline{\Omega})\,|\,(\nabla w\cdot\bm{n})|_{\partial\Omega}=0\}$ and ${\mathcal{L}}_{0}={\mathcal{L}}|_{D({\mathcal{L}}_{0})}$ , then Assumptions 1-(b) and 1-(c) obviously hold. For any $w\in D({\mathcal{L}}_{0})$ , assume that (8) is true for some $\bm{x}_{0}\in\partial\Omega$ . By using the boundary condition ( $\nabla w\cdot\bm{n})(\bm{x}_{0})=0$ and the reflective extension in classic analysis, we have that $\nabla w(\bm{x}_{0})=\bm{0}$ and the Hessian $\nabla^{2}w(\bm{x}_{0})$ is negative semi-definite, and consequently it holds ${\mathcal{L}}_{0}w(\bm{x}_{0})=\Delta w(\bm{x}_{0})\leq 0$ by using an orthogonal transformation of the coordinates. Thus, we can similarly show that Lemma 1-(i) is valid and consequently ${\mathcal{L}}_{0}$ generates a contraction semigroup $\{S_{{\mathcal{L}}_{0}}(t)\}_{t\geq 0}$ on $X$ (i.e., Lemma 1-(ii)).

Next we turn to the nonlocal diffusion operator case, where ${\mathcal{L}}$ is defined by (24) with the integrable kernel as discussed in the Example 2.5. For any $w\in C(\overline{\Omega})$ , there exists a unique $\widetilde{w}\in\mathcal{X}$ satisfying $\widetilde{w}|_{\overline{\Omega}}=w$ and $(\mathcal{N}_{c}\widetilde{w})|_{\Omega_{\delta}}=0$ ; in particular, we have

(38)

\widetilde{w}(\bm{x})=\dfrac{\int_{\Omega\cap B_{\delta}(\bm{x})}\gamma_{\delta}(|\bm{x}-\bm{y}|)w(\bm{y})\,\mathrm{d}\bm{y}}{\int_{\Omega\cap B_{\delta}(\bm{x})}\gamma_{\delta}(|\bm{x}-\bm{y}|)\,\mathrm{d}\bm{y}},\quad\bm{x}\in\Omega_{\delta}.

By defining $X=\{w|_{\overline{\Omega}}\,|\,w\in\widehat{X}\}$ with $\widehat{X}=\{w\in\mathcal{X}\,|\,(\mathcal{N}_{c}w)|_{\Omega_{\delta}}=0\}$ , we can regard $X=C(\overline{\Omega})$ as a linear subspace of $\mathcal{X}$ in the sense of isomorphism ( $X\simeq\widehat{X}$ ) using such Neumann extension. Assumption 1-(b) automatically comes from $D({\mathcal{L}}_{0})=C(\overline{\Omega})$ . For any $\lambda>0$ and $f\in C(\overline{\Omega})$ , it is also not hard to verify that the equation $\lambda w-{\mathcal{L}}_{0}w=f$ in $\overline{\Omega}$ has a unique solution in $C(\overline{\Omega})$ , similar to the discussion in the earlier example for the zero Dirichlet boundary condition, which implies Assumption 1-(c). For any $w\in C(\overline{\Omega})$ , assume that (8) holds for some $\bm{x}_{0}\in\Omega_{c}^{*}$ , then we know by (38) that the maximum of $w$ on $\widehat{\Omega}$ must also be attained on $\overline{\Omega}$ . Thus, according to the proof of Lemma 1-(i), we see that its conclusion still remains valid, and consequently we obtain Lemma 1-(ii).

3 Exponential time differencing for temporal approximation

In this section, we construct temporal approximations of the model equation (5) by using the exponential time differencing (ETD) method. We will begin with the equivalent form (11) to develop the MBP-preserving ETD schemes of first- and second-order, followed by the convergence analysis. For all discussions and theorems established in this section, we only focus on Case (C1) while the results are all valid for Case (C2) by removing the equations containing the boundary term $g$ and using the periodicity of the solution in the proofs.

3.1 ETD schemes and discrete MBPs

Given a time step size $\tau>0$ , we divide the time interval by $\{t_{n}=n\tau\}_{n\geq 0}$ . To establish the ETD schemes for solving the model problem (5), we focus the equivalent equation (11) on the interval $[t_{n},t_{n+1}]$ , or equivalently, $w(s,\bm{x})=u(t_{n}+s,\bm{x})$ satisfying the following system

(39)

\begin{dcases}w_{s}+\kappa w={\mathcal{L}}w+\mathcal{N}[w],&s\in(0,\tau],\ \bm{x}\in\Omega^{*},\\ w(s,\bm{x})=g(t_{n}+s,\bm{x}),&s\in[0,\tau],\ \bm{x}\in\Omega_{c}^{*},\\ w(0,\bm{x})=u(t_{n},\bm{x}),&\bm{x}\in\widehat{\Omega}.\end{dcases}

For the first-order (in time) approximation of (39), we set $\mathcal{N}[u(t_{n}+s)]\approx\mathcal{N}[u(t_{n})]$ to obtain the first-order ETD (ETD1) scheme: for $n\geq 0$ and given $v^{n}$ , find $v^{n+1}=w^{n}(\tau)$ with $w^{n}:[0,\tau]\to\mathcal{X}$ solving

(40)

\begin{dcases}w^{n}_{s}+\kappa w^{n}={\mathcal{L}}w^{n}+\mathcal{N}[v^{n}],&s\in(0,\tau],\ \bm{x}\in\Omega^{*},\\ w^{n}(s,\bm{x})=g(t_{n}+s,\bm{x}),&s\in[0,\tau],\ \bm{x}\in{\Omega_{c}^{*}},\\ w^{n}(0,\bm{x})=v^{n}(\bm{x}),&\bm{x}\in\widehat{\Omega},\end{dcases}

where $v^{n}$ represents an approximation of $u(t_{n})$ and $v^{0}=u_{0}$ is given. Similar to (15), it is easy to show that $w^{n}$ is uniquely defined on $[0,\tau]\times\widehat{\Omega}$ , thus $v^{n+1}$ is well-defined. We note that unlike the continuous-in-time case, different choices of $\kappa$ do lead to different discretization schemes. Indeed, ${\mathcal{L}}_{0}-\kappa\mathcal{I}$ serves as the generator of the semigroup $\{\mathrm{e}^{-\kappa t}S_{{\mathcal{L}}_{0}}(t)\}_{t\geq 0}$ to control the nonlinear term and then to stabilize the time discretizations. Thus, a suitable choice of $\kappa$ becomes very important in our design of ETD schemes.

Theorem 7 (Maximum bound principle of the ETD1 scheme).

Suppose that Assumptions 1 and 2, (13) and (14) hold. Then the ETD1 scheme preserves the discrete MBP unconditionally, i.e., for any time step size $\tau>0$ , the ETD1 solution satisfies $\|v^{n}\|\leq\beta$ for any $n\geq 0$ .

Proof.

Since $\|v^{0}\|\leq\beta$ , we just need to show that $\|v^{k}\|\leq\beta$ implies $\|v^{k+1}\|\leq\beta$ for any $k$ . We have $v^{k+1}=w^{k}(\tau)$ , where $w^{k}$ satisfies (40) with the superscript $n$ replaced by $k$ . Suppose there exists $(s^{*},\bm{x}^{*})\in(0,\tau]\times\widehat{\Omega}$ such that

(41)

w^{k}(s^{*},\bm{x}^{*})=\max_{0\leq s\leq\tau}\|w^{k}(s)\|.

Similar to the first part of the proof of Theorem 3, we can obtain $\kappa w^{k}(s^{*},\bm{x}^{*})\leq\mathcal{N}[v^{k}](\bm{x}^{*})$ , and then, $w^{k}(s^{*},\bm{x}^{*})\leq\beta$ . The lower bound $-\beta$ could be obtained similarly. Thus we obtain $\|w^{k}(s)\|\leq\beta$ for any $s\in[0,\tau]$ , and thus $\|v^{k+1}\|\leq\beta$ , which completes the proof. ∎

Next, we consider the higher-order approximations of the solution of (39). Let $P_{r}(s)$ be an interpolation of $\mathcal{N}[u(t_{n}+s)]$ with degree $r\geq 1$ on the times $\{s_{k}=\frac{k}{r}\tau\}_{k=0}^{r}$ . Then, we approximate $\mathcal{N}[u(t_{n}+s)]$ by $P_{r}(s)$ to obtain the higher-order ETD Runge–Kutta scheme: find $v^{n+1}=w^{n}(\tau)$ with $w^{n}:[0,\tau]\to\mathcal{X}$ solving

\begin{dcases}w^{n}_{t}+\kappa w^{n}={\mathcal{L}}w^{n}+P_{r}(s),&s\in(0,\tau],\ \bm{x}\in\Omega^{*},\\ w^{n}(s,\bm{x})=g(t_{n}+s,\bm{x}),&s\in[0,\tau],\ \bm{x}\in\Omega_{c}^{*},\\ w^{n}(0,\bm{x})=v^{n}(\bm{x}),&\bm{x}\in\widehat{\Omega},\end{dcases}

where $v^{n}$ is an approximation of $u(t_{n})$ with $v^{0}=u_{0}$ . More precisely, we have

P_{r}(s)=\sum_{k=0}^{r}\ell_{r,k}(s)\mathcal{N}[\tilde{v}^{n+\frac{k}{r}}],\quad s\in[0,\tau],

where $\{\ell_{r,k}(s)\}_{k=0}^{r}$ are the standard Lagrange basis functions associated with the times $\{s_{k}\}_{k=0}^{r}$ , and $\tilde{v}^{n+\frac{k}{r}}$ is an approximated value of $u(t_{n}+s_{k})$ , generated by the lower-order ETD schemes and $\tilde{v}^{n}=v^{n}$ . In the spirit of the proof of Theorem 7, we know that the discrete MBP would be preserved as long as the following condition is satisfied:

(42)

\|P_{r}(s)\|\leq\max\{\|\mathcal{N}[\tilde{v}^{n+\frac{k}{r}}]\|:0\leq k\leq r\},\quad\forall\,s\in[0,\tau],

with $\|\tilde{v}^{n+\frac{k}{r}}\|\leq\beta$ for all $k=0,1,\dots,r$ . However, the unique interpolation satisfying (42) corresponds to the case $r=1$ , that is, the linear interpolation

P_{1}(s)=\Big{(}1-\frac{s}{\tau}\Big{)}\mathcal{N}[\tilde{v}^{n}]+\frac{s}{\tau}\mathcal{N}[\tilde{v}^{n+1}],\quad s\in[0,\tau].

Approximating $\mathcal{N}[u(t_{n}+s)]$ in (39) by $P_{1}(s)$ , we obtain the second-order ETD Runge–Kutta (ETDRK2) scheme: find $v^{n+1}=w^{n}(\tau)$ with $w^{n}:[0,\tau]\to\mathcal{X}$ solving

(43)

\begin{dcases}w^{n}_{t}+\kappa w^{n}={\mathcal{L}}w^{n}+\Big{(}1-\frac{s}{\tau}\Big{)}\mathcal{N}[v^{n}]+\frac{s}{\tau}\mathcal{N}[\tilde{v}^{n+1}],&s\in(0,\tau],\ \bm{x}\in\Omega^{*},\\ w^{n}(s,\bm{x})=g(t_{n}+s,\bm{x}),&s\in[0,\tau],\ \bm{x}\in\Omega_{c}^{*},\\ w^{n}(0,\bm{x})=v^{n}(\bm{x}),&\bm{x}\in\widehat{\Omega},\end{dcases}

where $\tilde{v}^{n+1}$ is generated by the ETD1 scheme. It is worth noting that both ETD1 and ETDRK2 schemes are linear.

Theorem 8 (Maximum bound principle of the ETDRK2 scheme).

Suppose that Assumptions 1 and 2, (13) and (14) hold. Then the ETDRK2 scheme preserves the discrete MBP unconditionally, i.e., for any time step size $\tau>0$ , the ETDRK2 solution satisfies $\|v^{n}\|\leq\beta$ for any $n\geq 0$ .

Proof.

Again, we suppose $\|v^{k}\|\leq\beta$ for some $k$ . Similar to the proof of Theorem 7, for $(s^{*},\bm{x}^{*})\in(0,\tau]\times\Omega^{*}$ such that (41) holds, we can obtain

\kappa w^{k}(s^{*},\bm{x}^{*})\leq\Big{(}1-\frac{s^{*}}{\tau}\Big{)}\mathcal{N}[v^{k}](\bm{x}^{*})+\frac{s^{*}}{\tau}\mathcal{N}[\tilde{v}^{k+1}](\bm{x}^{*}),

and then, $w^{k}(s^{*},\bm{x}^{*})\leq\beta$ using $\|\tilde{v}^{k+1}\|\leq\beta$ by Theorem 7. The lower bound $-\beta$ could be obtained similarly. Therefore, $\|v^{k+1}\|\leq\beta$ , which completes the proof. ∎

As an application of Theorems 7 and 8, one can obtain the MBPs of the semi-discrete and fully discrete ETD schemes for the concrete examples presented in Sections 2.2 and 2.3.

Corollary 9.

For the space-continuous ETD1 and ETDRK2 schemes with ${\mathcal{L}}$ given by any of (21), (24), and (25), as well as the fully discrete ETD1 and ETDRK2 schemes with the (discretized) linear operator ${\mathcal{L}}={\mathcal{L}}_{h}$ defined by any of (27), (30), (31), or (32), it holds that

(i) if $f_{0}$ is given by (16) with even $p$ , $\kappa\geq\lambda p$ , and (14) is satisfied with $\beta=1$ , the solution to the ETD1 scheme (40) or the ETDRK2 scheme (43) satisfies $\|v^{n}\|\leq 1$ for any $n\geq 0$ ;

(ii) if $f_{0}$ is given by (18), $\kappa\geq\frac{\theta}{1-\rho^{2}}-\theta_{c}$ , and (14) is satisfied with $\beta=\rho$ , the solution to the ETD1 scheme (40) or the ETDRK2 scheme (43) satisfies $\|v^{n}\|\leq\rho$ for any $n\geq 0$ , where $\rho$ is the positive root of (19).

In summary, we have proved that the first- and second-order ETDRK schemes preserve the MBP of the underlying problem (5) unconditionally in the discrete sense. Meanwhile, we actually show that the classic ETDRK approximations, as defined in [11], with order greater than two fail to maintain the MBP unconditionally due to the lack of the property (42) for the interpolation polynomials with $r\geq 2$ . Apart from the ETDRK schemes, multistep approximations could be also utilized in the ETD temporal discretization [49]. The ETD-multistep scheme of order $r+1$ is generated from the extrapolation polynomial of degree $r\geq 1$ for the nonlinear term based on the approximated solutions at the previous $r$ time levels. Since the extrapolation polynomials with $r\geq 1$ cannot be bounded by the maximums and minimums of the extrapolated data, the ETD-multistep schemes with even order greater than one also fail to preserve the discrete MBP unconditionally.

Up till now, we only presented the differential forms of the ETD1 and ETDRK2 schemes by (40) and (43) since they are convenient for the theoretical analysis. In practice, we need formulas that can be more directly implemented for computations, in particular, the fully discrete schemes (in both time and space) with $X\simeq\mathbb{R}^{M^{d}}$ and ${\mathcal{L}}={\mathcal{L}}_{h}$ discussed in Section 2.3.2. Here we define a new operator ${\mathcal{L}}_{hc}:\partial X\to X$ as follows: for any $g\in\partial X$ and $\widehat{w}\in D({\mathcal{L}})$ with $\widehat{w}|_{\Omega_{c}^{*}}=g$ , let

(44)

{\mathcal{L}}_{hc}g=({\mathcal{L}}_{h}\widehat{w})|_{\Omega^{*}}-{\mathcal{L}}_{h0}(\widehat{w}|_{\Omega^{*}}).

It is easy to check that the right-hand side of (44) depends only on $g$ , so ${\mathcal{L}}_{hc}$ is well-defined through (44).

For a given boundary data $g:[0,T]\to\partial X$ and a function $u:[0,T]\to\mathcal{X}$ such that $u(t)|_{\Omega_{c}^{*}}=g(t)$ , we obtain from (39) and the definition of ${\mathcal{L}}_{hc}$ that

(45)

\begin{dcases}u_{t}+\kappa u={\mathcal{L}}_{h0}u+\mathcal{N}[u]+{\mathcal{L}}_{hc}g,&t\in(0,T],\ \bm{x}\in\Omega^{*},\\ u(0,\bm{x})=u_{0}(\bm{x}),&\bm{x}\in\widehat{\Omega},\end{dcases}

where ${\mathcal{L}}_{h0}$ is an $M^{d}$ -by- $M^{d}$ matrix. Let

(46)

{\mathcal{L}}_{\kappa,h}=\kappa I_{M^{d}}-{\mathcal{L}}_{h0}

and define the $\varphi$ -functions as follows:

\varphi_{0}(a)=\mathrm{e}^{-a},\quad\varphi_{1}(a)=\frac{1-\mathrm{e}^{-a}}{a},\quad\varphi_{2}(a)=\frac{a-1+\mathrm{e}^{-a}}{a^{2}},\qquad a\not=0.

Let $v^{0}=u_{0}$ . Then, we give the fully discrete ETD1 scheme by

(47)

v^{n+1}=\mathrm{e}^{-\tau{\mathcal{L}}_{\kappa,h}}v^{n}+\int_{0}^{\tau}\mathrm{e}^{-(\tau-s){\mathcal{L}}_{\kappa,h}}\big{\{}\mathcal{N}[v^{n}]+{\mathcal{L}}_{hc}g(t_{n}+s)\big{\}}\,\mathrm{d}s

or equivalently,

v^{n+1}=\varphi_{0}(\tau{\mathcal{L}}_{\kappa,h})v^{n}+\tau\varphi_{1}(\tau{\mathcal{L}}_{\kappa,h})\mathcal{N}[v^{n}]+\int_{0}^{\tau}\mathrm{e}^{-(\tau-s){\mathcal{L}}_{\kappa,h}}{\mathcal{L}}_{hc}g(t_{n}+s)\,\mathrm{d}s,

and the fully discrete ETDRK2 scheme by

\begin{dcases}\tilde{v}^{n+1}=\mathrm{e}^{-\tau{\mathcal{L}}_{\kappa,h}}v^{n}+\int_{0}^{\tau}\mathrm{e}^{-(\tau-s){\mathcal{L}}_{\kappa,h}}\big{\{}\mathcal{N}[v^{n}]+{\mathcal{L}}_{hc}g(t_{n}+s)\big{\}}\,\mathrm{d}s,\\ v^{n+1}=\mathrm{e}^{-\tau{\mathcal{L}}_{\kappa,h}}v^{n}+\int_{0}^{\tau}\mathrm{e}^{-(\tau-s){\mathcal{L}}_{\kappa,h}}\Big{\{}\Big{(}1-\frac{s}{\tau}\Big{)}\mathcal{N}[v^{n}]+\frac{s}{\tau}\mathcal{N}[\tilde{v}^{n+1}]+{\mathcal{L}}_{hc}g(t_{n}+s)\Big{\}}\,\mathrm{d}s\end{dcases}

or equivalently,

(48)

\begin{dcases}\tilde{v}^{n+1}=\varphi_{0}(\tau{\mathcal{L}}_{\kappa,h})v^{n}+\tau\varphi_{1}(\tau{\mathcal{L}}_{\kappa,h})\mathcal{N}[v^{n}]+\int_{0}^{\tau}\mathrm{e}^{-(\tau-s){\mathcal{L}}_{\kappa,h}}{\mathcal{L}}_{hc}g(t_{n}+s)\,\mathrm{d}s,\\ v^{n+1}=\tilde{v}^{n+1}+\tau\varphi_{2}(\tau{\mathcal{L}}_{\kappa,h})\big{\{}\mathcal{N}[\tilde{v}^{n+1}]-\mathcal{N}[v^{n}]\big{\}}.\end{dcases}

Note that the corresponding semigroup is given by the matrix exponential $S_{-{\mathcal{L}}_{\kappa,h}}(t)=\mathrm{e}^{t({\mathcal{L}}_{h0}-\kappa I)}$ , which again depends crucially on the choice of the constant $\kappa$ .

3.2 Convergence analysis of the ETD schemes

As an important application of the MBP, we now consider the convergence of the ETD1 and ETDRK2 schemes. From the practical point of view, we are mainly interested in the convergence of the fully discrete solution to the semi-discrete (in space) solution (with a fixed spatial mesh size) as the time step size goes to zero. Again, we have ${\mathcal{L}}={\mathcal{L}}_{h}$ and $S_{{\mathcal{L}}_{h0}}(t)=\mathrm{e}^{t{\mathcal{L}}_{h0}}$ as we discussed in Section 2.3.2. First, let us state the result for the ETD1 scheme (40).

Theorem 10.

Suppose that Assumptions 1 and 2, (13) and (14) are satisfied. For the fixed terminal time $T>0$ and spatial mesh size $h>0$ , assume that the exact solution $u$ to the space-discrete model equation (45) belongs to $C^{1}([0,T];\mathcal{X})$ and $\{v^{n}\}_{n\geq 0}$ is generated by the fully discrete ETD1 scheme (47) with $v^{0}=u_{0}$ . Then we have

(49)

\|v^{n}-u(t_{n})\|\leq C\mathrm{e}^{\kappa t_{n}}\tau,\quad\forall\,t_{n}\leq T

for any $\tau>0$ , where the constant $C>0$ depends on the $C^{1}([0,T];\mathcal{X})$ norm of $u$ , but independent of $\tau$ and $\kappa$ .

Proof.

We know that the ETD1 solution is given by $v^{n+1}=w^{n}(\tau)$ with the function $w^{n}:[0,\tau]\to\mathcal{X}$ solving

(50)

\begin{dcases}w^{n}_{s}+\kappa w^{n}={\mathcal{L}}_{h0}w^{n}+\mathcal{N}[v^{n}]+{\mathcal{L}}_{hc}g(t_{n}+s),&s\in(0,\tau],\ \bm{x}\in\Omega^{*},\\ w^{n}(0,\bm{x})=v^{n}(\bm{x}),&\bm{x}\in\widehat{\Omega}.\end{dcases}

Let $e_{1}^{n}=v^{n}-u(t_{n})$ . The difference between (50) and (45) yields

(51)

e_{1}^{n+1}=\mathrm{e}^{-\kappa\tau}\mathrm{e}^{\tau{\mathcal{L}}_{h0}}e_{1}^{n}+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\mathrm{e}^{(\tau-s){\mathcal{L}}_{h0}}\{\mathcal{N}[v^{n}]-\mathcal{N}[u(t_{n})]+R_{1}(s)\}\,\mathrm{d}s,

where $R_{1}(s)$ is the truncated error as

R_{1}(s)=\mathcal{N}[u(t_{n})]-\mathcal{N}[u(t_{n}+s)],\quad s\in[0,\tau].

Due to the MBP of (5), we have $\|u(t_{n})\|\leq\beta$ and $\|u(t_{n}+s)\|\leq\beta$ , and then, using Lemma 2-(ii), we derive

\|R_{1}(s)\|=\|N_{0}(u(t_{n}))-N_{0}(u(t_{n}+s))\|\leq 2\kappa\|u(t_{n})-u(t_{n}+s)\|\leq C_{1}\kappa\tau,\quad\forall\,s\in[0,\tau],

where the constant $C_{1}$ depends on the $C^{1}([0,T];\mathcal{X})$ norm of $u$ , but independent of $\tau$ and $\kappa$ . Similarly, since $\|v^{n}\|\leq\beta$ , we can obtain

(52)

\|\mathcal{N}[v^{n}]-\mathcal{N}[u(t_{n})]\|\leq 2\kappa\|v^{n}-u(t_{n})\|=2\kappa\|e_{1}^{n}\|.

Then, we derive from (51) that

	$\displaystyle\\|e_{1}^{n+1}\\|$	$\displaystyle\leq\mathrm{e}^{-\kappa\tau}\|\!\|\!\|\mathrm{e}^{\tau{\mathcal{L}}_{h0}}\|\!\|\!\|\\|e_{1}^{n}\\|$
		$\displaystyle\quad+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\|\!\|\!\|\mathrm{e}^{(\tau-s){\mathcal{L}}_{h0}}\|\!\|\!\|\big{\{}\\|\mathcal{N}[v^{n}]-\mathcal{N}[u(t_{n})]\\|+\\|R_{1}(s)\\|\big{\}}\,\mathrm{d}s$
		$\displaystyle\leq\mathrm{e}^{-\kappa\tau}\\|e_{1}^{n}\\|+(2\kappa\\|e_{1}^{n}\\|+C_{1}\kappa\tau)\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\,\mathrm{d}s$
		$\displaystyle=\mathrm{e}^{-\kappa\tau}\\|e_{1}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa}(2\kappa\\|e_{1}^{n}\\|+C_{1}\kappa\tau)$
		$\displaystyle=(2-\mathrm{e}^{\kappa\tau})\\|e_{1}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa\tau}\cdot C_{1}\kappa\tau^{2}$
(53)			$\displaystyle\leq(1+\kappa\tau)\\|e_{1}^{n}\\|+C_{1}\kappa\tau^{2},$

where in the last step we have used the fact that $1-\mathrm{e}^{-a}\leq a$ for any $a>0$ . By induction, we have

	$\displaystyle\\|e_{1}^{n}\\|$	$\displaystyle\leq(1+\kappa\tau)^{n}\\|e_{1}^{0}\\|+C_{1}\kappa\tau^{2}\sum_{k=0}^{n-1}(1+\kappa\tau)^{k}$
		$\displaystyle=(1+\kappa\tau)^{n}\\|e_{1}^{0}\\|+C_{1}\tau[(1+\kappa\tau)^{n}-1]$
		$\displaystyle\leq\mathrm{e}^{\kappa n\tau}\\|e_{1}^{0}\\|+C_{1}\mathrm{e}^{\kappa n\tau}\tau.$

By letting $C=C_{1}$ , we obtain (49) since $e_{1}^{0}=0$ and $n\tau=t_{n}$ . ∎

Remark 8.

We note that in the estimate (49), there is an exponential coefficient $\mathrm{e}^{\kappa t_{n}}$ which could be very large (e.g., in the case of Allen–Cahn equation, $f$ may include a negative power of a small parameter, which leads to large $\kappa$ ). The similar result is also obtained for the ETDRK2 scheme as shown later. Theoretically, this is inevitable due to the application of the Gronwall’s lemma. On the other hand, one may be able to further improve the error estimate by using the techniques in [34, 35].

Now, we turn to the convergence of the ETDRK2 scheme (43). Assume that $f_{0}$ is twice continuously differentiable and denote by $M^{f}_{2}$ the maximum of $|f_{0}^{\prime\prime}|$ on $[-\beta,\beta]$ .

Theorem 11.

Suppose that Assumptions 1 and 2, (13) and (14) are satisfied. For the fixed terminal time $T>0$ and spatial mesh size $h>0$ , assume that the exact solution $u$ to the space-discrete model equation (45) belongs to $C^{2}([0,T];\mathcal{X})$ and $\{v^{n}\}_{n\geq 0}$ is generated by the fully discrete ETDRK2 scheme (48) with $v^{0}=u_{0}$ . Then we have

(54)

\|v^{n}-u(t_{n})\|\leq C_{\kappa}\mathrm{e}^{\kappa t_{n}}\tau^{2},\quad\forall\,t_{n}\leq T

for any $\tau>0$ , where the constant $C_{\kappa}>0$ depends on $\kappa$ , $M^{f}_{2}$ , and the $C^{2}([0,T];\mathcal{X})$ norm of $u$ , but independent of $\tau$ .

Proof.

The proof strategy is quite similar to the case of the first-order scheme. Let $e_{2}^{n}=v^{n}-u(t_{n})$ , then we have

	$\displaystyle e_{2}^{n+1}$	$\displaystyle=\mathrm{e}^{-\kappa\tau}\mathrm{e}^{\tau{\mathcal{L}}_{h0}}e_{2}^{n}+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\mathrm{e}^{(\tau-s){\mathcal{L}}_{h0}}\Big{\{}\Big{(}1-\frac{s}{\tau}\Big{)}(\mathcal{N}[v^{n}]-\mathcal{N}[u(t_{n})])$
(55)			$\displaystyle\quad+\frac{s}{\tau}(\mathcal{N}[\tilde{v}^{n+1}]-\mathcal{N}[u(t_{n+1})])+R_{2}(s)\Big{\}}\,\mathrm{d}s,$

where $R_{2}(s)$ is the truncated error given by

R_{2}(s)=\Big{(}1-\frac{s}{\tau}\Big{)}\mathcal{N}[u(t_{n})]+\frac{s}{\tau}\mathcal{N}[u(t_{n+1})]-\mathcal{N}[u(t_{n}+s)],\quad s\in[0,\tau].

Using the MBP and the error estimates of the linear interpolation, we have

\|R_{2}(s)\|\leq(C_{2}\kappa+C_{3})\tau^{2},\quad\forall\,s\in[0,\tau],

where the constant $C_{2}$ depends on the $C^{2}([0,T];\mathcal{X})$ norm of $u$ and $C_{3}$ depends on $M^{f}_{2}$ and the $C^{2}([0,T];\mathcal{X})$ norm of $u$ ; both of them are independent of $\tau$ and $\kappa$ . From the last inequality in (53), we know

\|\tilde{v}^{n+1}-u(t_{n+1})\|\leq(1+\kappa\tau)\|v^{n}-u(t_{n})\|+C_{1}\kappa\tau^{2},

and then, using Lemma 2-(ii), we obtain

\|\mathcal{N}[\tilde{v}^{n+1}]-\mathcal{N}[u(t_{n+1})]\|\leq 2\kappa\|\tilde{v}^{n+1}-u(t_{n+1})\|\leq 2\kappa(1+\kappa\tau)\|e_{2}^{n}\|+2C_{1}\kappa^{2}\tau^{2}.

By combining with (52), we have, for any $s\in[0,\tau]$ ,

	$\displaystyle\Big{\\|}\Big{(}1-\frac{s}{\tau}\Big{)}(\mathcal{N}[v^{n}]-\mathcal{N}[u(t_{n})])+\frac{s}{\tau}(\mathcal{N}[\tilde{v}^{n+1}]-\mathcal{N}[u(t_{n+1})])\Big{\\|}$
	$\displaystyle\qquad\leq\Big{(}1-\frac{s}{\tau}\Big{)}2\kappa\\|e_{2}^{n}\\|+\frac{s}{\tau}[2\kappa(1+\kappa\tau)\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau^{2}]$
	$\displaystyle\qquad=2\kappa(1+\kappa s)\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau s.$

Then, we obtain from (55) that

	$\displaystyle\\|e_{2}^{n+1}\\|$	$\displaystyle\leq\mathrm{e}^{-\kappa\tau}\\|e_{2}^{n}\\|+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}[2\kappa(1+\kappa s)\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau s+(C_{2}\kappa+C_{3})\tau^{2}]\,\mathrm{d}s$
		$\displaystyle=\mathrm{e}^{-\kappa\tau}\\|e_{2}^{n}\\|+[2\kappa\\|e_{2}^{n}\\|+(C_{2}\kappa+C_{3})\tau^{2}]\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\,\mathrm{d}s$
		$\displaystyle\quad+(2\kappa^{2}\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau)\int_{0}^{\tau}s\mathrm{e}^{-\kappa(\tau-s)}\,\mathrm{d}s$
		$\displaystyle=\mathrm{e}^{-\kappa\tau}\\|e_{2}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa}[2\kappa\\|e_{2}^{n}\\|+(C_{2}\kappa+C_{3})\tau^{2}]$
		$\displaystyle\quad+\frac{\mathrm{e}^{-\kappa\tau}-1+\kappa\tau}{\kappa^{2}}(2\kappa^{2}\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau)$
		$\displaystyle=(\mathrm{e}^{-\kappa\tau}+2\kappa\tau)\\|e_{2}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa}(C_{2}\kappa+C_{3})\tau^{2}+\frac{\mathrm{e}^{-\kappa\tau}-1+\kappa\tau}{\kappa^{2}}\cdot 2C_{1}\kappa^{2}\tau$
		$\displaystyle\leq\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}\\|e_{2}^{n}\\|+(C_{2}\kappa+C_{3}+C_{1}\kappa^{2})\tau^{3},$

where we have used the inequality $1-a\leq\mathrm{e}^{-a}\leq 1-a+\frac{a^{2}}{2}$ for any $a>0$ . Finally, by induction, we obtain

	$\displaystyle\\|e_{2}^{n}\\|$	$\displaystyle\leq\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{n}\\|e_{2}^{0}\\|+(C_{1}\kappa^{2}+C_{2}\kappa+C_{3})\tau^{3}\sum_{k=0}^{n-1}\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{k}$
		$\displaystyle\leq\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{n}\\|e_{2}^{0}\\|+\Big{(}C_{1}\kappa+C_{2}+\frac{C_{3}}{\kappa}\Big{)}\tau^{2}\Big{[}\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{n}-1\Big{]}$
		$\displaystyle\leq\mathrm{e}^{\kappa n\tau}\\|e_{2}^{0}\\|+\Big{(}C_{1}\kappa+C_{2}+\frac{C_{3}}{\kappa}\Big{)}\mathrm{e}^{\kappa n\tau}\tau^{2}.$

By letting $C_{\kappa}=C_{1}\kappa+C_{2}+C_{3}/\kappa$ , we obtain (54). ∎

By specifying the linear operator by any one of Examples 2.7–2.10 and the nonlinear one by either Example 2.1, 2.2, or 2.3, one can obtain the convergence results of corresponding concrete problems.

3.3 Energy stability of the ETD schemes for gradient flow models

Next we show the application of the MBP and convergence of the ETD schemes to gradient flow models, a class of important examples of the model equation (5). We only consider the periodic or time-independent Dirichlet boundary condition and also assume that the linear operator ${\mathcal{L}}$ is dissipative on $L^{2}(\Omega)$ in the sense that the inner product $(w,{\mathcal{L}}w)_{L^{2}(\Omega)}\leq 0$ for any $w\in D({\mathcal{L}})$ , which can be satisfied by a large number of operators such as those given in Examples 2.4–2.10.

Phase-field (or diffuse-interface) models are usually derived as the gradient flows with respect to the energy functional

E[u]=-\frac{1}{2}(u,{\mathcal{L}}u)_{L^{2}(\Omega)}+\int_{\Omega}F(u(\bm{x}))\,\mathrm{d}\bm{x}

with $F:\mathbb{R}\to\mathbb{R}$ such that $f_{0}=-F^{\prime}$ . Some simple calculations give us the energy law under the periodic or time-independent Dirichlet boundary condition:

E[u(t_{2})]\leq E[u(t_{1})],\quad\forall\,t_{2}\geq t_{1}\geq 0.

There have been numerous researches devoted to the energy stable numerical methods for the phase-field models, see, e.g., [55, 65, 89, 98] and the references in the recent review [16].

An interesting problem is whether the energy law can be preserved by the proposed ETD discretization schemes. In our recent work [20], we investigated the MBP-preserving ETD schemes for the nonlocal Allen–Cahn equation with periodic boundary condition, namely, the equation (5) with the linear operator given by (24) and nonlinear function defined as (17). We concluded that the solution to the ETD1 scheme decreases the energy along the time steps and the energy for the ETDRK2 scheme is uniformly bounded by the initial energy plus a constant.

For the more general case we consider in this work, the ETD1 and ETDRK2 schemes for (5) still satisfy the energy stability. Since the proof could be conducted in a quite similar way as done in [20], we state the result directly as follows and leave the details of the proof to interested readers.

Proposition 12.

Suppose that Assumptions 1–2, (13) and (14) hold. Then

(i) the solution $\{v^{n}\}_{n\geq 0}$ to the ETD1 scheme (40) satisfies

E[v^{n+1}]\leq E[v^{n}],\quad\forall\,n\geq 0,

for any $\tau>0$ , i.e., the ETD1 scheme is unconditionally energy stable;

(ii) the solution $\{v^{n}\}_{n\geq 0}$ to the ETDRK2 scheme (43) satisfies

E[v^{n}]\leq E[v^{0}]+\widehat{C},\quad\forall\,t_{n}\leq T,

for any $\tau\in(0,1]$ , where the constant $\widehat{C}=\widehat{C}(|\Omega|,T,\kappa)$ is independent of $\tau$ , i.e., the energy is uniformly bounded.

4 Some extensions

In the framework we have established, although the function spaces $\mathcal{X}$ and $X$ consists of only real scalar-valued functions, the main results on the MBP could be extended to some complex scalar-valued equations, or more generally, real vector-valued ones. In addition, the equations on real matrix-valued fields also share some similar characteristics. The MBP for the case of systems implies the existence of the invariant regions of the solution [2, 63, 85, 91]. One can find studies on the invariant-region-preserving numerical methods for classic reaction-diffusion systems (see, e.g., [37, 72]) and hyperbolic systems of conservation laws (see, e.g., [43, 44, 54]). For simplicity, we focus our analysis on the space-continuous setting (i.e., Case (D1)) in this section.

4.1 Extension to complex scalar-valued and real vector-valued equations

The Ginzburg–Landau model for superconductivity [12, 19] is one of the popular models describing the phase transition of the superconducting material under effects of magnetic and electric fields. For simplicity, we only consider the equation with respect to the order parameter without electric effect, that is, the electric potential vanishing in the Ginzburg–Landau model. Let $\phi:[0,T]\times\overline{\Omega}\to\mathbb{C}$ be a complex-valued order parameter satisfying

(56)

\phi_{t}=(\nabla+\mathrm{i}\bm{A})^{2}\phi+(1-|\phi|^{2})\phi,\quad t\in(0,T],\ \bm{x}\in\Omega,

subject to the initial condition

\phi(0,\cdot)=\phi_{0},\quad\bm{x}\in\overline{\Omega}

and either the Dirichlet, periodic, or homogeneous natural boundary condition

(\nabla+\mathrm{i}\bm{A})\phi\cdot\bm{n}=0,\quad t\in(0,T],\ \bm{x}\in\partial\Omega,

where $\bm{A}\in\mathbb{R}^{d}$ is a given magnetic potential and $|\phi|$ denotes the modulus of $\phi$ . The MBP of the solution to (56) with the natural boundary condition was proved in [12] in the sense of weak solution as follows.

Proposition 13.

If $|\phi_{0}(\bm{x})|\leq 1$ for a.e. $\bm{x}\in\overline{\Omega}$ , then it holds $|\phi(t,\bm{x})|\leq 1$ for a.e. $t\in[0,T]$ and $\bm{x}\in\overline{\Omega}$ .

If we let $\psi=\mathrm{e}^{\mathrm{i}\bm{A}\cdot\bm{x}}\phi$ , simple calculations give us $\Delta\psi=\mathrm{e}^{\mathrm{i}\bm{A}\cdot\bm{x}}(\nabla+\mathrm{i}\bm{A})^{2}\phi$ and $|\psi|=|\phi|$ . Thus, the equation (56) is equivalent to

(57)

\psi_{t}=\Delta\psi+(1-|\psi|^{2})\psi,\quad t\in(0,T],\ \bm{x}\in\Omega,

subject to the initial condition $\psi(0,\bm{x})=\mathrm{e}^{\mathrm{i}\bm{A}\cdot\bm{x}}\phi_{0}(\bm{x})$ and corresponding boundary conditions. Noting that $\psi$ and $\phi$ have the same modulus, the MBP (Proposition 13) is also valid for (57). Since a complex number can be viewed as an element in $\mathbb{R}^{2}$ in the sense of isomorphism, the complex-valued equation (57) is actually equivalent to a real vector-valued equation

(58)

\bm{u}_{t}=\Delta\bm{u}+(1-|\bm{u}|^{2})\bm{u},\quad t\in(0,T],\ \bm{x}\in\Omega

with $|\cdot|$ denoting the standard Euclidean norm, where $\bm{u}:[0,T]\times\overline{\Omega}\to\mathbb{R}^{m}$ ( $m=2$ for the Ginzburg–Landau model) is subject to the initial condition

\bm{u}(0,\cdot)=\bm{u}_{0},\quad\bm{x}\in\overline{\Omega}

and either the Dirichlet boundary condition

\bm{u}=\bm{g},\quad t\in[0,T],\ \bm{x}\in\partial\Omega

with $\bm{g}\in C([0,T]\times\partial\Omega;\mathbb{R}^{m})$ , the periodic boundary condition, or the homogeneous Neumann boundary condition

\nabla\bm{u}\cdot\bm{n}=\bm{0},\quad t\in[0,T],\ \bm{x}\in\partial\Omega.

Similar to the Allen–Cahn equation, the vector-valued equation (58) could be also regarded as the $L^{2}$ gradient flow with respect to the energy functional

(59)

E[\bm{u}]=\int_{\Omega}\Big{(}\frac{1}{2}|\nabla\bm{u}|_{F}^{2}+\frac{1}{4}(|\bm{u}|^{2}-1)^{2}\Big{)}\,\mathrm{d}\bm{x},

where $|\nabla\bm{u}|_{F}$ denotes the Frobenius norm of the Jacobian matrix $\nabla\bm{u}$ . The solution to (58) decreases the energy (59) along with the time under either the time-independent Dirichlet, the homogeneous Neumann, or the periodic boundary condition.

Introducing the stabilizing constant $\kappa>0$ as before, the equation (58) is then equivalent to

\bm{u}_{t}+\kappa\bm{u}=\Delta\bm{u}+\bm{N}_{0}(\bm{u}),\quad t\in(0,T],\ \bm{x}\in\Omega,

where $\bm{N}_{0}(\bm{\xi})=\kappa\bm{\xi}+\bm{f}_{0}(\bm{\xi})$ and $\bm{f}_{0}(\bm{\xi})=(1-|\bm{\xi}|^{2})\bm{\xi}$ is the vector-valued analogue of the scalar function (17). Corresponding to (13) and Lemma 2, we choose $\kappa\geq 2$ and then have the following lemma on the vector-valued function $\bm{N}_{0}$ .

Lemma 14.

Denote $Y_{1}=\{\bm{\xi}\in\mathbb{R}^{m}\,|\,|\bm{\xi}|\leq 1\}$ . It holds that

(i) $|\bm{N}_{0}(\bm{\xi})|\leq\kappa$ for any $\bm{\xi}\in Y_{1}$ ;

(ii) $|\bm{N}_{0}(\bm{\xi}_{1})-\bm{N}_{0}(\bm{\xi}_{2})|\leq 2\kappa|\bm{\xi}_{1}-\bm{\xi}_{2}|$ for any $\bm{\xi}_{1},\bm{\xi}_{2}\in Y_{1}$ .

Proof.

(i) For $\bm{\xi}\in Y_{1}$ , we note that

|\bm{N}_{0}(\bm{\xi})|\leq\kappa|\bm{\xi}|+|\bm{f}_{0}(\bm{\xi})|=\kappa|\bm{\xi}|+f_{0}(|\bm{\xi}|)=N_{0}(|\bm{\xi}|),

where $N_{0}$ is the scalar function (12). The rest follows the proof of Lemma 2-(i) with $\beta=1$ .

(ii) The Jacobian matrix of $\bm{N}_{0}$ at $\bm{\xi}\in\mathbb{R}^{m}$ is given by

\nabla\bm{N}_{0}(\bm{\xi})=(\kappa+1-|\bm{\xi}|^{2})I_{m}-2(\bm{\xi}\otimes\bm{\xi}),

where $\otimes$ denotes the tensor product. The $m$ eigenvalues of $\nabla\bm{N}_{0}(\bm{\xi})$ are given by

\lambda_{1}(\bm{\xi})=\kappa+1-3|\bm{\xi}|^{2},\quad\lambda_{2}(\bm{\xi})=\lambda_{3}(\bm{\xi})=\cdots=\lambda_{m}(\bm{\xi})=\kappa+1-|\bm{\xi}|^{2}.

Since $\kappa\geq 2$ , for any $\bm{\xi}\in Y_{1}$ we have

0\leq\lambda_{1}(\bm{\xi})\leq\lambda_{m}(\bm{\xi})\leq\kappa+1<2\kappa.

Thus, using the mean-value theorem, we obtain (ii). ∎

Let $\mathcal{Y}=C(\overline{\Omega};\mathbb{R}^{m})$ the space of continuous $\mathbb{R}^{m}$ -valued functions defined on $\overline{\Omega}$ equipped with the supremum norm

\|\bm{w}\|_{\mathcal{Y}}=\max_{\bm{x}\in\overline{\Omega}}|\bm{w}(\bm{x})|,\quad\forall\,\bm{w}\in{\mathcal{Y}}.

By using Banach’s fixed-point theorem as done in the proof of Theorem 3, it is easy to show the existence and MBP of the problem (58) with any integer $m\geq 2$ as follows. The main tools are the uniform ellipticity of the Laplace operator (see, e.g., [28]), the properties of the nonlinear term given by Lemma 14, and the inequality

(60)

\Delta|\bm{u}|^{2}=2|\nabla\bm{u}|_{F}^{2}+2\bm{u}\cdot\Delta\bm{u}\geq 2\bm{u}\cdot\Delta\bm{u}.

Proposition 15.

If it holds that

(61a)		$\|\bm{u}_{0}(\bm{x})\|\leq 1,\quad\forall\,\bm{x}\in\overline{\Omega},$
the equation (58) with either periodic, homogeneous Neumann, or Dirichlet boundary condition subject to
(61b)		$\|\bm{g}(t,\bm{x})\|\leq 1,\quad\forall\,t\in[0,T],\ \forall\,\bm{x}\in\partial\Omega$

has a unique solution $\bm{u}\in C([0,T];\mathcal{Y})$ and it satisfies $\|\bm{u}(t)\|_{\mathcal{Y}}\leq 1$ for any $t\in[0,T]$ .

This proposition actually implies that the closed unit ball in $\mathbb{R}^{m}$ is an invariant region of the solution to the equation (58). Moreover, according to Corollary 14.8-(b) in [91], the closed unit ball is the smallest invariant region.

Next we show that, the ETD1 and ETDRK2 schemes for time integration of the equation (58) both preserve the discrete MBP unconditionally. Let $\bm{v}^{0}=\bm{u}_{0}$ . For the case of Dirichlet boundary condition, the ETD1 and ETDRK2 solutions are then given by $\bm{v}^{n+1}=\bm{w}^{n}(\tau)$ with $\bm{w}^{n}:[0,\tau]\to\mathcal{Y}$ such that

(62)

\begin{dcases}\bm{w}^{n}_{s}+\kappa\bm{w}^{n}=\Delta\bm{w}^{n}+\widehat{\bm{N}}_{0}(s,\bm{v}^{n}),&s\in(0,\tau],\ \bm{x}\in\Omega,\\ \bm{w}^{n}(s,\bm{x})=\bm{g}(t_{n}+s,\bm{x}),&s\in[0,\tau],\ \bm{x}\in\partial\Omega,\\ \bm{w}^{n}(0,\bm{x})=\bm{v}^{n}(\bm{x}),&\bm{x}\in\overline{\Omega},\end{dcases}

where

\widehat{\bm{N}}_{0}(s,\bm{v}^{n})=\begin{dcases}\bm{N}_{0}(\bm{v}^{n}),&\text{for ETD1},\\ \Big{(}1-\frac{s}{\tau}\Big{)}\bm{N}_{0}(\bm{v}^{n})+\frac{s}{\tau}\bm{N}_{0}(\tilde{\bm{v}}^{n+1}),&\text{for ETDRK2}\end{dcases}

with $\tilde{\bm{v}}^{n+1}$ generated by the ETD1 scheme. For the cases of periodic or homogeneous Neumann boundary condition, the ETD1 and ETDRK2 solutions are still given by the system (62) with small modifications of removing its second equation and using corresponding properties of the solution on the boundary.

Theorem 16.

Assume that the stabilizing constant $\kappa\geq 2$ and (61) holds. Then the ETD1 and ETDRK2 schemes of the equation (58) both preserve the discrete MBP unconditionally, i.e., for any time step size $\tau>0$ , the solutions satisfy $\|\bm{v}^{n}\|_{\mathcal{Y}}\leq 1$ for any $n\geq 0$ .

Proof.

Since $\|\bm{v}^{0}\|_{\mathcal{Y}}\leq 1$ , we just need to show that $\|\bm{v}^{k}\|_{\mathcal{Y}}\leq 1$ implies $\|\bm{v}^{k+1}\|_{\mathcal{Y}}\leq 1$ for any $k$ . We have $\bm{v}^{k+1}=\bm{w}^{k}(\tau)$ , where $\bm{w}^{k}$ satisfies (62) with the superscript $n$ replaced by $k$ . Taking the dot product of the first equation in (62) with $\bm{w}^{k}$ and using the fact from (60) that $\Delta|\bm{w}^{k}|^{2}\geq 2\bm{w}^{k}\cdot\Delta\bm{w}^{k}$ , we obtain

(63)

\frac{1}{2}(|\bm{w}^{k}|^{2})_{s}+\kappa|\bm{w}^{k}|^{2}\leq\frac{1}{2}\Delta|\bm{w}^{k}|^{2}+|\widehat{\bm{N}}_{0}(s,\bm{v}^{k})|\,|\bm{w}^{k}|.

Suppose there exists $(s^{*},\bm{x}^{*})\in(0,\tau]\times\widehat{\Omega}$ such that

|\bm{w}^{k}(s^{*},\bm{x}^{*})|=\max_{0\leq s\leq\tau}\|\bm{w}^{k}(s)\|_{\mathcal{Y}}.

Since $|\bm{w}^{k}(s,\bm{x})|^{2}$ is a real scalar-valued function, we have $(|\bm{w}^{k}|^{2})_{s}\geq 0$ at $(s^{*},\bm{x}^{*})$ . If $\bm{x}^{*}\in\Omega$ , we have $\Delta|\bm{w}^{k}(s^{*},\bm{x}^{*})|^{2}\leq 0$ . If $\bm{x}^{*}\in\partial\Omega$ , we have $|\bm{w}^{k}(s^{*},\bm{x}^{*})|\leq 1$ for the case of Dirichlet boundary condition (the proof is then in fact completed in this case) or $\Delta|\bm{w}^{k}(s^{*},\bm{x}^{*})|^{2}\leq 0$ for the cases of periodic and homogeneous Neumann boundary conditions. Then we obtain from (63) that $\kappa|\bm{w}^{k}(s^{*},\bm{x}^{*})|\leq|\widehat{\bm{N}}_{0}(s^{*},\bm{v}^{k}(\bm{x}^{*}))|$ . Since $\|\bm{v}^{k}\|_{\mathcal{Y}}\leq 1$ , according to Lemma 14-(i), for both ETD1 and ETDRK2 schemes, we always have $|\widehat{\bm{N}}_{0}(s^{*},\bm{v}^{k}(\bm{x}^{*}))|\leq\kappa$ , and thus $|\bm{w}^{k}(s^{*},\bm{x}^{*})|\leq 1$ . Then we have $\|\bm{v}^{k+1}\|_{\mathcal{Y}}\leq 1$ , which completes the proof. ∎

Remark 9.

For the space-discrete version of the equation (58), an essential condition for the MBP to hold is that $\Delta_{h}$ , the spatial discretization of the operator $\Delta$ , satisfies $\Delta_{h}|\bm{u}|^{2}\geq 2\bm{u}\cdot\Delta_{h}\bm{u}$ .

Similar to the scalar-valued problem, here we present the fully discrete ETD1 and ETDRK2 schemes for practical computations. With ${\mathcal{L}}=\Delta$ , we still use the notations ${\mathcal{L}}_{hc}$ and ${\mathcal{L}}_{\kappa,h}$ as defined by (44) and (46) respectively, then the fully discrete ETDRK2 scheme reads

\begin{dcases}\tilde{\bm{v}}^{n+1}=\varphi_{0}(\tau{\mathcal{L}}_{\kappa,h})\bm{v}^{n}+\tau\varphi_{1}(\tau{\mathcal{L}}_{\kappa,h})\bm{N}_{0}(\bm{v}^{n})+\int_{0}^{\tau}\mathrm{e}^{-(\tau-s){\mathcal{L}}_{\kappa,h}}{\mathcal{L}}_{hc}\bm{g}(t_{n}+s)\,\mathrm{d}s,\\ \bm{v}^{n+1}=\tilde{\bm{v}}^{n+1}+\tau\varphi_{2}(\tau{\mathcal{L}}_{\kappa,h})[\bm{N}_{0}(\tilde{\bm{v}}^{n+1})-\bm{N}_{0}(\bm{v}^{n})],\end{dcases}

and the fully discrete ETD1 scheme is given by the first step of the ETDRK2 scheme. The schemes for cases of homogeneous Neumann and periodic boundary conditions could be given in the similar way.

4.2 Extension to real matrix-valued equations

In [79], a problem of finding the stationary points of an energy for orthogonal matrix-valued functions was studied. Since the orthogonality constraint is non-trivial to enforce, a penalty term is added to the energy to offer a relaxed (phase-field or diffuse-interface) formulation. The gradient flow for such energy reads

(64)

U_{t}=\Delta U+U(I_{m}-U^{T}U),\quad t\in(0,T],\ \bm{x}\in\Omega,

where $U:[0,T]\times\overline{\Omega}\to\mathbb{R}^{m\times m}$ is subject to the initial condition

U(0,\cdot)=U_{0},\quad\bm{x}\in\overline{\Omega}

and either homogeneous Dirichlet, periodic, or homogeneous Neumann boundary condition. Denote by $|\cdot|_{2}$ the matrix $2$ -norm and by $\mathbb{R}^{m\times m}_{s}$ the set of all real symmetric $m$ -by- $m$ matrices.

Let us define $\mathcal{Z}=C(\overline{\Omega};\mathbb{R}^{m\times m}_{s})$ the space of continuous $\mathbb{R}^{m\times m}_{s}$ -valued functions defined on $\overline{\Omega}$ equipped with the supremum norm

\|W\|_{\mathcal{Z}}=\max_{\bm{x}\in\overline{\Omega}}|W(\bm{x})|_{2},\quad\forall\,W\in\mathcal{Z}.

Similarly to $X$ for the scalar-value equation case, we then define $Z$ , as well as $D({\mathcal{L}}_{0})$ and ${\mathcal{L}}_{0}$ , in accordance with ${\mathcal{L}}=\Delta$ and the respective boundary conditions. Then, we can show that ${\mathcal{L}}_{0}:D({\mathcal{L}}_{0})\to Z$ satisfies the matrix-valued analogue of Lemma 1.

Lemma 17.

For all $\lambda>0$ and all $W\in D({\mathcal{L}}_{0})$ , it holds

(65)

\|(\lambda\mathcal{I}-\Delta)W\|_{\mathcal{Z}}\geq\lambda\|W\|_{\mathcal{Z}},

and thus ${\mathcal{L}}_{0}$ generates a contraction semigroup $\{S_{{\mathcal{L}}_{0}}(t)\}_{t\geq 0}$ on $Z$ , i.e., $|\!|\!|S_{{\mathcal{L}}_{0}}(t)|\!|\!|\leq 1$ .

Proof.

First, for any diagonal matrix $W(\bm{x})=\mathop{\operator@font diag}\nolimits\{w_{i}(\bm{x}):1\leq i\leq m\}$ , there exists $\bm{x}_{0}\in\Omega$ (for the homogeneous Dirchlet boundary condition) or $\bm{x}_{0}\in\overline{\Omega}$ (for the periodic or homogeneous Neumann boundary condition) and $i_{0}$ such that

\|W\|_{\mathcal{Z}}=|W(\bm{x}_{0})|_{2}=|w_{i_{0}}(\bm{x}_{0})|=\max_{\bm{x}\in\overline{\Omega}}|w_{i_{0}}(\bm{x})|.

Since $|w_{i_{0}}(\bm{x})|^{2}$ is a real scalar-valued function, we have

2w_{i_{0}}(\bm{x}_{0})\Delta w_{i_{0}}(\bm{x}_{0})\leq 2w_{i_{0}}(\bm{x}_{0})\Delta w_{i_{0}}(\bm{x}_{0})+2|\nabla w_{i_{0}}(\bm{x}_{0})|^{2}=\Delta|w_{i_{0}}(\bm{x}_{0})|^{2}\leq 0.

Then, for any $\lambda>0$ , we have

	$\displaystyle\lambda\|W(\bm{x}_{0})\|_{2}^{2}$	$\displaystyle\leq\lambda\|w_{i_{0}}(\bm{x}_{0})\|^{2}-w_{i_{0}}(\bm{x}_{0})\Delta w_{i_{0}}(\bm{x}_{0})$
(66)			$\displaystyle=w_{i_{0}}(\bm{x}_{0})\cdot(\lambda\mathcal{I}-\Delta)w_{i_{0}}(\bm{x}_{0})\leq\|W(\bm{x}_{0})\|_{2}\|(\lambda\mathcal{I}-\Delta)W(\bm{x}_{0})\|_{2},$

which implies (65) for any diagonal matrix $W$ .

Next, for any $W\in D({\mathcal{L}}_{0})$ , let $\bm{x}_{0}$ be the point such that

\|W\|_{\mathcal{Z}}=|W(\bm{x}_{0})|_{2}=\max_{\bm{x}\in\overline{\Omega}}|W(\bm{x})|_{2}.

Since $W(\bm{x}_{0})$ is symmetric, there exists an orthonormal matrix $O$ such that $\widehat{W}=O^{T}W(\bm{x}_{0})O$ is diagonal. We then derive from (66) that

\lambda|W(\bm{x}_{0})|_{2}=\lambda|\widehat{W}|_{2}\leq|(\lambda\mathcal{I}-\Delta)\widehat{W}|_{2}=|(\lambda\mathcal{I}-\Delta)W(\bm{x}_{0})|_{2}\leq\|(\lambda\mathcal{I}-\Delta)W\|_{\mathcal{Z}},

which leads to (65). ∎

Introducing a stabilizing constant $\kappa>0$ , the equation (64) is equivalent to

(67)

U(t+\tau)=\mathrm{e}^{-\kappa\tau}S_{{\mathcal{L}}_{0}}(\tau)U(t)+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}S_{{\mathcal{L}}_{0}}(\tau-s)\mathcal{N}_{0}(U(t+s))\,\mathrm{d}s,

where

(68)

\mathcal{N}_{0}(Q)=\kappa Q+Q(I_{m}-Q^{T}Q),\quad Q\in\mathbb{R}^{m\times m}.

We again require $\kappa\geq 2$ , and then obtain the following lemma on the matrix-valued function $\mathcal{N}_{0}$ .

Lemma 18.

Denote ${\mathcal{M}}_{1}=\{Q\in\mathbb{R}^{m\times m}_{s}\,|\,|Q|_{2}\leq 1\}$ . It holds that

(i) $|\mathcal{N}_{0}(Q)|_{2}\leq\kappa$ for any $Q\in{\mathcal{M}}_{1}$ ;

(ii) $|\mathcal{N}_{0}(Q_{1})-\mathcal{N}_{0}(Q_{2})|_{2}\leq 2\kappa|Q_{1}-Q_{2}|_{2}$ for any $Q_{1},Q_{2}\in{\mathcal{M}}_{1}$ .

Proof.

(i) Since any real symmetric matrix can be diagonalized orthonormally, and $\mathcal{N}_{0}(Q)$ is also diagonal for any diagonal matrix $Q\in{\mathcal{M}}_{1}$ , the property (i) is the direct consequence of Lemma 2-(i).

(ii) Since $\mathbb{R}^{m\times m}$ is identical to $\mathbb{R}^{m^{2}}$ in the sense of isomorphism, the matrix-valued function $\mathcal{N}_{0}:\mathbb{R}^{m\times m}\to\mathbb{R}^{m\times m}$ defined by (68) could be regarded as a vector-valued mapping $\mathbb{R}^{m^{2}}\to\mathbb{R}^{m^{2}}$ , whose Jacobian matrix gives the matrix derivative of $\mathcal{N}_{0}$ . In this sense, the matrix derivative of $\mathcal{N}_{0}$ at $Q\in\mathbb{R}^{m\times m}_{s}$ is given by [71, Theorem 4]

\mathrm{D}\mathcal{N}_{0}(Q)=(\kappa+1)I_{m^{2}}-(Q^{2}\otimes I_{m}+I_{m}\otimes Q^{2}+Q\otimes Q).

Denote by $\{\mu_{j}\}_{j=1}^{m}$ the eigenvalues of $Q$ . Then the eigenvalues of $\mathrm{D}\mathcal{N}_{0}(Q)$ , denoted by $\{\lambda_{ij}\}_{i,j=1}^{m}$ , are given by [50, Theorem 4.2.12]

\lambda_{ij}=\kappa+1-(\mu_{i}^{2}+\mu_{j}^{2}+\mu_{i}\mu_{j}),\quad 1\leq i,j\leq m.

For any $Q\in{\mathcal{M}}_{1}$ , it holds $0\leq\mu_{i}^{2}+\mu_{j}^{2}+\mu_{i}\mu_{j}\leq 3$ . Since $\kappa\geq 2$ , we have

0\leq\lambda_{ij}\leq\kappa+1<2\kappa,\quad 1\leq i,j\leq m.

Thus, we obtain the property (ii) by using the mean-value theorem. ∎

By conducting the similar analysis as done in Section 2.1 and [20], we can prove the MBP for the matrix-valued equation (64).

Proposition 19.

If $U_{0}(\bm{x})$ is symmetric and $|U_{0}(\bm{x})|_{2}\leq 1$ for any $\bm{x}\in\overline{\Omega}$ , then the equation (64) with either homogeneous Dirichlet, periodic, or homogeneous Neumann boundary condition has a unique solution $U\in C([0,T];\mathcal{Z})$ and it satisfies $\|U(t)\|_{\mathcal{Z}}\leq 1$ for any $t\in[0,T]$ .

Proof.

Setting $t=0$ in (67) gives us

U(\tau)=\mathrm{e}^{-\kappa\tau}S_{{\mathcal{L}}_{0}}(\tau)U_{0}+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}S_{{\mathcal{L}}_{0}}(\tau-s)\mathcal{N}_{0}(U(s))\,\mathrm{d}s,\quad\tau\geq 0.

Denote $Z_{1}=\{W\in Z\,|\,\|W\|_{\mathcal{Z}}\leq 1\}$ . For a fixed $t_{1}>0$ and a given $V\in C([0,t_{1}];Z_{1})$ , let us define $W:[0,t_{1}]\to Z$ by

(69)

W(\tau)=\mathrm{e}^{-\kappa\tau}S_{{\mathcal{L}}_{0}}(\tau)U_{0}+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}S_{{\mathcal{L}}_{0}}(\tau-s)\mathcal{N}_{0}(V(s))\,\mathrm{d}s,\quad\tau\in[0,t_{1}].

Obviously, $W$ is uniquely defined and

\|W(\tau)\|_{\mathcal{Z}}\leq\mathrm{e}^{-\kappa\tau}|\!|\!|S_{{\mathcal{L}}_{0}}(\tau)|\!|\!|\|U_{0}\|_{\mathcal{Z}}+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}|\!|\!|S_{{\mathcal{L}}_{0}}(\tau-s)|\!|\!|\|\mathcal{N}_{0}(V(s))\|_{\mathcal{Z}}\,\mathrm{d}s.

According to Lemmas 17 and 18-(i), we derive

\|W(\tau)\|_{\mathcal{Z}}\leq\mathrm{e}^{-\kappa\tau}+\Big{(}\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\,\mathrm{d}s\Big{)}\kappa=\mathrm{e}^{-\kappa\tau}+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa}\cdot\kappa=1,\quad\tau\in[0,t_{1}],

which means $W\in C([0,t_{1}];Z_{1})$ .

Then, for $t_{1}<\kappa^{-1}\ln 2$ , similar to the proof of Theorem 3 and using Lemma 18-(ii), the mapping $\mathcal{A}:C([0,t_{1}];Z_{1})\to C([0,t_{1}];Z_{1})$ defined by $\mathcal{A}[V]=W$ according to (69) is a contraction. We then can conclude, similarly to the earlier scalar case, a unique solution $U(t)\in Z_{1}$ to the matrix-valued equation (64) exists on the time interval $[0,t_{1}]$ and can be further extended to $[0,T]$ . This completes the proof. ∎

By using Lemma 18, it also can be shown that both the ETD1 and ETDRK2 schemes for time integration of (64) again preserve the discrete MBP unconditionally when the stabilizing constant $\kappa\geq 2$ .

5 Numerical experiments

There exists a large amount of literature comparing and showing the excellent performance of the ETD schemes in numerical simulations for local continuum and nonlocal models [20, 55, 57, 58, 99, 105, 107]. The practical efficiency of the fully discrete ETD schemes (i.e., the finite dimensional cases) depends highly on the implementation of the actions of the operator/matrix exponentials. We first review existing algorithms for computing the products of the matrix exponentials with vectors, and then we present some detailed experimental results.

5.1 Implementations of matrix exponentials

Let us recall the fully discrete ETDRK2 scheme (48). There exist many efficient algorithms for computing the matrix functions $\varphi_{\gamma}(\tau{\mathcal{L}}_{\kappa,h})$ , $\gamma=0,1,2$ and their products with vectors.

When the spatial domain of the problem is regular, for instance, a rectangle $\prod_{j=1}^{d}(a_{j},b_{j})$ , and the matrix ${\mathcal{L}}_{\kappa,h}$ has some certain special structure, fast Fourier transform (FFT) based algorithms are adequate to calculate the above products of $\varphi$ -functions with vectors. This case arises commonly in the models for material sciences. For example, the matrices ${\mathcal{L}}_{h}$ given in Examples 2.7 and 2.8 are symmetric Toeplitz matrices for Case (C1) or circulant matrices for Case (C2). As we know, the product of a circulant matrix with a vector actually gives a circulant convolution, which could be implemented by using the FFT. For Toeplitz matrices, their products with a vector can be calculated by the sine or cosine transform. Alternatively, one can expand a Toeplitz matrix to a large circulant matrix and again make use of FFT for fast implementations.

When the shape of the spatial domain is arbitrary or the matrix ${\mathcal{L}}_{\kappa,h}$ does not possess a certain special structure, it is usually difficult to develop fast algorithms for matrix exponentials and their products with vector. In [45, 74, 75], many methods are surveyed for computing the exponential of a matrix, such as Taylor series, ODE solver, inverse Laplace transform, matrix decomposition and so on. The Matlab built-in command expm(A) computes the matrix exponential of $A$ by using a scaling and squaring algorithm with a Padé approximation. The performance of these methods often depends on the target problems and some of them are appropriate only for certain special problems. In the past two decades, Krylov subspace method based on Arnoldi or Lanczos iterations has become a powerful tool for computing the products of matrix exponentials with vectors [38, 46, 47, 90], especially for large-scale sparse problems. More recently, this method also has been further combined with incomplete orthogonalization [39] and adaptive time stepping [40, 77].

Remark 10.

In the ETD1 scheme, by approximating $\mathrm{e}^{\tau{\mathcal{L}}_{\kappa,h}}\approx\mathcal{I}+\tau{\mathcal{L}}_{\kappa,h}$ , one can obtain the first-order semi-implicit scheme, which is linear just as the ETD1 scheme and also preserves the MBP [88, 93]. In general, the ETD1 scheme is slightly more time-consuming than the classic semi-implicit schemes but with the added benefits of being more accurate [55, 57]. The fully-implicit backward Euler scheme naturally possesses good numerical stability (no extra stabilization term is needed), but it is not linear and usually more time-consuming due to the need of nonlinear iterations at each time step. A nonlinear Crank–Nicolson scheme was proposed and proved to conditionally preserve the MBP in [51]. To our best knowledge, there is so far no other linear second-order (in time) scheme like the ETDRK2 scheme which is unconditionally MBP-preserving. Furthermore, the ETD schemes preserve the exponential behavior of the linear operator, which is often crucial in the practical simulations of stiff systems (i.e., the spectral radius $\rho({\mathcal{L}}_{\kappa,h})\gg 1$ ).

5.2 Experimental results

Two examples will be tested to show the MBP-preserving property and numerical performance of the ETD methods. The first example focuses on the scalar equation (5) with the nonlinear function (18) consisting of logarithmic terms. The second example solves the vector-valued equation (58) to show some interesting evolutions of vortices in composite domains. In all experiments, the ETDRK2 scheme with uniform time step size $\tau=0.01$ is used for time integration.

Example 5.1.

Consider the scalar equation (5) of the unknown function $u:\Omega\subset\mathbb{R}^{2}\to\mathbb{R}$ in the domain $\Omega=(0,2\pi)\times(0,2\pi)$ with ${\mathcal{L}}=0.01\Delta$ and $f_{0}$ given by (18), i.e., the negative of the derivative of Flory–Huggins potential. The maximum bound principle plays an important role in this case since the equation consists of the logarithmic terms which will yield complex numbers if the value of the solution is located out of the interval $(-1,1)$ . The periodic and homogeneous Neumann boundary conditions are considered.

We use the uniform rectangular mesh with size $h=2\pi/512$ for partition of the domain in this example and a random data ranging from $-0.9$ to $0.9$ is generated on the mesh as the initial configuration of $u$ . The spatial discretization is done through the central finite difference as Example 2.7 so that the FFT can be used here for fast calculations needed in the ETDRK2 scheme. We set $\theta=0.8$ and $\theta_{c}=1.6$ in (18). According to Example 2.2, the positive root of $f_{0}(\rho)=0$ is $\rho\approx 0.9575$ and the stabilizing constant is thus chosen as $\kappa=8.02$ .

Fig. 1 shows the configurations of the approximate solution $u$ at $t=1$ , $5$ , $8$ , $30$ , $100$ , and $160$ subject to the periodic boundary condition. The simulation results illustrate the dynamics beginning with a random state and towards the homogeneous steady state of $u\equiv-\rho$ , which is reached after about $t=166$ in our simulations. The evolution of the energy is plotted in Fig. 3-(left) and that of the supremum norm of $u$ in Fig. 4-(left). The simulation results subject to the homogeneous Neumann boundary condition are presented in Fig. 2, where the same homogeneous steady state is reached after about $t=547$ . Fig. 3-(right) and Fig. 4-(right) show the corresponding evolutions of the energy and the supremum norm respectively. We observe that the energy decreases monotonically under both boundary conditions and the MBP is also preserved perfectly so that the solution is always located in the interval $(-1,1)$ .

Refer to caption — Fig. 1: The simulated solution $u$ subject to periodic boundary condition at $t=1$ , $5$ , $8$ , $30$ , $100$ , and $160$ respectively (left to right and top to bottom) in Example 5.1.

Example 5.2.

Consider the vector-valued equation (58) of the unknown vector field $\bm{u}:\Omega\subset\mathbb{R}^{2}\to\mathbb{R}^{2}$ with $\Delta$ replaced by $0.005\Delta$ . We test the domain

\Omega=\{(x,y)\in\mathbb{R}^{2}\,|\,\sqrt{x^{2}+y^{2}}<1\}\setminus\{(x,y)\in\mathbb{R}^{2}\,|\,\sqrt{(x-0.2)^{2}+y^{2}}\leq 0.5\},

that is, a region inside the unit disk but outside a circle with a shifted center. The Dirichlet boundary condition is set to be

\bm{u}(x,y)=\begin{dcases}(x,y),&\text{if }\sqrt{x^{2}+y^{2}}=1,\\ (2x-0.4,-2y),&\text{if }\sqrt{(x-0.2)^{2}+y^{2}}=0.5,\end{dcases}

i.e., the values of the vector field $\bm{u}$ on the outside boundary are always fixed to be a unit vector in the direction of $(x,y)$ and on the inside boundary be a unit vector in the direction of $(x,-y)$ . Thus the winding number of the boundary is $2$ .

We adopt the $C^{0}$ finite element spatial discretization with piecewise linear basis functions on triangular meshes and the mass-lumping technique in this example. The stabilizing constant is set to be $\kappa=2$ . PHIPM [77] is used for computing linear combinations of the products of the $\varphi$ -functions with vectors in the ETDRK2 scheme. We generate a triangular mesh with $2210$ nodes and $4158$ elements for the domain $\Omega$ and the initial configuration of the vector field $\bm{u}$ on the interior nodes with the fixed magnitude $0.9$ but random directions according to a uniform distribution. Fig. 5 shows the simulated vector fields at $t=0.1$ , $0.5$ , $1.2$ , $2.5$ , $15$ , and $100$ respectively. We observe that the initial disordered state quickly transitions into a more orderly structure which then asymptotically evolves to a steady state. The obtained steady state of the vector field $\bm{u}$ contain two vortices/defects which are symmetric with respect to the $x$ -axis. The evolution of the energy is plotted in the left graph of Fig. 6 for this example and we see the energy decreases monotonically in time. The right figure in Fig. 6 presents the evolution of the maximum value of $|\bm{u}|$ over the interior nodes and it demonstrates again that the MBP is perfectly preserved.

6 Concluding remarks

In this work, we establish an abstract mathematical framework for studying MBPs of a class of semilinear parabolic equations subject to a variety of boundary conditions, as well as unconditionally MBP-preserving temporal approximation schemes, ETD1 and ETDRK2, based on the exponential integrator with a suitably chosen generator of a contraction semigroup. The motivation is to reveal essential characteristics of the semilinear parabolic equations having the MBPs and to analyze the fundamental conditions under which the ETD schemes unconditionally preserve the MBPs of the underlying problems. We conclude that, to ensure the MBP of the model equation (5) and the ETD schemes, a crucial condition is the dissipation property of the linear operator stated in Assumption 1 and the sign change property of the nonlinear term stated in Assumption 2. The main results presented in this paper significantly generalize those in [20] in many aspects. The existence of MBP-preserving ETDRK schemes of even higher-order defined in other approaches is an open question and remains as one of our future works.

Extensions of the framework to the cases of complex-valued, real vector-valued and matrix-valued equations in the space-continuous setting are also carried out by considering the nonlinear operators taking the double-well-like forms and either Dirichlet, homogeneous Neumann, or periodic boundary condition, where the MBPs are proposed with respect to the vector and matrix $2$ -norms, respectively. We also note that whether the matrix-valued equation (64) with nonhomogeneous Dirichlet boundary condition has the MBP with respect to the matrix $2$ -norm is still an open question and needs to be further studied. The main difficulty resides in that the matrix $2$ -norm of the matrix case behaves essentially differently from the absolute value in the scalar case or the magnitude in the vector case. Moreover, it is worth noting that the matrix-valued equation (64) is derived in [79] as the $L^{2}$ gradient flow under some energy functional with the nonlinear part as a penalty term for matrix-valued fields which do not take orthogonal matrix values. The penalty term takes the form $|I_{m}-U^{T}U|_{F}^{2}$ , measuring the difference between a matrix $U$ and an orthogonal matrix in the sense of Frobenius norm ( $F$ -norm), so it seems more natural to consider the MBP with respect to the $F$ -norm. Actually, it is not hard to prove that, if the $F$ -norm of the initial data is not greater than $\sqrt{m}$ , then nor is the solution of the equation (64) at any time. However, unlike the case of $2$ -norm, the boundedness of the $F$ -norm of the solution is not sufficient to bound the nonlinear part in the splitting form; thus, whether this $F$ -norm-based MBP could be preserved by the ETD schemes is still an open question.

It is also worth mentioning that, apart from the ETD method, the integrating factor (IF) method is also a widely-used temporal integration method based on exponential integrators. While the ETD method only approximates the nonlinear terms as mentioned above, the IF method often applies numerical quadrature rules to the whole integrand. The IF method was introduced by Lawson [64] to solve ODE systems with large Lipschitz constants and applied to some problems with stiff linear part and nonstiff nonlinear part, such as the reaction-diffusion problems [56, 68, 76] and the advection-diffusion problems [53, 70, 106]. Due to highly different decaying rates of the exponential integrator components, the ETD schemes are usually more accurate than the IF ones for highly stiff systems. It is also the case that, if the nonlinear function is given by a constant (e.g., $f[u]\equiv c$ ), the ETD schemes can produce the exact solution to (1) while only approximate solutions are obtained by the IF schemes. One may ask whether, similar to ETD schemes, the IF schemes could preserve the MBPs. Related to this question, the authors of [53] focused on the property of strong stability-preserving (SSP) [42] for the IF Runge–Kutta (IFRK) schemes. SSP is a stronger stability than the MBP considered here. In fact, if we weaken the assumptions given in [53] appropriately, we could obtain MBP-preserving IFRK schemes under some suitable constraint on the time step size. Recall that a scheme is SSP means that if the linear operator ${\mathcal{L}}$ satisfies

(70)

\|\mathrm{e}^{\tau{\mathcal{L}}}\|\leq 1,\quad\forall\,\tau>0,

which is exactly Lemma 1-(ii), and there exists some $\tau_{0}>0$ such that the nonlinear mapping $f$ satisfies

(71)

\|w+\tau f[w]\|\leq\|w\|,\quad\forall\,w\in\mathcal{X},\ \tau\leq\tau_{0},

then the solution always satisfies $\|v^{n+1}\|\leq\|v^{n}\|$ for any $\tau\leq\tau_{0}$ . Instead of (71), if one makes an assumption on $f$ as follows:

(72)

\|w+\tau f[w]\|\leq\beta,\quad\forall\,\text{$w\in\mathcal{X}$ with $\|w\|\leq\beta$ and $\tau\leq\tau_{0}$},

which is actually equivalent to Lemma 2-(i) with $\kappa=1/\tau$ , With slight modifications to the proof for the IFRK schemes in [53], one can also conclude that the IFRK schemes, which are SSP under the conditions (70) and (71), are also MBP-preserving under the assumptions (70) and (72) (or say, Assumptions 1 and 2). Note that such property of MBP-preserving is only conditional in the sense that some constraint on the time step size is necessary. An open question is whether the MBP could be preserved unconditionally by the IF schemes as the ETD schemes if an appropriate stabilizer is introduced and this also remains as one of our future works.

Acknowledgements

The authors would like to thank the anonymous reviewers for their constructive comments and suggestions, which have helped us greatly improve this work. The authors are also grateful to Dr. Xiaochuan Tian of University of Texas at Austin and Dr. Zhi Zhou of The Hong Kong Polytechnic University for many valuable discussions.

References

[1] S. M. Allen and J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metall., 27 (1979), pp. 1085–1095.
[2] H. Amann, Invariant sets and existence theorems for semilinear parabolic and elliptic systems, J. Math. Anal. Appl., 65 (1978), pp. 432–467.
[3] S. Badia and A. Hierro, On discrete maximum principles for discontinuous Galerkin methods, Comput. Methods Appl. Mech. Engrg., 286 (2015), pp. 107–122.
[4] G. Beylkin, J. M. Keiser, and L. Vozovoi, A new class of time discretization schemes for the solution of nonlinear PDEs, J. Comput. Phys., 147 (1998), pp. 362–387.
[5] J. H. Bramble and B. E. Hubbard, New monotone type approximations for elliptic problems, Math. Comp., 18 (1964), pp. 349–367.
[6] J. H. Brandts, S. Korotov, and M. Křížek, The discrete maximum principle for linear simplicial finite element approximations of a reaction-diffusion problem, Linear Algebra Appl., 429 (2008), pp. 2344–2357.
[7] E. Burman and A. Ern, Discrete maximum principle for Galerkin approximations of the Laplace operator on arbitrary meshes, C. R. Math. Acad. Sci. Paris, Ser. I 338 (2004), pp. 641–646.
[8] E. Chasseigne, The Dirichlet problem for some nonlocal diffusion equations, Differential Integral Equations, 20 (2007), pp. 1389–1404.
[9] P. G. Ciarlet, Discrete maximum principle for finite-difference operators, Aequations Math., 4 (1970), pp. 338–352.
[10] P. G. Ciarlet and P. A. Raviart, Maximum principle and uniform convergence for the finite element method, Comput. Methods Appl. Mech. Engrg., 2 (1973), pp. 17–31.
[11] S. M. Cox and P. C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys., 176 (2002), pp. 430–455.
[12] Q. Du, Global existence and uniqueness of solutions of the time-dependent Ginzburg–Landau model for superconductivity, Appl. Anal., 53 (1994), pp. 1–17.
[13] Q. Du, Discrete gauge invariant approximations of a time-dependent Ginzburg–Landau model of superconductivity, Math. Comp., 67 (1998), pp. 965–986.
[14] Q. Du, Numerical approximations of the Ginzburg–Landau models for superconductivity, J. Math. Phys., 46 (2005), Art. 095109.
[15] Q. Du, Nonlocal Modeling, Analysis, and Computation, Vol. 94 of CBMS-NSF regional conference series in Applied Mathematics, SIAM, 2019.
[16] Q. Du and X. B. Feng, The phase field method for geometric moving interfaces and their numerical approximations, in Geometric Partial Differential Equations - Part I, Handbook of Numerical Analysis, Vol. 21, Elsevier, 2020.
[17] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou, Analysis and approximation of nonlocal diffusion problems with volume constraints, SIAM Rev., 54 (2012), pp. 667–696.
[18] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou, A nonlocal vector calculus, nonlocal volume-constrained problems, and nonlocal balance laws, Math. Models Methods Appl. Sci., 23 (2013), pp. 493–540.
[19] Q. Du, M. Gunzburger, and J. Peterson, Analysis and approximation of Ginzburg–Landau models for superconductivity, SIAM Rev., 34 (1992), pp. 54–81.
[20] Q. Du, L. Ju, X. Li, and Z. H. Qiao, Maximum principle preserving exponential time differencing schemes for the nonlocal Allen–Cahn equation, SIAM J. Numer. Anal., 57 (2019), pp. 875–898.
[21] Q. Du, R. Nicolaides, and X. Wu, Analysis and convergence of a covolume approximation of the Ginzburg–Landau model of superconductivity, SIAM J. Numer. Anal., 35 (1998), pp. 1049–1072.
[22] Q. Du, Y. Z. Tao, X. C. Tian, and J. Yang, Asymptotically compatible discretization of multidimensional nonlocal diffusion models and approximation of nonlocal Green’s functions, IMA J. Numer. Anal., 39 (2019), pp. 607–625.
[23] Q. Du and X. B. Yin, A conforming DG method for linear nonlocal models with integrable kernels, J. Sci. Comput., 80 (2019), pp. 1913–1935.
[24] Q. Du and W.-X. Zhu, Analysis and applications of the exponential time differencing schemes and their contour integration modifications, BIT Numer. Math., 45 (2005), pp. 307–328.
[25] S. W. Duo, H. W. van Wyk, and Y. Z. Zhang, A novel and accurate finite difference method for the fractional Laplacian and the fractional Poisson problem, J. Comput. Phys., 355 (2018), pp. 233–252.
[26] S. W. Duo and Y. Z. Zhang, Accurate numerical methods for two and three dimensional integral fractional Laplacian with applications, Comput. Methods Appl. Mech. Engrg., 355 (2019), pp. 639–662.
[27] K.-J. Engel and R. Nagel, One-Parameter Semigroups for Linear Evolution Equations, Graduate Texts in Mathematics, Vol. 194, Springer-Verlag, New York, 2000.
[28] L. C. Evans, Partial Differential Equations, American Mathematical Society, Providence, Rhode Island, 2000.
[29] L. C. Evans, H. M. Soner, and P. E. Souganidis, Phase transitions and generalized motion by mean curvature, Commun. Pure Appl. Math., 45 (1992), pp. 1097–1123.
[30] I. Faragó and R. Horváth, Discrete maximum principle and adequate discretizations of linear parabolic problems, SIAM J. Sci. Comput., 28 (2006), pp. 2313–2336.
[31] I. Faragó, J. Karátson, and S. Korotov, Discrete maximum principles for nonlinear parabolic PDE systems, IMA J. Numer. Anal., 32 (2012), pp. 1541–1573.
[32] I. Faragó, S. Korotov, and T. Szabó, On continuous and discrete maximum principles for elliptic problems with the third boundary condition, Appl. Math. Comput., 219 (2013), pp. 7215–7224.
[33] I. Faragó, S. Korotov, and T. Szabó, On continuous and discrete maximum-minimum principles for reaction-diffusion problems with the Neumann boundary condition, Applications of mathematics, pp. 34–44, Czech. Acad. Sci., Prague, 2015.
[34] X. B. Feng and A. Prohl, Numerical analysis of the Allen–Cahn equation and approximation for mean curvature flows, Numer. Math., 94 (2003), pp. 33–65.
[35] X. B. Feng and A. Prohl, Error analysis of a mixed finite element method for the Cahn–Hilliard equation, Numer. Math., 99 (2004), pp. 47–84.
[36] X. Fernández-Real and X. Ros-Oton, Boundary regularity for the fractional heat equation, Rev. R. Acad. Cienc. Exactas Fís. Nat. Ser. A Mat., 110 (2016), pp. 49–64.
[37] M. Frittelli, A. Madzvamuse, I. Sgura, and C. Venkataraman, Preserving invariance properties of reaction-diffusion systems on stationary surfaces, IMA J. Numer. Anal., 39 (2019), pp. 235–270.
[38] E. Gallopoulous and Y. Saad, Efficient solution of parabolic equations by Krylov approximation methods, SIAM J. Sci. Comput., 13 (1992), pp. 1236–1264.
[39] S. Gaudreault and J. A. Pudykiewicz, An efficient exponential time integration method for the numerical solution of the shallow water equations on the sphere, J. Comput. Phys., 322 (2016), pp. 827–848.
[40] S. Gaudreault, G. Rainwater, and M. Tokman, KIOPS: A fast adaptive Krylov subspace solver for exponential integrators, J. Comput. Phys., 372 (2018), pp. 236–255.
[41] V. Ginzburg and L. Landau, On the theory of superconductivity, Zh. Eksperim. Teor. Fiz., 20 (1950), pp. 1064–1082.
[42] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
[43] J.-L. Guermond, B. Popov, and L. Saavedra, Second-order invariant domain preserving ALE approximation of hyperbolic systems, J. Comput. Phys., 401 (2020), Art. 108927.
[44] J.-L. Guermond, B. Popov, and I. Tomas, Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems, Comput. Methods Appl. Mech. Engrg., 347 (2019), pp. 143–175.
[45] N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, PA, 2008.
[46] M. Hochbruck and C. Lubich, On Krylov subspace approximations to the matrix exponential operator, SIAM J. Numer. Anal., 34 (1997), pp. 1911–1925.
[47] M. Hochbruck, C. Lubich, and H. Selhofer, Exponential integrators for large systems of differential equations, SIAM J. Sci. Comput., 19 (1998), pp. 1552–1574.
[48] M. Hochbruck and A. Ostermann, Explicit exponential Runge–Kutta methods for semilinear parabolic problems, SIAM J. Numer. Anal., 43 (2005), pp. 1069–1090.
[49] M. Hochbruck and A. Ostermann, Exponential integrators, Acta Numer., 19 (2010), pp. 209–286.
[50] R. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge University Press, Cambridge, 1991.
[51] T. L. Hou, T. Tang, and J. Yang, Numerical analysis of fully discretized Crank–Nicolson scheme for fractional-in-space Allen–Cahn equations, J. Sci. Comput., 72 (2017), pp. 1214–1231.
[52] Y. H. Huang and A. Oberman, Numerical methods for the fractional Laplacian: a finite difference-quadrature approach, SIAM J. Numer. Anal., 52 (2014), pp. 3056–3084.
[53] L. Isherwood, Z. J. Grant, and S. Gottlieb, Strong stability preserving integrating factor Runge–Kutta methods, SIAM J. Numer. Anal., 56 (2018), pp. 3276–3307.
[54] Y. Jiang and H. L. Liu, Invariant-region-preserving DG methods for multi-dimensional hyperbolic conservation law systems, with an application to compressible Euler equations, J. Comput. Phys., 373 (2018), pp. 385–409.
[55] L. Ju, X. Li, Z. H. Qiao, and H. Zhang, Energy stability and error estimates of exponential time differencing schemes for the epitaxial growth model without slope selection, Math. Comp., 87 (2018), pp. 1859–1885.
[56] L. Ju, X. F. Liu, And W. Leng, Compact implicit integration factor methods for a family of semilinear fourth-order parabolic equations, Discrete Contin. Dyn. Syst. Ser. B, 19 (2014), pp. 1667–1687.
[57] L. Ju, J. Zhang, and Q. Du, Fast and accurate algorithms for simulating coarsening dynamics of Cahn–Hilliard equations, Comput. Mater. Sci., 108 (2015), pp. 272–282.
[58] L. Ju, J. Zhang, L. Y. Zhu, and Q. Du, Fast explicit integration factor methods for semilinear parabolic equations, J. Sci. Comput., 62 (2015), pp. 431–455.
[59] A. Karafiat, Discrete maximum principle in parabolic boundary-value problems, Ann. Polon. Math., 53 (1991), pp. 253–265.
[60] J. Karàtson and S. Korotov, Some discrete maximum principles arising for nonlinear elliptic finite element problems, Comput. Math. Appl., 70 (2015), pp. 2732–2741.
[61] J. Karàtson, S. Korotov, and M. Křížek, On discrete maximum principles for nonlinear elliptic problems, Math. Comput. Simulation, 76 (2007), pp. 99–108.
[62] D. R. Kassoy and J. Poland, The thermal explosion confined by a constant temperature boundary, SIAM J. Appl. Math., 39 (1980), pp. 412–430.
[63] H. J. Kuiper, Invariant sets for nonlinear elliptic and parabolic systems, SIAM J. Math. Anal., 11 (1980), pp. 1075–1103.
[64] J. D. Lawson, Generalized Runge–Kutta processes for stable systems with large Lipschitz constants, SIAM J. Numer. Anal., 4 (1967), pp. 372–380.
[65] D. Li, Z. H. Qiao, and T. Tang, Characterizing the stabilization size for semi-implicit Fourier-spectral method to phase field equations, SIAM J. Numer. Anal., 54 (2016), pp. 1653–1681.
[66] R. H. Li, Z. Y. Chen, and W. Wu, Generalized Difference Methods for Differential Equations: Numerical Analysis of Finite Volume Methods, Marcel Dekker, Inc., 2000.
[67] Z. L. Li and K. Ito, Maximum principle preserving schemes for interface problems with discontinuous coefficients, SIAM J. Sci. Comput., 23 (2001), pp. 339–361.
[68] X. F. Liu and Q. Nie, Compact integration factor methods for complex domains and adaptive mesh refinement, J. Comput. Phys., 229 (2010), pp. 5692–5706.
[69] H. L. Liu and H. Yu, Maximum-principle-satisfying third order discontinuous Galerkin schemes for Fokker–Planck equations, SIAM J. Sci. Comput., 36 (2014), pp. A2296–A2325.
[70] D. Lu and Y. T. Zhang, Computational complexity study on Krylov integration factor WENO method for high spatial dimension convection-diffusion problems, J. Sci. Comput., 73 (2017), pp. 980–1027.
[71] E. C. MacRae, Matrix derivatives with an application to an adaptive linear decision problem, Ann. Statist., 2 (1974), pp. 337–346.
[72] C. Mastroserio and M. Montrone, Invariant regions and asymptotic behaviour for the numerical solution of reaction-diffusion systems by a class of alternating direction methods, Calcolo, 21 (1984), pp. 269–279.
[73] M. E. Mincsovics and T. L. Horváth, On the differences of the discrete weak and strong maximum principles for elliptic operators, Large-scale scientific computing, pp. 614–621, Lecture Notes in Comput. Sci., 7116, Springer, Heidelberg, 2012.
[74] C. Moler and C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, SIAM Rev., 20 (1978), pp. 801–836.
[75] C. Moler and C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev., 45 (2003), pp. 3–49.
[76] Q. Nie, F. Y. M. Wan, Y. T. Zhang, and X. F. Liu, Compact integration factor methods in high spatial dimensions, J. Comput. Phys., 227 (2008), pp. 5238–5255.
[77] J. Niesen and W. M. Wright, Algorithm 919: A Krylov subspace algorithm for evaluating the $\phi$ -functions appearing in exponential integrators, ACM Trans. Math. Soft., 38 (2008), Art. 22.
[78] R. H. Nochetto and W. J. Zhang, Discrete ABP estimate and convergence rates for linear elliptic equations in non-divergence form, Found. Comput. Math., 18 (2018), pp. 537–593.
[79] B. Osting and D. Wang, A diffusion generated method for orthogonal matrix-valued fields, Math. Comp., 89 (2020), pp. 515–550.
[80] E.-M. Ouhabaz, $L^{\infty}$ -contractivity of semigroups generated by sectorial forms, J. London Math. Soc., 46 (1992), pp. 529–542.
[81] C. V. Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, New York, 1992.
[82] D.-Y. Peng and D. B. Robinson, A new two-constant equation of state, Ind. Eng. Chem. Fundamen., 15 (1976), pp. 59–64.
[83] H. S. Price, Monotone and oscillation matrices applied to finite difference approximations, Math. Comp., 22 (1968), pp. 489–516.
[84] Z. H. Qiao and S. Y. Sun, Two-phase fluid simulation using a diffuse interface model with Peng–Robinson equation of state, SIAM J. Sci. Comput., 36 (2014), pp. B708–B728.
[85] R. Redheffer and W. Walter, Invariant sets for systems of partial differential equations. I. Parabolic equations, Arch. Rational Mech. Anal., 67 (1978), pp. 41–52.
[86] W. Rudin, Principles of Mathematical Analysis, Third Edition, McGraw-Hill Book Co., New York, 1976.
[87] S. G. Samko, A. A. Kilbas, and O. I. Marichev, Fractional Integrals and Derivatives, Gordon and Breach Science Publishers, Yverdon, 1993.
[88] J. Shen, T. Tang, and J. Yang, On the maximum principle preserving schemes for the generalized Allen–Cahn equation, Commun. Math. Sci., 14 (2016), pp. 1517–1534.
[89] J. Shen and X. F. Yang, Numerical approximations of Allen–Cahn and Cahn–Hilliard equations, Discrete Contin. Dyn. Syst., 28 (2010), pp. 1669–1691.
[90] R. B. Sidje, Expokit: Software package for computing matrix exponentials, ACM Trans. Math. Soft., 24 (1998), pp. 130–156.
[91] J. Smoller, Shock Waves and Reaction-Diffusion Equations, Fundamental Principles of Mathematical Sciences, Vol. 258, Second Edition, Springer-Verlag, New York, 1994.
[92] P. Stehlík and J. Volek, Maximum principles for discrete and semidiscrete reaction-diffusion equation, Discrete Dyn. Nat. Soc., 2015, Art. 791304.
[93] T. Tang and J. Yang, Implicit-explicit scheme for the Allen–Cahn equation preserves the maximum principle, J. Comput. Math., 34 (2016), pp. 471–481.
[94] H. Tian, L. Ju, and Q. Du, Nonlocal convection-diffusion problems and finite element approximations, Comput. Methods Appl. Mech. Engrg., 289 (2015), pp. 60–78.
[95] X. C. Tian and Q. Du, Analysis and comparison of different approximations to nonlocal diffusion and linear peridynamic equations, SIAM J. Numer. Anal., 51 (2013), pp. 3458–3482.
[96] X. C. Tian, Q. Du, and M. Gunzburger, Asymptotically compatible schemes for the approximation of fractional Laplacian and related nonlocal diffusion problems on bounded domains, Adv. Comput. Math., 42 (2016), pp. 1363–1380.
[97] R. S. Varga, On a discrete maximum principle, SIAM J. Numer. Anal., 3 (1966), pp. 355–359.
[98] C. Wang, S. M. Wise, and J. S. Lowengrub, An energy-stable and convergent finite-difference scheme for the phase field crystal equation, SIAM J. Numer. Anal., 47 (2009), pp. 2269–2288.
[99] X. Q. Wang, L. Ju, and Q. Du, Efficient and stable exponential time differencing Runge–Kutta methods for phase field elastic bending energy models, J. Comput. Phys., 316 (2016), pp. 21–38.
[100] J. Yang, Q. Du, and W. Zhang, Uniform $L^{p}$ -bound of the Allen–Cahn equation and its numerical discretization, Int. J. Numer. Anal. Mod., 18 (2018), pp. 213–227.
[101] P. Yang, T. Xiong, J. M. Qiu, and Z. F. Xu, High order maximum principle preserving finite volume method for convection dominated problems, J. Sci. Comput., 67 (2016), pp. 795–820.
[102] E. G. Yanik, A discrete maximum principle for collocation methods, Comput. Math. Appl., 14 (1987), pp. 459–464.
[103] E. G. Yanik, Sufficient conditions for a discrete maximum principle for high order collocation methods, Comput. Math. Appl., 17 (1989), pp. 1431–1434.
[104] G. W. Yuan and Y. L. Yu, Existence of solution of a finite volume scheme preserving maximum principle for diffusion equations, Numer. Meth. Part. Diff. Eq., 34 (2018), pp. 80–96.
[105] J. Zhang, C. B. Zhou, Y. G. Wang, L. Ju, Q. Du, X. B. Chi, D. S. Xu, D. X. Chen, Y. Liu, and Z. Liu, Extreme-scale phase field simulations of coarsening dynamics on the Sunway Taihulight supercomputer, in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC’16), 2016, Art. 4.
[106] S. Zhao, J. Ovadia, X. F. Liu, Y. T. Zhang, and Q, Nie, Operator splitting implicit integration factor methods for stiff reaction-diffusion-advection systems, J. Comput. Phys., 230 (2011), pp. 5996–6009.
[107] L. Y. Zhu, L. Ju, and W. D. Zhao, Fast high-order compact exponential time differencing Runge–Kutta methods for second-order semilinear parabolic equations, J. Sci. Comput., 67 (2016), pp. 1043–1065.

	$\displaystyle\\|w(t)-\tilde{w}(t)\\|$	$\displaystyle\leq\int_{0}^{t}\mathrm{e}^{-\kappa(t-s)}\|\!\|\!\|S_{{\mathcal{L}}_{0}}(t-s)\|\!\|\!\|\\|\mathcal{N}[v(s)]-\mathcal{N}[\tilde{v}(s)]\\|\,\mathrm{d}s$
		$\displaystyle\leq 2\kappa\int_{0}^{t}\mathrm{e}^{-\kappa(t-s)}\\|v(s)-\tilde{v}(s)\\|\,\mathrm{d}s$
		$\displaystyle\leq 2(1-\mathrm{e}^{-\kappa t_{1}})\\|v-\tilde{v}\\|_{C([0,t_{1}];\mathcal{X})},\quad t\in[0,t_{1}],$

	$\displaystyle\\|e_{1}^{n+1}\\|$	$\displaystyle\leq\mathrm{e}^{-\kappa\tau}\|\!\|\!\|\mathrm{e}^{\tau{\mathcal{L}}_{h0}}\|\!\|\!\|\\|e_{1}^{n}\\|$
		$\displaystyle\quad+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\|\!\|\!\|\mathrm{e}^{(\tau-s){\mathcal{L}}_{h0}}\|\!\|\!\|\big{\{}\\|\mathcal{N}[v^{n}]-\mathcal{N}[u(t_{n})]\\|+\\|R_{1}(s)\\|\big{\}}\,\mathrm{d}s$
		$\displaystyle\leq\mathrm{e}^{-\kappa\tau}\\|e_{1}^{n}\\|+(2\kappa\\|e_{1}^{n}\\|+C_{1}\kappa\tau)\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\,\mathrm{d}s$
		$\displaystyle=\mathrm{e}^{-\kappa\tau}\\|e_{1}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa}(2\kappa\\|e_{1}^{n}\\|+C_{1}\kappa\tau)$
		$\displaystyle=(2-\mathrm{e}^{\kappa\tau})\\|e_{1}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa\tau}\cdot C_{1}\kappa\tau^{2}$
(53)			$\displaystyle\leq(1+\kappa\tau)\\|e_{1}^{n}\\|+C_{1}\kappa\tau^{2},$

	$\displaystyle\\|e_{1}^{n}\\|$	$\displaystyle\leq(1+\kappa\tau)^{n}\\|e_{1}^{0}\\|+C_{1}\kappa\tau^{2}\sum_{k=0}^{n-1}(1+\kappa\tau)^{k}$
		$\displaystyle=(1+\kappa\tau)^{n}\\|e_{1}^{0}\\|+C_{1}\tau[(1+\kappa\tau)^{n}-1]$
		$\displaystyle\leq\mathrm{e}^{\kappa n\tau}\\|e_{1}^{0}\\|+C_{1}\mathrm{e}^{\kappa n\tau}\tau.$

	$\displaystyle\\|e_{2}^{n+1}\\|$	$\displaystyle\leq\mathrm{e}^{-\kappa\tau}\\|e_{2}^{n}\\|+\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}[2\kappa(1+\kappa s)\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau s+(C_{2}\kappa+C_{3})\tau^{2}]\,\mathrm{d}s$
		$\displaystyle=\mathrm{e}^{-\kappa\tau}\\|e_{2}^{n}\\|+[2\kappa\\|e_{2}^{n}\\|+(C_{2}\kappa+C_{3})\tau^{2}]\int_{0}^{\tau}\mathrm{e}^{-\kappa(\tau-s)}\,\mathrm{d}s$
		$\displaystyle\quad+(2\kappa^{2}\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau)\int_{0}^{\tau}s\mathrm{e}^{-\kappa(\tau-s)}\,\mathrm{d}s$
		$\displaystyle=\mathrm{e}^{-\kappa\tau}\\|e_{2}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa}[2\kappa\\|e_{2}^{n}\\|+(C_{2}\kappa+C_{3})\tau^{2}]$
		$\displaystyle\quad+\frac{\mathrm{e}^{-\kappa\tau}-1+\kappa\tau}{\kappa^{2}}(2\kappa^{2}\\|e_{2}^{n}\\|+2C_{1}\kappa^{2}\tau)$
		$\displaystyle=(\mathrm{e}^{-\kappa\tau}+2\kappa\tau)\\|e_{2}^{n}\\|+\frac{1-\mathrm{e}^{-\kappa\tau}}{\kappa}(C_{2}\kappa+C_{3})\tau^{2}+\frac{\mathrm{e}^{-\kappa\tau}-1+\kappa\tau}{\kappa^{2}}\cdot 2C_{1}\kappa^{2}\tau$
		$\displaystyle\leq\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}\\|e_{2}^{n}\\|+(C_{2}\kappa+C_{3}+C_{1}\kappa^{2})\tau^{3},$

	$\displaystyle\\|e_{2}^{n}\\|$	$\displaystyle\leq\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{n}\\|e_{2}^{0}\\|+(C_{1}\kappa^{2}+C_{2}\kappa+C_{3})\tau^{3}\sum_{k=0}^{n-1}\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{k}$
		$\displaystyle\leq\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{n}\\|e_{2}^{0}\\|+\Big{(}C_{1}\kappa+C_{2}+\frac{C_{3}}{\kappa}\Big{)}\tau^{2}\Big{[}\Big{(}1+\kappa\tau+\frac{(\kappa\tau)^{2}}{2}\Big{)}^{n}-1\Big{]}$
		$\displaystyle\leq\mathrm{e}^{\kappa n\tau}\\|e_{2}^{0}\\|+\Big{(}C_{1}\kappa+C_{2}+\frac{C_{3}}{\kappa}\Big{)}\mathrm{e}^{\kappa n\tau}\tau^{2}.$

Maximum bound principles for a class of semilinear parabolic equations and exponential time differencing schemes

Abstract

keywords:

AMS:

1 Introduction

2 Maximum bound principle of the model equation

2.1 Abstract framework

Assumption 1 (Requirements on ℒ{\mathcal{L}} and ℒ0{\mathcal{L}}_{0}).

Assumption 2 (Requirements on ff).

Remark 1.

Lemma 1.

Proof.

Remark 2.

Lemma 2.

Proof.

Theorem 3.

Proof.

Remark 3.

2.2 Examples of the nonlinear function f0f_{0}

Example 2.1.

Example 2.2.

Example 2.3.

2.3 Examples of the linear operator ℒ{\mathcal{L}}

2.3.1 Infinite dimensional examples

Example 2.4 (Second-order elliptic differential operator).

Remark 4.

Example 2.5 (Nonlocal diffusion operator).

Example 2.6 (Fractional Laplace operator).

Corollary 4.

Remark 5.

2.3.2 Finite dimensional examples

Proposition 5.

Example 2.7 (Central difference operator for the Laplacian).

Remark 6.

Example 2.8 (Quadrature-based difference operator for the nonlocal diffusion).

Example 2.9 (Finite difference operator for the fractional Laplacian).

Example 2.10 (Finite element operator for the Laplacian).

Remark 7.

Corollary 6.

2.4 Discussion on homogeneous Neumann boundary conditions

3 Exponential time differencing for temporal approximation

3.1 ETD schemes and discrete MBPs

Theorem 7 (Maximum bound principle of the ETD1 scheme).

Proof.

Theorem 8 (Maximum bound principle of the ETDRK2 scheme).

Proof.

Corollary 9.

3.2 Convergence analysis of the ETD schemes

Theorem 10.

Proof.

Remark 8.

Theorem 11.

Proof.

3.3 Energy stability of the ETD schemes for gradient flow models

Proposition 12.

4 Some extensions

4.1 Extension to complex scalar-valued and real vector-valued equations

Proposition 13.

Lemma 14.

Proof.

Proposition 15.

Theorem 16.

Proof.

Remark 9.

4.2 Extension to real matrix-valued equations

Lemma 17.

Proof.

Lemma 18.

Proof.

Proposition 19.

Proof.

5 Numerical experiments

5.1 Implementations of matrix exponentials

Remark 10.

5.2 Experimental results

Example 5.1.

Example 5.2.

6 Concluding remarks

Acknowledgements

References

Assumption 1 (Requirements on ${\mathcal{L}}$ and ${\mathcal{L}}_{0}$ ).

Assumption 2 (Requirements on $f$ ).

2.2 Examples of the nonlinear function $f_{0}$

2.3 Examples of the linear operator ${\mathcal{L}}$