$hp$-FEM for the fractional heat equation

TL;DR

This paper develops an $hp$-Finite Element method combined with a Caffarelli-Silvestre extension and $hp$-Discontinuous Galerkin timestepping to solve the fractional heat equation, proving exponential convergence.

Contribution

It introduces a novel discretization scheme for the fractional heat equation using $hp$-FEM and DG timestepping, with a rigorous convergence proof.

Findings

Exponential convergence of the proposed method.

Effective discretization of nonlocal fractional operators.

Rigorous theoretical framework for the scheme.

Abstract

We consider a time dependent problem generated by a nonlocal operator in space. Applying a discretization scheme based on $hp$-Finite Elements and a Caffarelli-Silvestre extension we obtain a semidiscrete semigroup. The discretization in time is carried out by using $hp$-Discontinuous Galerkin based timestepping. We prove exponential convergence for such a method in an abstract framework for the discretization in the original domain $\Omega$.

Equations

$$\dot{u}(t)+{\mathcal{L}}^{s}u(t)$$

$$u(\cdot ,t)$$

$${\mathcal{L}}^{s}u:=\sum _{j=0}^{\infty }{\mu }_j^s(u,{\phi }_j{)}_{L^2(\Omega )}{\phi }_j,$$

$$-\mathrm{div}\left(y^{\alpha }A\nabla \mathcal{U}\right)+y^{\alpha }c\mathcal{U}$$

$$d_s\mathrm{tr}\dot{\mathcal{U}}+{\partial }_{\nu }^{\alpha }\mathcal{U}$$

$$\mathcal{U}$$

$${\partial }_{\nu }^{\alpha }\mathcal{U}:=-\underset{y\to 0^{+}}{\lim }{y}^{\alpha }{\partial }_y\mathcal{U}(\cdot ,y),\text{and }\mathrm{tr}\mathcal{U}:=\mathcal{U}(\cdot ,0)$$

$$L^2(y^{\alpha },D)$$

$$H^1(y^{\alpha },D)$$

$${\mathring{H}}^1(y^{\alpha },D)$$

$$\displaystyle\mathcal{A}(\mathscr{U},\mathscr{V}):=\int_{\mathcal{C}}{y^{\alpha}\big{(}A\nabla\mathscr{U}\big{)}\cdot\nabla\mathscr{V}+y^{\alpha}c\,\mathscr{U}\mathscr{V}}.$$

$${\Vert u\Vert }_{[X_0,X_1{]}_{\theta ,2}}^2$$

$${\left[X_0,X_1\right]}_{\theta ,2}$$

$${\tilde{H}}^s(\Omega ):={\left[L^2(\Omega ),H_0^1(\Omega )\right]}_{s,2},\text{and}{H}^s(\Omega )$$

$${\mathbb{H}}^s(\Omega ):=\left\{u\in L^2(\Omega ):\sum _{j=0}^{\infty }{\mu }_j^s|{\left(u,{\phi }_j\right)}_{L^2(\Omega )}{|}^2<\infty \right\}.$$

$${\mathcal{N}}_{\Omega }:=\dim ({\mathbb{V}}_h^{\mathcal{X}})\text{and}{\mathcal{N}}_{\mathcal{Y}}:=\dim ({\mathbb{V}}_h^{\mathcal{Y}}).$$

$${\left({\mathcal{L}}_h^{s}u,v\right)}_{L^2(\Omega )}$$

$$\mathcal{A}({\mathcal{D}}_{h}u,{\mathcal{V}}_h)$$

$$\mathrm{tr}{\mathcal{D}}_{h}u$$

$$\mathcal{A}(\mathcal{D}u,\mathcal{V})$$

$$\mathrm{tr}\mathcal{D}u$$

$${\Vert {\left(\lambda \mathrm{I}+{\mathcal{L}}_h^s\right)}^{-1}f\Vert }_{L^2(\Omega )}\le \frac{M}{1+\left|\lambda \right|}{\Vert f\Vert }_{L^2(\Omega )}$$

$${\left(\lambda d_s\mathrm{tr}{\mathcal{U}}_{\lambda },\mathrm{tr}{\mathcal{V}}_h\right)}_{L^2(\Omega )}+\mathcal{A}({\mathcal{U}}_{\lambda },{\mathcal{V}}_h)$$

$$\mathrm{Re}(\lambda )d_s{\Vert \mathrm{tr}{\mathcal{U}}_{\lambda }\Vert }^2+\mathcal{A}({\mathcal{U}}_{\lambda },{\mathcal{U}}_{\lambda })\le \left|\lambda d_s{\Vert \mathrm{tr}{\mathcal{U}}_{\lambda }\Vert }^2+\mathcal{A}({\mathcal{U}}_{\lambda },{\mathcal{U}}_{\lambda })\right|$$

$${\Vert u\Vert }_{{\mathbb{V}}_h^{\mathcal{X}}}:={\Vert {\mathcal{D}}_{h}u\Vert }_{{\mathring{H}}^1(y^{\alpha },\mathcal{C})},$$

$$c_1{\Vert u\Vert }_{{\mathbb{V}}_h^{\mathcal{X}}}^2\le {\left({\mathcal{L}}_h^{s}u,u\right)}_{L^2(\Omega )}\le c_2{\Vert u\Vert }_{{\mathbb{V}}_h^{\mathcal{X}}}^2.$$

$$c_3{\Vert u\Vert }_{{\tilde{H}}^s(\Omega )}^2\le {\left({\mathcal{L}}_h^{s}u,u\right)}_{L^2(\Omega )}.$$

$${\Vert u\Vert }_{{\tilde{H}}^s(\Omega )}^2$$

$${\left({\mathcal{L}}_h^{s}u,u\right)}_{L^2(\Omega )}$$

$${\dot{u}}_h+{\mathcal{L}}_h^s{u}_h$$

$$u_h(0)$$

$$u(t)$$

$${\mathbb{V}}_{h,\beta }^{\mathcal{X}}:={\left[\left({\mathbb{V}}_h^{\mathcal{X}},{\Vert \cdot \Vert }_{L^2(\Omega )}\right),\left({\mathbb{V}}_h^{\mathcal{X}},{\Vert \cdot \Vert }_{{\mathbb{V}}_h^{\mathcal{X}}}\right)\right]}_{\beta ,2}.$$

$${\Vert u_{h,0}\Vert }_{{\mathbb{V}}_{h,\beta }^{\mathcal{X}}}\le C_{stab}<\infty \text{ and }{\Vert u_0-u_{h,0}\Vert }_{L^2(\Omega )}\le C_{approx}{e}^{-b{\mathcal{N}}_{\Omega }^{\mu }},$$

$$\mathcal{S}:=\mathbb{C}\setminus \left[\left\{z_0+z:\left|\mathrm{Arg}(z)\right|\le \frac{\pi }{8},\mathrm{Re}(z)\ge 0\right\}\cup B_{{\epsilon }_0}(0)\right].$$

Full text

1 Introduction

For stationary fractional diffusion, numerical techniques have recently been proposed that provide exponential convergence of the error with respect to the computational effort, [BMN*+*18, BMS19]. The construction is based on $hp$-Finite Elements on appropriate geometric meshes. The purpose of the present article is to generalize these techniques to the time dependent setting. We consider the discretization of the time dependent problem (2.1), generated by a fractional power of an elliptic operator. The spatial discretization of the nonlocal operator is based on a reformulation using the Caffarelli-Silvestre extension, for which an $hp$-Finite Element discretization (FEM) is employed. The discretization in time is then carried out by a Discontinuous Galerkin method in the spirit of [SS00] of either fixed order or in its $hp$ version. Our analysis hinges on two conditions, one related to stable liftings of the initial condition and the second one related to the ability to approximate solutions of singularly perturbed problems.

After establishing the abstract framework, we work out the case of $hp$-FEM in the special case of 1D or 2D with analytic data and geometry and show that the basic Assumptions 3.5 and 3.9 are satisfied on appropriate geometric meshes. The reduction of scope to smooth geometries and at most 2D mainly is done to keep the presentation to a reasonable length; we expect that it is possible to establish the assumptions of the abstract framework also for the case of polygons or $\Omega\subseteq{\mathbb{R}}^{d}$, $d>2$.

Discretization schemes for the same model problem have already appeared in the literature. In [BLP17], the approximation is done by applying numerical quadrature to the Dunford-Taylor representation of the solution and using a low-order finite element method in space. The idea of treating the extension problem via finite elements is already well established for the case of elliptic problems, e.g. [NOS15] for the low-order FEM or [MPSV18] as well as [BMN*+*18] for using $hp$-based discretizations. The use of an extension problem in order to discretize a time-dependent problem was used in [NOS16], focusing on low order finite elements and time-stepping schemes, but allowing also for fractional time derivatives. In the context of wave equations, such a discretization was recently analyzed in [BO19].

When dealing with parabolic problems, it is well-known that, if the initial condition does not satisfy certain compatibility conditions, so called startup singularities form. They need to be accounted for in the numerical method. We rigorously prove that, as long as the meshes are designed in a proper way, our discretization scheme delivers exponential convergence rate for the spatial discretization and optimal convergence rate in time, i.e., optimal order for fixed order timestepping like implicit Euler and exponential convergence for the $hp$-DG based method.

The paper is structured as follows: Section 2 presents the model problem and the functional analytic setting. In Section 3, we then perform a first discretization step with respect to the spatial variables. This yields a continuous in time/discrete in space approximation. In order to prove exponential convergence for this discretization, we take a small detour in Section 3.1 to analyze an auxiliary elliptic problem. This problem will allow us to lift a representation formula from the domain $\Omega$ to the extended cylinder $\Omega\times\mathbb{R}_{+}$ while allowing to reuse the techniques developed in [BMN*+*18]. These preparations then allow us to prove exponential convergence for the space discretization in Section 3.2. The discretization in time is then carried out in Section 4 yielding a fully discrete scheme. This scheme was implemented and Section 5 confirms the exponential convergence. The appendices provide results that could not readily be cited from the literature: Appendix A generalizes results on $hp$-FEM for singularly perturbed problems to the case of complex perturbation parameters. Appendix B is concerned with the lifting of piecewise polynomials in $\Omega$ to piecewise polynomials on the cylinder $\Omega\times\mathbb{R}_{+}$ in a stable way.

We also would like to point out that using the Caffarelli-Silvestre extension is not the only approach to discretize the nonlocal operator which is able to yield an exponentially convergent scheme. We mention schemes based on sinc-quadrature and the Balakrishnan or Riesz-Dunford formulations of the fractional Laplacian (see [BLP17]). We expect that it is possible to combine such a scheme with $hp$-FEM in the space discretization and by combining [BLP17] with the techniques laid out in this paper it should be possible to show exponential convergence.

We close with a remark on notation. We write $A\lesssim B$ to mean there exists a constant $C>0$, which is independent of the main quantities of interest, i.e., mesh size or polynomial degree used, etc., such that $A\leq CB$. We write $A\sim B$ to mean $A\lesssim B$ and $B\lesssim A$. The exact dependencies of the implied constant is specified in the context.

2 Model problem

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded Lipschitz domain. We consider the following model problem for $s\in(0,1)$:

[TABLE]

with initial condition $u(0)=u_{0}$ and right-hand side $f:\Omega\times\mathbb{R}_{+}\to\mathbb{R}$. We assume that the initial condition and right-hand side are analytic but do not require any compatibility or boundary conditions.

The operator $\mathcal{L}u:=-\operatorname{div}(A\nabla u)+cu$ is a linear, elliptic and self-adjoint differential operator, where we assume that A\in L^{\infty}\big{(}\Omega,\mathbb{R}^{d\times d}\big{)} is uniformly SPD in $\Omega$ and $c\in L^{\infty}(\Omega)$ satisfies $c\geq 0$. The fractional power $\mathcal{L}^{s}$ is defined using the spectral decomposition

[TABLE]

where $(\mu_{j},\varphi_{j})_{j\in\mathbb{N}_{0}}$ are eigenvalues and eigenfunctions of the operator $\mathcal{L}$ with homogeneous Dirichlet boundary conditions.

Using the Caffarelli-Silvestre extension one can localize the nonlocal operator $\mathcal{L}^{s}$ and rewrite (2.1) in the following form with $\alpha:=1-2s$:

[TABLE]

Here $\mathcal{C}$ denotes the cylinder $\Omega\times\mathbb{R}_{+}$, $d_{s}:=2^{\alpha}\Gamma(1-s)/\Gamma(s)$. The lateral boundary is defined as $\partial_{L}\mathcal{C}:=\partial\Omega\times\mathbb{R}_{+}$ and

[TABLE]

is the conormal derivative and boundary trace at $y=0$ respectively. The connection to $u$ is then given by $\operatorname{tr}{\mathscr{U}}(t)=u(t)$.

In order to treat this extended problem, we introduce the following weighted Sobolev spaces:

[TABLE]

The space $\mathring{H}^{1}(y^{\alpha},\mathcal{C})$ is equipped with the norm $\displaystyle\left\|\mathscr{U}\right\|_{\mathring{H}^{1}(y^{\alpha},\mathcal{C})}^{2}:=\int_{\mathcal{C}}{y^{\alpha}\left|\nabla\mathscr{U}\right|^{2}}$.

We also define the bilinear form corresponding to the weak form of (2.3a) as:

[TABLE]

Throughout this paper, we will make use of fractional Sobolev and interpolation spaces. We define for two Banach spaces $X_{1}\subseteq X_{0}$ with continuous embedding and $\theta\in(0,1)$:

[TABLE]

For the endpoints $\theta\in\{0,1\}$ we set $[X_{0},X_{1}]_{0,2}:=X_{0}$ and $[X_{0},X_{1}]_{1,2}:=X_{1}$. Fractional Sobolev spaces with and without zero boundary conditions are defined as

[TABLE]

The boundary condition in (2.1) is understood in the sense of $u(t)\in\widetilde{H}^{s}(\Omega)$ for all $t>0$. That is, for $s<1/2$ no boundary condition is imposed, while for $s>1/2$ it is imposed in the sense of traces. For $s=1/2$ the boundary condition is imposed as membership in the Lions-Magenes space, often also denoted $H^{1/2}_{00}(\Omega)$.

Sometimes it is useful to work with a different scale of spaces, characterized using the eigendecomposition of $\mathcal{L}$, as

[TABLE]

For $s\in[0,1]$, the spaces coincide, i.e., $\widetilde{H}^{s}(\Omega)=\mathbb{H}^{s}(\Omega)$ with equivalent norms.

We consider the discretization in two separate steps. We semidiscretize in space and subsequently discretize in time, i.e.,

discretize in space using tensor product $hp$-FEM in $\Omega$ and the artificial variable $y$, 2. 2.

discretize in time by a discontinuous Galerkin method.

3 Discretization in space – the semidiscrete scheme

In this section we investigate the convergence of a semidiscrete semigroup to the solution of (2.1). We consider finite dimensional subspaces $\mathbb{V}_{h}^{\mathcal{X}}\subseteq H_{0}^{1}(\Omega)$ and $\{0\}\neq\mathbb{V}_{h}^{\mathcal{Y}}\subseteq H^{1}(y^{\alpha},\mathbb{R}_{+})$, and set $\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}:=\mathbb{V}_{h}^{\mathcal{X}}\otimes\mathbb{V}_{h}^{\mathcal{Y}}\subseteq\mathring{H}^{1}(y^{\alpha},\mathcal{C})$ as our approximation space. We keep most of our analysis as general as possible, but will provide concrete examples on how to implement these spaces in Sections 3.1.1 and 3.1.2. Throughout the paper, we will write

[TABLE]

While we will give a detailed construction of $\mathbb{V}_{h}^{\mathcal{Y}}$ later on, for now we just assume that there exists $v\in\mathbb{V}_{h}^{\mathcal{Y}}$ with $v(0)=1$ in order to be able to solve Dirichlet problems.

We define the Galerkin approximation $\mathcal{L}^{s}_{h}:\mathbb{V}_{h}^{\mathcal{X}}\to\mathbb{V}_{h}^{\mathcal{X}}$ to the operator $\mathcal{L}^{s}$ via the relation:

[TABLE]

where $\mathscr{D}_{h}:\mathbb{V}_{h}^{\mathcal{X}}\to\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ denotes the solution to the following “lifting problem”:

[TABLE]

We also introduce the notation $\mathscr{D}$ for the solution to

[TABLE]

Remark 3.1**.**

*We note that by [NOS15, Proposition 2.5] and the ellipticity of $\mathcal{A}$, the operator $\mathscr{D}$ is bounded with respect to the $\widetilde{H}^{s}(\Omega)\to\mathring{H}^{1}(y^{\alpha},\mathcal{C})$-norm. For $hp$-FEM spaces, it is non-trivial to show that $\mathscr{D}_{h}$ is bounded, especially on anisotropic meshes. See Appendix B for a related result in a simplified setting.  *

Theorem 3.2**.**

The operator $-\mathcal{L}_{h}^{s}$ is the generator of an analytic semigroup on $\left(\mathbb{V}_{h}^{\mathcal{X}},\left\|\cdot\right\|_{L^{2}(\Omega)}\right)$.

Proof.

The operator $\mathcal{L}_{h}^{s}$ is symmetric due to the symmetry of $\mathcal{A}$. By [Paz83, Section 2.5, Theorem 5.2], it remains to show the estimate

[TABLE]

for $\operatorname{Re}(\lambda)\geq 0$ and a constant $M$ that is independent of $u$ and $\lambda$. It is easy to see that $\left(\lambda\mathrm{I}+\mathcal{L}_{h}^{s}\right)^{-1}f=\operatorname{tr}{\mathscr{U}_{\lambda}}$ where $\mathscr{U}_{\lambda}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ solves

[TABLE]

Existence of the inverse follows from the coercivity of the bilinear form on the left-hand side. The a priori estimate follows by testing with $\mathscr{V}_{h}:=\mathscr{U}_{\lambda}$ to get:

[TABLE]

We use the continuity of the trace operator $\left\|u\right\|_{L^{2}(\Omega)}\lesssim\left\|\mathscr{U}_{\lambda}\right\|_{\mathring{H}^{1}(y^{\alpha},\mathcal{C})}$ and distinguish between the cases $|\lambda|$ large and $|\lambda|$ small to get the desired estimate. ∎

Lemma 3.3**.**

If we equip the space $\mathbb{V}_{h}^{\mathcal{X}}$ with the norm

[TABLE]

the operator $\mathcal{L}_{h}^{s}$ is elliptic, i.e.,

[TABLE]

We also have the following estimate of the $\widetilde{H}^{s}(\Omega)$-norm:

[TABLE]

The constants $c_{i}$ are independent of the spaces $\mathbb{V}_{h}^{\mathcal{X}}$ and $\mathbb{V}_{h}^{\mathcal{Y}}$ and depend only on $\Omega$, $\alpha$, and $\mathcal{L}$.

Proof.

By the trace estimate [NOS15, Proposition 2.5], we get

[TABLE]

On the other hand we get:

[TABLE]

The operator $\mathcal{L}_{h}^{s}$ gives rise to the semidiscrete problem posed in $\mathbb{V}_{h}^{\mathcal{X}}$:

[TABLE]

where $\Pi_{L^{2}}:L^{2}(\Omega)\to\mathbb{V}_{h}^{\mathcal{X}}$ denotes the $L^{2}$-orthogonal projection and $u_{h,0}\in\mathbb{V}_{h}^{\mathcal{X}}$ denotes some approximation to the initial condition.

By Duhamel’s principle, $u$ and $u_{h}$ can be written as

[TABLE]

where $\mathcal{E}:\mathbb{R}_{+}\to\mathscr{B}(L^{2}(\Omega),L^{2}(\Omega))$ and $\mathcal{E}_{h}:\mathbb{R}_{+}\to\mathscr{B}(\mathbb{V}_{h}^{\mathcal{X}},\mathbb{V}_{h}^{\mathcal{X}})$ are the semigroups generated by $-\mathcal{L}^{s}$ and $-\mathcal{L}^{s}_{h}$ respectively.

When considering the discrete flow for initial conditions without compatibility conditions, the right spaces will be the following:

Definition 3.4**.**

Let $\beta\in(0,1)$. Recall that the space $\mathbb{V}_{h}^{\mathcal{X}}$ is equipped with the norm $\left\|u\right\|_{\mathbb{V}_{h}^{\mathcal{X}}}:=\left\|\mathscr{D}_{h}u\right\|_{\mathring{H}^{1}(y^{\alpha},\mathcal{C})}$. We define the interpolation spaces

[TABLE]

We employ the convention $\left\|\cdot\right\|_{\mathbb{V}^{\mathcal{X}}_{h,0}}=\left\|\cdot\right\|_{L^{2}(\Omega)}$ and $\left\|\cdot\right\|_{\mathbb{V}^{\mathcal{X}}_{h,1}}=\left\|\cdot\right\|_{\mathbb{V}_{h}^{\mathcal{X}}}$ for the endpoints.

Throughout this paper, we will work with abstract spaces $\mathbb{V}_{h}^{\mathcal{X}}$. Exponential convergence of the numerical method relies on the following Assumptions 3.5 and 3.9:

Assumption 3.5**.**

There exist constants $\beta$, $b$, $\mu>0$, such that for all $u_{0}$ that are analytic on a fixed neighborhood $\widetilde{\Omega}$ of $\overline{\Omega}$, there exists a function $u_{h,0}\in\mathbb{V}_{h}^{\mathcal{X}}$ and constants $C_{stab},C_{approx}>0$ such that

[TABLE]

where \mathcal{N}_{\Omega}:=\operatorname{dim}\big{(}\mathbb{V}_{h}^{\mathcal{X}}\big{)}.

When considering the Riesz-Dunford representation of $u$, the contour lies in the set of values for which $\mathcal{L}-z$ is elliptic. Therefore we consider the set of complex numbers up to a cone containing the part of the positive real axis for which $\mathcal{L}-z$ is no longer elliptic.

Definition 3.6**.**

With the Poincaré constant $C_{P}$ of $\Omega$ and fixed $0<\varepsilon_{0}<z_{0}\leq\min\left(\frac{1}{2C_{\text{P}}},1\right)^{2}$, we define

[TABLE]

Remark 3.7**.**

*The set $\mathscr{S}$ is chosen in such a way that it contains the contour $\mathcal{C}$ used in the proof of Theorem 3.23. Namely, it contains the rays $\{re^{i\pi/4}\,|\,r>0\}$, $\{re^{-i\pi/4}\,|\,r>0\}$ as well as the circular arc $\{r_{0}e^{i\theta}\,|\,\theta\in(-\pi/4,\pi/4)\}$, with $\varepsilon_{0}<r_{0}<z_{0}$ connecting the two rays. The ball $B_{\varepsilon_{0}}(0)$ is removed in order to avoid problems at [math] when dividing by $z$. See Figure 3.1.  *

Definition 3.8**.**

A function $f:[0,T]\to L^{\infty}(\Omega)$ is said to be uniformly analytic if:

1. (i)

For all $t\in[0,T]$, $f(t)$ is analytic in a fixed neighborhood $\widetilde{\Omega}$ of $\overline{\Omega}$, 2. (ii)

there exist constants $C_{f},\gamma_{f}>0$, the analyticity constants of $f$, such that for all $t\in[0,T]$ and $p\in\mathbb{N}_{0}$,

[TABLE]

The second assumption we have to make is that for a certain class of singularly perturbed elliptic problems, the solution can be approximated exponentially well. We formalize this as follows.

Assumption 3.9**.**

A function space $\mathbb{V}_{h}^{\mathcal{X}}$ is said to resolve down to the scale $\varepsilon>0$ if for all $z\in\mathscr{S}$ with $\left|z\right|^{-1/2}\geq\varepsilon$ and for all functions $f$ that are analytic on a fixed neighborhood $\widetilde{\Omega}$ of $\overline{\Omega}$, the solutions to the elliptic problem

[TABLE]

can be approximated exponentially well from it. That is, there exist constants $C(f)$, $b$ and $\mu>0$ such that

[TABLE]

where $\mathcal{N}_{\Omega}:=\operatorname{dim}(\mathbb{V}_{h}^{\mathcal{X}})$. The constant $C(f)$ may depend only on $\widetilde{\Omega}$, the analyticity constants of $f$, on $A$, $c$, $\Omega$, $z_{0}$ and $\varepsilon_{0}$, while the constants $b$ and $\mu$ depend only on $A$, $c$, $\widetilde{\Omega}$, $\Omega$, $z_{0}$ and $\varepsilon_{0}$ Most notably the constants are independent of $z$, $\varepsilon$, and $\mathcal{N}_{\Omega}$.

For simplicity of notation, we assume that the constants $b$ and $\mu$ in Assumptions 3.5 and 3.9 coincide. All our results will hold for general spaces $\mathbb{V}_{h}^{\mathcal{X}}$, as long as they resolve specific scales. We will later provide a concrete example of constructing such spaces in $1D$ and 2D, see also [BMN*+*18].

The next lemma collects some facts about the time evolution. These results are well-known for the case of the heat equation, and their proof easily carries over to our setting.

Lemma 3.10**.**

The following statements hold for the continuous and the semidiscrete problems:

1. (i)

The maps $t\mapsto u(t)$ and $t\mapsto u_{h}(t)$ are in C\big{(}[0,\infty),L^{2}(\Omega)\big{)}. 2. (ii)

For all $t>0$ and $\ell\in\mathbb{N}_{0}$, $\beta\in(0,1)$ and $\gamma\in[0,1]$ such that $2\ell+\gamma-\beta\geq 0$,

[TABLE]

provided that the right-hand side is finite. In the discrete setting, these estimates read as

[TABLE]

provided that the right-hand side is finite. 3. (iii)

Set $w_{h}:=\int_{0}^{t}{\mathcal{E}_{h}(\tau)\Pi_{L^{2}}f(t-\tau)\,d\tau}$. Then the following estimates hold:

[TABLE]

Proof.

Statement (i) is one of the defining properties of a $C_{0}$-semigroup; thus it follows from Theorem 3.2. The statements of (ii) follow by considering the eigen decomposition of $u$, see [Tho06, Lemma 3.2]. The proof of (iii) is an energy argument. ∎

Corresponding to the operator $\mathcal{L}^{s}_{h}$, we define the Ritz approximation $\Pi_{h}:\operatorname{dom}(\mathcal{L}^{s})\to\mathbb{V}_{h}^{\mathcal{X}}$ via

[TABLE]

(Note: unlike in the heat equation case, the operator $\Pi_{h}$ is not a projection). Since the bilinear form on the left-hand side is elliptic by Lemma 3.3 and $\left(\mathcal{L}^{s}u,v\right)_{L^{2}(\Omega)}$ is a linear functional in $v$, $\Pi_{h}u$ exists and is well defined. (Since $\mathbb{V}_{h}^{\mathcal{X}}$ is finite dimensional we do not have to worry about the norms involved.)

Lemma 3.11**.**

Let $u$ solve (2.1), and $u_{h}$ solve (3.5). Define $\rho:=u-\Pi_{h}u$ and $\theta:=\Pi_{h}u-u_{h}$. Then $\theta$ satisfies the following semidiscrete equation for all $t>0$:

[TABLE]

Proof.

Straightforward computation, see [Tho06, Equation (1.27)]. ∎

The following proposition holds:

Proposition 3.12**.**

Let $u$ solve (2.1), and $u_{h}$ solve (3.5). Define $\rho:=u-\Pi_{h}u$ and $\theta:=\Pi_{h}u-u_{h}$. Then the following estimates hold for all $t>0$:

[TABLE]

Proof.

These estimates are well known for the case of the heat equation. Similar results and techniques can be found, for example, in [Tho06, Chapter 3]. The use of the backward parabolic problem goes back at least to [LR82]. For completeness, we provide a proof in Appendix C. ∎

The previous results mean that it is sufficient to analyze the behavior of the Ritz approximation when applied to $u$. We start this endeavor by showing that the Ritz approximation is quasi-optimal.

Lemma 3.13**.**

Let $u\in\operatorname{dom}(\mathcal{L}^{s})$, and let $\mathscr{D}u$ denotes its lifting to $\mathring{H}^{1}(y^{\alpha},\mathcal{C})$ defined in (3.3). Then the following estimate holds:

[TABLE]

Proof.

We set $u_{h}:=\Pi_{h}u$, and show Galerkin orthogonality $\mathcal{A}(\mathscr{D}u-\mathscr{D}_{h}u_{h},\mathscr{V}_{h})=0$ for all $\mathscr{V}_{h}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$. We first note that $\mathcal{A}(\mathscr{D}u,\mathscr{V}_{h})$ and $\mathcal{A}(\mathscr{D}_{h}u_{h},\mathscr{V}_{h})$ depend only on the trace of $\mathscr{V}_{h}$: By the definition of the liftings (see (3.2a) and (3.3a) respectively), we have for $\mathscr{W}_{h}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ with $\operatorname{tr}{\mathscr{W}_{h}}=\operatorname{tr}{\mathscr{V}_{h}}$:

[TABLE]

Therefore, we get by inserting the definition of $u_{h}=\Pi_{h}u$ and (3.6):

[TABLE]

since $\operatorname{tr}{\left(\mathscr{D}_{h}\operatorname{tr}{\mathscr{V}_{h}}-\mathscr{D}\operatorname{tr}{\mathscr{V}_{h}}\right)}=0$ and (3.3) holds. The approximation result then follows easily from the boundedness of the trace operator and the ellipticity of $\mathcal{A}$. ∎

The combination of Proposition 3.12 and Lemma 3.13 shows that we need to study the best approximation of $\mathscr{U}(t)=\mathscr{D}[u(t)]$ in the space $\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$. This will be done in the next sections.

3.1 A related elliptic problem

In this section, we analyze a family of elliptic problems that will allow us to pass from the function $u\in\widetilde{H}^{s}(\Omega)$ to $\mathscr{U}\in\mathring{H}^{1}(y^{\alpha},\mathcal{C})$.

Instead of using the more intuitive lifting $\mathscr{D}_{h}$, we use one in the form of a Neumann problem. This is done so as to be able to reuse the techniques developed in [BMN*+*18] instead of having to analyze a Dirichlet problem from scratch.

Definition 3.14**.**

Let $\lambda>0$ be fixed. For $f\in L^{2}(\Omega)$, we define the solution operator $\mathcal{G}^{\lambda}f$ by:

[TABLE]

Lemma 3.15**.**

The following stability estimate holds:

[TABLE]

The implied constant depends only on $c$, $A$, and $\Omega$ but is independent of $\lambda$ and $f$.

Proof.

We note that

[TABLE]

Inserting the definition of $\mathcal{G}^{\lambda}$ gives:

[TABLE]

Remark 3.16**.**

*This “damping property” of the factor $\lambda^{-1/2}$ in (3.13) is the main motivation for considering such operators, compared to the more intuitive $\lambda=0$ case, which is the operator analyzed in [BMN*+*18]. It will allow us to better control the behavior of $\mathscr{U}$ for small times $t$ by choosing $\lambda\sim 1/t$, see Section 3.2. It is also the operator which needs to be inverted when discretizing using a implicit Euler timestepping scheme, where $\lambda^{-1}$ is the timestep size, see Section 4. We also point out the strong relation of the operator $\mathcal{G}^{\lambda}$ to the resolvent $(\lambda+\mathcal{L})^{-1}$, see the proof of Theorem 3.2.  *

3.1.1 Discretization of the extended variable $y$

$hp$-fem in 1d:

In this section, we introduce the basics of $hp$-Finite Elements in 1D. This will provide us with the discretization scheme for the extended variable $y$. Additionally, it will serve as a model construction for satisfying Assumptions 3.5 and 3.9.

We introduce the notion of a geometrically refined mesh. For a grading factor $0<\sigma<1$ and $L\in\mathbb{N}$ layers, the geometric mesh on the domain $(-1,1)$ refined towards $-1$, denoted by $\mathcal{T}^{L}_{(\mathbf{-1},1)}:=(x_{i})_{i=0}^{L+1}$ is given by

[TABLE]

Analogously we define the geometric mesh refined towards $1$ and denote it by $\mathcal{T}^{L}_{(-1,\mathbf{1})}$, and the mesh geometrically refined towards both endpoints $\mathcal{T}_{(\mathbf{-1},\mathbf{1})}^{L}$ with nodes at

[TABLE]

In general, triangulations on $(a,b)$, for example denoted by $\mathcal{T}_{(\mathbf{a},b)}^{L}$ are obtained by an affine mapping of $\mathcal{T}_{(\mathbf{-1},1)}^{L}$ etc.

Let $\mathcal{T}$ be a triangulation of a domain $\Omega$. For a polynomial degree distribution $\mathbf{r}\in\mathbb{N}_{0}^{\left|\mathcal{T}\right|}$, we define the space of piecewise polynomials

[TABLE]

For the discontinuous case, we define:

[TABLE]

To simplify the notation, we write $\SS^{p,1}(\mathcal{T}):=\SS^{(p,\dots,p),1}(\mathcal{T})$ for the case of constant polynomial degree $p$, and analogously for $\SS^{p,0}(\mathcal{T})$.

We also sometimes need to impose Dirichlet conditions on the boundary. We write

[TABLE]

The space $\mathbb{V}_{h}^{\mathcal{Y}}$:

We now give the precise construction for the space $\mathbb{V}_{h}^{\mathcal{Y}}$. It is based on an $hp$-FEM on a graded mesh. The details are laid out in the next definition.

Definition 3.17**.**

Fix $\mathcal{Y}>0$. Let $\mathcal{T}_{(\mathbf{0},\mathcal{Y})}^{L}$ be a geometric mesh on $(0,\mathcal{Y})$, refined towards [math] with $L$ layers and a grading factor $\sigma\in(0,1)$, i.e., given by the nodes $\{0,\mathcal{Y}\,\sigma^{i}\,|\,i=0,\ldots,L\}$. Assume that $\mathcal{Y}\sim L$. Let $\mathbb{V}_{h}^{\mathcal{Y}}:=\SS^{\mathbf{r},1}(\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})})\cap\{u:u(\mathcal{Y})=0\}$ be the space of piecewise polynomials with degree distribution vector $\mathbf{r}$ which vanish at the endpoint $\mathcal{Y}$.

Using the eigenpairs $(\varphi_{j},\mu_{j})_{j=0}^{\infty}$ from (2.2), we have the following representation of $\mathscr{U}:=\mathcal{G}^{\lambda}f$:

[TABLE]

Here $\psi_{j}$ are the functions from [BMN*+*18, Formula (4.2)]. They satisfy the differential equation:

[TABLE]

Lemma 3.18**.**

The coefficients $u_{j}$ satisfy the follwing a priori estimate:

[TABLE]

Proof.

From the definition, we get by multiplying with $u_{j}$:

[TABLE]

which implies $\lambda\left|u_{j}\right|\leq|\!\left(f,\varphi_{j}\right)_{L^{2}(\Omega)}\!|$. Inserting this knowledge gives:

[TABLE]

Lemma 3.19**.**

Let $\Pi_{\mathcal{Y}}$ denote the Galerkin projection onto the space $H_{0}^{1}(\Omega)\otimes\mathbb{V}_{h}^{\mathcal{Y}}$ for the problem (3.12). Then the following estimate holds for all $f\in L^{2}(\Omega)$:

[TABLE]

Proof.

We follow the argument of [BMN*+*18, Secs. 4, 5]. By Galerkin orthogonality, we are only concerned with proving an estimate for the best approximation to $\mathcal{G}^{\lambda}f$. The functions $\psi_{j}$ all decay exponentially for $y\to\infty$. We can bound

[TABLE]

where $\lambda_{1}>0$ denotes the smallest eigenvalue of the operator $\mathcal{L}$ on $\Omega$, see [NOS15, Lemma 3.3] for details; the proof can be taken verbatim, just replacing the definition of the coefficients $u_{j}$. It is thus sufficient to study the approximation on the finite cylinder $\Omega\times(0,\mathcal{Y})$.

We define the weights $\omega_{\beta,\gamma}:=y^{\beta}e^{\gamma y}$, and the weighted $L^{2}$-norms

[TABLE]

We note that the function $\mathcal{G}^{\lambda}u$ satisfies the following a priori estimates:

[TABLE]

Again, this follows [BMN*+*18, Theorem 4.7] verbatim, only plugging in the stronger estimate for the coefficients $u_{j}$ to gain the extra factor $\lambda^{-1/2}$. This in turn implies that $\mathcal{G}^{\lambda}f$ is in some Banach-space valued countably normed spaces. Invoking the interpolation operator $\Pi_{y,\{\mathcal{Y}\}}^{\mathbf{r}}$ from [BMN*+*18, Section 5.5.1] then shows the stated result. ∎

3.1.2 Discretization in $x$

In this section, we study the discretization error due to the choice of space $\mathbb{V}_{h}^{\mathcal{X}}$. We will show that the requirement that $\mathbb{V}_{h}^{\mathcal{X}}$ resolves appropriate scales (see Assumption 3.9) suffices to show exponential convergence.

Before we prove an approximation result for $\mathcal{G}^{\lambda}$, we need the following result on the solution of singularly perturbed problems, generalizing the theory developed in, e.g., [Mel97, Mel02] (for real singular perturbation parameters) to the case where the right hand side is itself the solution to a singularly perturbed problem:

Lemma 3.20**.**

Let $\varepsilon>0$ and $z\in\mathscr{S}$ with $\operatorname{Re}(z)\geq 0$. Assume that the space $\mathbb{V}_{h}^{\mathcal{X}}$ resolves the scales $\varepsilon$ and $\left|z\right|^{-1/2}$, as defined in Assumption 3.9. Let $u_{z}\in H_{0}^{1}(\Omega)$ denote the solution to $(\mathcal{L}-z)u_{z}=z\,f$, where $f\in L^{2}(\Omega)$ is analytic on $\overline{\Omega}$. Let $u\in H_{0}^{1}(\Omega)$ solve

[TABLE]

Then the following best approximation result holds:

[TABLE]

The constant $C$ depends on $\mathscr{S}$, the constants of analyticity of $f$, and the constants from Assumption 3.9 but not on $\varepsilon$ or $z$. $b$ and $\mu$ are as in Assumption 3.9.

Proof.

We make the ansatz $u=\nu u_{z}-w$, for $\nu\in\mathbb{C}$ and some function $w\in H_{0}^{1}(\Omega)$. Plugging this decomposition into (3.14) and using the PDE for $u_{z}$, we get $\nu=\frac{1}{1+\varepsilon^{2}z}$ and that $w$ solves

[TABLE]

Since we assumed $\operatorname{Re}(z)\geq 0$, the coefficient $\nu$ is bounded independently of $\varepsilon$ and $z$. We also compute

[TABLE]

which shows that $\frac{\varepsilon^{2}z}{1+\varepsilon^{2}z}$ is also uniformly bounded.

Since we assumed that the mesh resolves the scale $\varepsilon$, we can apply Assumption 3.9 to $w$ to get the estimate:

[TABLE]

We also assumed that the mesh resolves the scale $\left|z\right|^{-1/2}$. Thus we get an exponential approximation property for $u_{z}$ in the $\left|z\right|^{-1/2}$ weighted norm. In order to get the estimate in the $\varepsilon$-weighted norm, we note that for $\varepsilon<\left|z\right|^{-1/2}$ we get the estimate trivially. For $\varepsilon>\left|z\right|^{-1/2}$ we note that

[TABLE]

This means we can approximate $\nu\,u_{z}$ in the $\varepsilon$-weighted norm at an exponential rate, which concludes the proof. ∎

We now employ the decoupling strategy of [BMN*+*18]. Let $(v_{i})_{i=0}^{\mathcal{M}}\subseteq\mathbb{V}_{h}^{\mathcal{Y}}$ denote a basis with the following properties:

[TABLE]

for coefficients $\kappa_{i}>0$. Since the bilinear forms are SPD, such a basis exists. On $\Omega$, we define the bilinear forms

[TABLE]

and note that the following norm equivalence holds on $H_{0}^{1}(\Omega)\otimes\mathbb{V}_{h}^{\mathcal{Y}}$ for all $\mathscr{V}:=\sum_{i=0}^{\mathcal{M}}{\mathscr{V}_{i}v_{i}}$:

[TABLE]

(3.17) shows that estimates in the $\mathring{H}^{1}(y^{\alpha},\mathcal{C})$ norm can also be obtained from bounds on each component in the corresponding $\kappa_{i}$-weighted $H^{1}$-norm.

The bilinear forms $a_{\kappa_{i}}$ correspond to singularly perturbed problems for small $\kappa_{i}$. We want to apply Assumption 3.9. For this we need bounds on the $\kappa_{i}$ as well as on $v_{i}(0)$.

Lemma 3.21**.**

Let $h_{\operatorname{min}}>0$ denote the smallest element size in $\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})}$ and $p$ the maximal polynomial degree used for $\mathbb{V}_{h}^{\mathcal{Y}}$. Then the eigenpairs $(v_{i},\kappa_{i})_{i=0}^{\mathcal{M}}$ of (3.15) satisfy:

[TABLE]

Proof.

By definition we have $1=d_{s}\lambda v_{i}(0)^{2}+\int_{0}^{\mathcal{Y}}{y^{\alpha}|v_{i}^{\prime}|^{2}}=\kappa_{i}^{-1}\int_{0}^{\mathcal{Y}}{y^{\alpha}|v_{i}|^{2}}$, or $\kappa_{i}=\int_{0}^{\mathcal{Y}}{y^{\alpha}|v_{i}|^{2}}$. By the proof of [BMN*+*18, Lemma B.2] we can estimate

[TABLE]

[BMN*+*18, Lemma B.1] provides $|v_{i}(0)|^{2}\leq\mathcal{Y}^{1-\alpha}/(1-\alpha^{2})\left\|v_{i}^{\prime}\right\|^{2}_{L^{2}(y^{\alpha},(0,\mathcal{Y}))}$. This and the inverse estimate from [BMN*+*18, Lemma B.3] yield

[TABLE]

To see (3.19), we calculate:

[TABLE]

Lemma 3.22**.**

Let $u\in L^{2}(\Omega)$ be either holomorphic in $\overline{\Omega}$ or the solution to the singularly perturbed problem $-z^{-1}\mathcal{L}u+u=f$ with $f$ holomorphic on $\overline{\Omega}$ and $z\in\mathscr{S}$ with $\operatorname{Re}(z)\geq 0$. Assume that $\mathbb{V}_{h}^{\mathcal{X}}$ resolves the scales ${\left|z\right|}^{-1/2}$ and $\sqrt{\kappa_{i}}$ for all $i=0,\dots,\mathcal{M}$.

Then the following best approximation result holds:

[TABLE]

where $\mu$ is the exponent for $\mathbb{V}_{h}^{\mathcal{X}}$ in Assumption 3.33. The constant $b$ depends on the domain of holomorphy of $u$ or $f$.

Proof.

By Lemma 3.19, it is sufficient to consider a semidiscrete functions $\mathscr{U}^{h}_{y}:=\Pi_{\mathcal{Y}}\mathcal{G}^{\lambda}u\in H_{0}^{1}(\Omega)\otimes\mathbb{V}_{h}^{\mathcal{Y}}$ and their approximation in $\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$. Using the basis $(v_{j})_{j=0}^{\mathcal{M}}$, the function $\mathscr{U}^{h}_{y}=:\sum_{i=0}^{\mathcal{M}}{U_{i}v_{i}}$ from Lemma 3.19 solves:

[TABLE]

This is just the weak formulation of the singularly perturbed problems from Assumption 3.9, with $\varepsilon=\sqrt{\kappa_{i}}$. Note that $|v_{i}(0)|\lesssim\lambda^{-1/2}$ by (3.19). Since we assumed that the scales are resolved, we can either apply Assumption 3.9 directly (if $u$ is holomorphic on $\overline{\Omega}$) or apply Lemma 3.20 (if $u$ solves $-z^{-1}\mathcal{L}u+u=f$) to get the following estimate for the best approximations $\Pi U_{i}\in\mathbb{V}_{h}^{\mathcal{X}}$:

[TABLE]

the norm equivalence (3.17) then concludes the proof. ∎

3.2 Returning to the semidiscretization

We are now in a position to show exponential convergence for the best approximation (and thus also the Ritz approximation) of the exact solution $\mathscr{U}$. We first consider positive times $t$ bounded away from [math]. In this regime, our finite element mesh is assumed to resolve the pertinent scales. The smaller times, for which the scales are not resolved, are treated separately later on.

Theorem 3.23**.**

Let $t\geq t_{0}>0$ be fixed. Let $u_{0}$ be analytic on a fixed neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$ (but we do not assume boundary conditions, i.e., $u_{0}\notin\widetilde{H}^{s}(\Omega)$ is allowed), and assume homogeneous right-hand side, i.e., $f=0$. Also assume that for a chosen “high frequency” cutoff $z_{\text{hf}}>z_{0}>0$, the space $\mathbb{V}_{h}^{\mathcal{X}}$ resolves down to the scale

[TABLE]

where $h_{\min}$ and $p$ are the minimum element size and maximal polynomial degree of $\mathbb{V}_{h}^{\mathcal{Y}}$. Then, for each $\ell\in\mathbb{N}_{0}$, there exists a function $\mathscr{V}_{h}(t)\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ such that the following estimate holds:

[TABLE]

The implied constant depends on $\Omega$,$\widetilde{\Omega}$, $s$, the constants of analyticity of $u_{0}$, $z_{0}$, and the constants from Assumption 3.9, but is independent of $t$, $t_{0}$ and $z_{\text{hf}}$. The rate $b_{2}$ also depends on the mesh grading for $y$. The rate $b_{1}$ depends in addition on the constants from Assumption 3.9. The constant $\gamma$ can be chosen to depend on $s$ only.

Proof.

Since we assumed homogeneous right hand side, we only need to investigate $\mathscr{U}=\mathcal{E}(t)u_{0}$. We use the representation of $\mathcal{E}(t)u_{0}$ via the Riesz-Dunford calculus (following what is done in [BLP17, Section 2]) to write:

[TABLE]

where $\mathcal{C}$ is the following contour consisting of three segments:

[TABLE]

and $z^{s}:=e^{s\log(z)}$ with the logarithm defined with the branch cut along the negative real axis. The parameter $r_{0}\in(0,z_{0})$ is fixed such that the whole path lies in the domain of ellipticity $\mathscr{S}$, as defined in Definition 3.6; see Figure 3.1.

By adding the term $\frac{1}{t}d_{s}\operatorname{tr}{\mathscr{U}}$ to both sides of (2.3b), we get that $\mathscr{U}$ solves

[TABLE]

Using the operator $\mathcal{G}^{1/t}$, we can therefore write the function $\mathscr{U}$ as

[TABLE]

or using the Riesz-Dunford calculus:

[TABLE]

For the derivatives, a similar formula holds:

[TABLE]

Hence, we have to study integrals of the form

[TABLE]

and their best approximation, paying attention to the dependence on $t$.

If $\left|z\right|\in(\varepsilon_{0},z_{\text{hf}})$, we obtain from Lemma 3.22 (since we require Assumption 3.9 to hold down to scale $z_{\text{hf}}^{-1/2}$)

[TABLE]

for some function $\widehat{\mathscr{V}}_{h}(z)\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$. If we pick $\widehat{\mathscr{V}}_{h}(z)$ as the Galerkin approximation, we get continuous dependence on $z$. On $\mathcal{C}_{2}$, we can therefore estimate:

[TABLE]

The more interesting case are the paths $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$. We focus on $\mathcal{C}_{1}$, and consider two cases, namely, $\left|z\right|\leq z_{\text{hf}}$ and $\left|z\right|>z_{\text{hf}}$. In the first case, the mesh resolves the scales down to $z_{\text{hf}}^{-1/2}\leq\left|z\right|^{-1/2}$, and we can apply Lemma 3.22. Setting $\gamma:=\cos(\pi\,s/4)$ we estimate:

[TABLE]

Making the substitution $\gamma t\,r^{s}=:y$, we get:

[TABLE]

We need to consider the case $m=0$ separately, as the integrand then has a singularity at $r=0$. Splitting the integration we get:

[TABLE]

For $m>0$, we do not get the logarithmic growth for small times, since:

[TABLE]

Overall, this gives the estimate:

[TABLE]

In the case $\left|z\right|>z_{\text{hf}}$, we set $\widehat{\mathscr{V}}_{h}:=0$ and use the stability estimate (3.13) together with the uniform stability of the operator $(z-\mathcal{L})^{-1}z$ (see Lemma A.2). For $m>0$, we estimate:

[TABLE]

For $m=0$, the same calculation can be done, but picking up an extra logarithmic term from the integral where $y=z_{hf}^{s}t\lesssim 1$.

The same argument can be repeated for $\mathcal{C}_{3}$. The stated estimates then follow easily by setting $m=0$ and $m=1$ to estimate $\mathscr{U}$ (this term involves the logarithmic contributions) and $m=\ell$ and $m=\ell+1$ to estimate higher derivatives. ∎

For small $t<t_{0}$, we cannot hope to retain exponential convergence, as it would require our mesh to resolve infinitely small scales. Instead, we rely on on our ability to control the behavior of the solution near $t=0$ using some smoothness of $u_{0}$.

Lemma 3.24**.**

Let $u_{0}\in H^{\theta}(\Omega)$ for $0<\theta<1/2$, and assume homogeneous right hand-side, i.e., $f=0$. For all $\ell\in\mathbb{N}_{0}$, the following estimate holds for $t>0$:

[TABLE]

The constant depends on $\Omega$, $\theta$, $s$ and the coefficients $A$, $c$.

Proof.

For simplicity we assume additionally $\theta\leq 2s$. We note that for $\theta\in(0,1/2)$, the spaces $\widetilde{H}^{\theta}(\Omega)$ and $H^{\theta}(\Omega)$ coincide with equivalent norms (see [Tri06, Section 1.11.6] or [McL00, Theorem 3.33, Theorem B.9, Theorem 3.40]).

Hence, we get $u_{0}\in\widetilde{H}^{\theta}(\Omega)$. By Lemma 3.10, this implies for $\ell\in\mathbb{N}_{0}$:

[TABLE]

We write $\mathscr{U}(t)=\mathscr{D}u(t)$ using the lifting operator from (3.3). Since the lifting on the continuous level is bounded (see Remark 3.1), we can estimate:

[TABLE]

As a final step before showing convergence of the semidiscrete approximation, we remove the restriction to homogeneous right-hand sides $f$. This is a simple consequence of the previous results and Duhamel’s principle.

Corollary 3.25**.**

Let $t_{0}>0$ and $\delta>0$ be fixed. Let $u_{0}$ be analytic on $\overline{\Omega}$ and assume that $f$ is $\ell$ times continuously differentiable with respect to $t$ such that the functions $f^{(j)}$, $j=0,\ldots,\ell$ are uniformly analytic in the sense of Definition 3.8.

Assume that $\mathbb{V}_{h}^{\mathcal{X}}$ resolves scales down to (3.20). Then, for each $\ell\in\mathbb{N}_{0}$, there exists a function $\mathscr{V}_{h}(t)\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ such that the following estimates holds for all $t\in(0,T)$:

[TABLE]

The implied constant depends on the end time $T$, $\Omega$, the data $u_{0}$, the constants of analyticity of $f^{(j)}$, $\delta$, and the implied constants in Lemma 3.22, e.g., the mesh grading factor. It is independent of $t$, $t_{0}$,$z_{\text{hf}}$, $\mathcal{N}_{\Omega}$ or $\mathcal{N}_{\mathcal{Y}}$. For $\ell=0$ and $\ell=1$ we can explicitly give $C(T)\lesssim\max(1,T)$.

Proof.

For $f=0$, this is just a collection of Lemma 3.23 and 3.24. For $f\neq 0$ we write

[TABLE]

(see [Paz83, Section 4.2, Corollary 2.5] for the derivative of Duhamel’s formula). The terms involving only $\mathcal{E}(t)$ are already covered by the results for the homogeneous problem. For fixed $\tau\in(0,t)$, the integrand in the last term corresponds to solving the homogeneous problem with initial condition $f(t-\tau)$ (or $\dot{f}(t-\tau)$ in the case of $\dot{\mathscr{U}}$). This means we can also apply Lemmas 3.23 and 3.24, only picking up an extra power of $t$ due to the additional integration in $\tau$. This gives the stated estimate for $\ell=0$ and $\ell=1$.

For higher derivatives, we proceed by induction and see that we can write $\mathscr{U}^{(\ell)}$ as

[TABLE]

All the terms can be estimated as before, where we estimate $t^{-j}\leq C(T)t^{-\ell}$ and only keep the dominant terms. ∎

Theorem 3.26**.**

Assume that $u_{0}$ is analytic and $f$ is uniformly analytic on a fixed neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$ (in the sense of Definition 3.8). Let $\mathbb{V}_{h}^{\mathcal{Y}}$ be given by Definition 3.17. Fix $t_{0}>0$, $\delta>0$, and set $z_{\text{hf}}:=t_{0}^{-1/s}L^{1/s}$, where $L$ is the number of layers used for constructing the geometric mesh $\mathbb{V}_{h}^{\mathcal{Y}}$. Let the space $\mathbb{V}_{h}^{\mathcal{X}}$ resolve the scales down to (3.20), and let Assumption 3.5 hold for the initial condition. Then the following estimate holds:

[TABLE]

Proof.

We just collect all the previous results, most notably Proposition 3.12 and Corollary 3.25. Since we only need the best approximation estimate on $\mathscr{U}$ and $\dot{\mathscr{U}}$, we keep the dependence on the time $t$ explicit. The error due to the different initial conditions is exponentially small by assumption. ∎

We can also obtain estimates in the energy norm or pointwise in time:

Theorem 3.27**.**

Assume that $u_{0}$ is analytic, $f$ and $\dot{f}$ are uniformly analytic on a neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$, and that $u_{h,0}\in\mathbb{V}^{\mathcal{X}}_{h,\beta}$ is as in Assumption 3.5.

Let $L$ denote the number of layers used for $\mathbb{V}_{h}^{\mathcal{Y}}$, set $t_{0}:=e^{-L}$, and $z_{\text{hf}}:=t_{0}^{-1/s}L^{1/s}$ and assume that the space $\mathbb{V}_{h}^{\mathcal{X}}$ resolves the scales down to (3.20).

Set $M:=\min(L,\operatorname{dim}(\mathbb{V}_{h}^{\mathcal{X}})^{\mu})$ with $\mu>0$ from Assumptions 3.5 and 3.9.

Then there exists a constant $b$, independent of $L$, $p$ and the specific choice of $\mathbb{V}_{h}^{\mathcal{X}}$, i.e. depending only on the constants from Assumptions 3.5 and 3.9 such that the following estimate holds:

[TABLE]

Proof.

Without loss of generality, we may assume $\beta s<1/2$. Fix $t_{1}>0$ to be chosen later. We consider two regimes, $t\in(0,t_{1})$ and $t\geq t_{1}$. For $t\leq t_{1}$, we use the stability estimates of Lemma 3.10 (ii), together with the insight that $u_{0}\in\widetilde{H}^{\beta s}(\Omega)$ for $\beta s<1/2$ which was already used in Lemma 3.24.

We start with the energy norm estimate and use Lemma 3.10 to get:

[TABLE]

For the pointwise estimate, we write $u(t)=u_{0}+\int_{0}^{t}{\dot{u}(\tau)d\tau}$ and $u_{h}(t)=u_{h,0}+\int_{0}^{t}{\dot{u}(\tau)d\tau}$ to get:

[TABLE]

For larger times $t>t_{1}$, we can establish the following bound by using (3.11) and plugging in the results on the best approximation from Corollary 3.25.

[TABLE]

Or, since $\tau>t_{1}$:

[TABLE]

Setting $t_{1}\sim e^{-\frac{b^{\prime}}{2}M}$ we get the stated exponential convergence with rate $b:=-b^{\prime}\beta/2$ after using Assumption 3.5 to estimate the error due to approximating the initial condition. ∎

3.3 Example of a space $\mathbb{V}_{h}^{\mathcal{X}}$: $hp$-FEM in 1D and 2D

In this section, we give an exemplary construction for $\mathbb{V}_{h}^{\mathcal{X}}$ given a simpler model problem in one or two space dimension using $hp$-Finite Elements meeting our requirements. In other words, $\mathbb{V}_{h}^{\mathcal{X}}$ satisfies Assumptions 3.5 and 3.9.

Assumption 3.28**.**

The domain $\Omega\subseteq\mathbb{R}^{d}$ for $d=1,2$ has analytic boundary. Also, the coefficient functions $A$ and $c$ are analytic on a neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$.

In 1D, we have already introduced $hp$-FEM spaces. For analytic 2D geometries they are given in the following Definition 3.29. We follow [MS98] and [BMN*+*18]; see also [Mel02, Definition 2.4.1]. We first introduce the (shape regular) reference mesh.

Definition 3.29** (reference mesh).**

Denote by $\widehat{S}:=(0,1)^{2}$ the reference square, and let \mathcal{T}_{\Omega}:=\big{\{}K_{i}\big{\}}_{i=0}^{\left|\mathcal{T}_{\Omega}\right|} be a mesh of curved quadrilaterals with bijective element maps $F_{K}:\overline{\widehat{S}}\to\overline{K}$ satisfying

1. (M1)

The elements $K_{i}$ partition $\Omega$, i.e., $\bigcup_{K_{i}\in\mathcal{T}}\overline{K}=\overline{\Omega}$; 2. (M2)

for $i\neq j$, $\overline{K}_{i}\cap\overline{K}_{j}$ is either empty, a vertex or an entire edge; 3. (M3)

the element maps $F_{K}:\widehat{S}\to K$ are analytic diffeomorphisms; 4. (M4)

the common edge of two neighboring elements $K_{i}$, $K_{j}$ has the same parametrization from both sides, i.e., if $\gamma_{ij}$ is the common edge with endpoints $P_{1}$, $P_{2}$, then for $P\in\gamma_{i,j}$ we have

[TABLE]

Definition 3.30** (anisotropic geometric mesh).**

Given a reference mesh $\mathcal{T}_{\Omega}$. Let $K_{i}$, $i=0,\dots,n<\left|\mathcal{T}_{\Omega}\right|$ be the elements at the boundary, and assume that the left edge $e:=\{0\}\times(0,1)$ is mapped to $\partial\Omega$, i.e., $F_{K_{i}}(e)\subseteq\partial\Omega$ and $F_{K_{i}}(\partial S\setminus\overline{e})=\emptyset$ for $i=0,\dots,n$. Assume that the remaining elements satisfy $\overline{K_{i}}\cap\partial\Omega=\emptyset$, $i=n+1,\dots\left|\mathcal{T}_{\Omega}\right|$.

For $L\in\mathbb{N}$ and a mesh grading factor $\sigma\in(0,1)$, we subdivide the reference square

[TABLE]

The anisotropic geometric mesh $\mathcal{T}_{\Omega}^{L}$ is then given by the push-forwards of the refinements plus the unrefined interior elements:

[TABLE]

Definition 3.31**.**

In one dimension, for $\Omega=(a,b)$, the anisotropic mesh is defined as in Section 3.1.1 by $\mathcal{T}_{\Omega}^{L}:=\mathcal{T}_{(\mathbf{a},\mathbf{b})}^{L}$. The reference mesh is given by the single element \mathcal{T}_{\Omega}:=\big{\{}(a,b)\big{\}}.

We are now able to define the space $\mathbb{V}_{h}^{\mathcal{X}}$ using these meshes.

Definition 3.32** ($\mathbb{V}_{h}^{\mathcal{X}}$ via $hp$-FEM).**

Let $\mathcal{T}_{\Omega}^{L}$ be an anisotropic geometric mesh refined towards $\partial\Omega$ and fix $p\in\mathbb{N}$. We write \mathcal{Q}^{p}:=\operatorname{span}_{0\leq i_{1},\,\dots,i_{d}\leq p}\big{\{}x_{1}^{i_{1}}\,\dots\,x_{d}^{i_{d}}\big{\}} for the space of tensor product polynomials and set

[TABLE]

This choice of approximation space will prove suitable to satisfy Assumptions 3.5 and 3.9. We start with the fact that we can resolve certain scales:

Theorem 3.33**.**

Let $\mathcal{T}^{L}_{\Omega}$ be an anisotropic mesh on $\Omega$ that is geometrically refined towards $\partial\Omega$ with grading factor $\sigma\in(0,1)$ and $L$ layers.

Then $\mathbb{V}_{h}^{\mathcal{X}}$ defined in (3.27) resolves the scales down to $\sigma^{L}$, i.e., there exist constants $C,b>0$, such that for $z\in\mathscr{S}$ with ${\left|z\right|}^{-1/2}>\sigma^{L}$ and every $f$ which is analytic on a neighborhood $\widetilde{\Omega}$ of $\Omega$, the solution $u_{z}$ to $(\mathcal{L}-z)u=zf$ can be approximated by $v_{h}\in\mathbb{V}_{h}^{\mathcal{X}}$ satisfying

[TABLE]

The constant $b$ depends only on $\sigma$ and $\widetilde{\Omega}$. The constant $C$ also depends on the constants of analyticity of $f$.

Proof.

We consider two cases. For $\left|z\right|<2z_{0}$, the problems are not actually singularly perturbed and standard results for $hp$-FEM can be applied. We thus only focus on the case $\left|z\right|>2z_{0}$.

From the definition, we get

[TABLE]

Defining $\varepsilon:=\left|z\right|^{-1/2}$ and $\zeta:=\varepsilon^{2}(c-z)$, the problem can be rewritten as

[TABLE]

We make the following observations:

1. (i)

$\left|\varepsilon^{2}z\right|\sim 1$ and therefore $\|\widetilde{f}\|\sim\|f\|$ for any norm, 2. (ii)

since we assumed $\left|z\right|>2z_{0}$, by slightly decreasing the opening angle, we may ignore the shift in the definition of $\mathcal{S}$ and assume that $\left|\pi-\operatorname{Arg}(\zeta)\right|\geq\delta>0$; see Figure 3.2.

This means that we can apply the results from Appendix A, most notably Theorem A.5. ∎

The $hp$-FEM spaces can also approximate the initial conditions at an exponential rate. But more importantly, they can do so in a way that is stable with respect to the non-standard $\mathbb{V}^{\mathcal{X}}_{h,\beta}$ norm. We start with a simple lemma.

Lemma 3.34**.**

Given $u\in\mathbb{V}_{h}^{\mathcal{X}}$, for any function $\mathscr{V}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ with $\operatorname{tr}{\mathscr{V}}=u$, we can estimate

[TABLE]

In other words, $\mathscr{V}=\mathscr{D}_{h}u$ has “minimal energy”.

Proof.

We compute for $\mathscr{W}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ with $\operatorname{tr}{\mathscr{W}}=0$:

[TABLE]

where we used $\mathcal{A}(\mathscr{D}_{h}u,\mathscr{W})=0$ for $\mathscr{W}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ with $\operatorname{tr}{\mathscr{W}}=0$ by the definition of the lifting. Setting $\mathscr{W}:=\mathscr{D}_{h}u-\mathscr{V}$ then shows the estimate $\left\|u\right\|_{\mathbb{V}_{h}^{\mathcal{X}}}\leq\left\|\mathscr{V}\right\|_{\mathring{H}^{1}(y^{\alpha},\mathcal{C})}$. ∎

Working with anisotropic meshes for the discrete liftings imposes additional difficulties. Instead, we split the lifting process into two steps, first we lift using the shape-regular reference triangulation (but ignoring boundary conditions), and then we use a cutoff procedure to correct the boundary conditions on the anisotropic geometric mesh. This cutoff operator is constructed in the following Lemma.

Lemma 3.35**.**

Let $\mathcal{T}^{L}_{\Omega}$ denote an anisotropic geometric mesh with reference mesh $\mathcal{T}_{\Omega}$. Given $\ell\in\mathbb{N}_{0}$, $\ell\leq L$, there exists a bounded linear operator $\mathscr{C}_{\ell}:\SS^{p,1}(\mathcal{T}_{\Omega})\to\SS^{p,1}_{0}(\mathcal{T}^{L}_{\Omega})$ such that

[TABLE]

and for all $0\leq\beta<1/2$

[TABLE]

Proof.

We fix a layer of thickness $\sigma^{\ell}$ around the boundary and pick a piecewise affine function $\chi$ such that $\chi=1$ on all elements outside of this layer. This can be easily done working on the reference patches. This leads to a function which only has non-vanishing gradient on this layer, and there satisfies the estimate

[TABLE]

Working on the reference element, we note that the Gauss Lobatto in one direction satisfies the following stability estimate, also on anisotropic elements $K$:

[TABLE]

(see [BM97, Rem. 13.5 and (13.27)] for the 1D case. 2D then follows by a tensor product argument, see also [HS03, Eqn. (15)]). Combining such element-wise Gauss-Lobatto interpolants on each element, we get a global stable interpolation operator $\Pi_{p}$.

We then define $\mathscr{C}_{\ell}v:=\Pi_{p}(\chi v)$. We note that since $\chi\,v$ is a piecewise polynomial of degree $p+1$, we get the stability estimates:

[TABLE]

From the estimate on $\nabla\chi$ we then immediately get (3.28).

The approximation estimate follows from the fact that $\Pi_{p}$ reproduces $v$ and is $L^{2}$ stable and the following estimate taken from [LMWZ10, Lemma 2.1]. Since $1-\chi$-vanishes outside of the strip size $\sigma^{\ell}$:

[TABLE]

Lemma 3.36**.**

Assume that the triangulation $\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})}$ used for the discretization in $y$ satisfies $\sigma^{L}\lesssim p_{x}^{-2}$, where $p_{x}$ denotes the (maximal) polynomial degree used for $\mathbb{V}_{h}^{\mathcal{X}}$.

Let $u_{0}$ be analytic in a neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$, and let $0\leq\beta<1/2$ and $\varepsilon>0$. Then there exists a function $u_{h,0}\in\mathbb{V}_{h}^{\mathcal{X}}$ :

[TABLE]

In other words, if the number of refinement layers $L\sim p_{x}$, then $\mathbb{V}_{h}^{\mathcal{X}}$ satisfies Assumption 3.5 with $\mu:=1/(d+1)$.

Proof.

Since $u_{0}$ is analytic, we do not need to approximate any boundary layers or singularities. What we do need to take care of is the fact that our FEM space has homogeneous boundary conditions, while $u_{0}$ does not.

We will construct the lifting in two steps: First, we approximate and lift in a space without boundary conditions and then we will perform a cutoff procedure.

For $\varepsilon\in(0,1/2-\beta)$, Let $\widehat{u}_{0}\in\SS^{p,1}(\mathcal{T}_{\Omega})$ be the $H^{\max(s,\beta+\varepsilon)}(\Omega)$-best approximation of $u_{0}$. (Note that we work on the shape-regular grid $\mathcal{T}_{\Omega}$ and do not impose boundary conditions.) By standard results, we have $\left\|u_{0}-\widehat{u}_{0}\right\|_{L^{2}(\Omega)}\lesssim e^{-b\mathbb{N}_{\Omega}^{\mu}}$ and by Lemma B.4 we can lift this function to $\widehat{\mathscr{U}}_{0}\in H^{1}(y^{\alpha},D)$ such that $\widehat{\mathscr{U}}_{0}(\cdot,y)\in\SS^{p,1}(\mathcal{T}_{\Omega})$ and $\widehat{\mathscr{U}}_{0}(x,\cdot)\in\SS^{1,1}(\mathcal{T}_{\mathcal{Y}})$ for all $y\in(0,\mathcal{Y})$ and $x\in\Omega$.

Using the cutoff operator from Lemma 3.35 we then define $u_{h,0}:=\mathscr{C}_{L}(\widehat{u}_{0})$ and the piecewise constant function $\mathscr{V}(t):\mathbb{R}_{+}\to\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$

[TABLE]

Since $\mathscr{C}_{\ell}$ is $L^{2}$ and $H^{1}$ stable and is applied on each $y$-slice, we get stability in the $H^{1}(y^{\alpha},D)$-norm, i.e. \left\|\mathscr{V}(t)\right\|_{H^{1}(y^{\alpha},D)}\lesssim\sigma^{-\ell/2}\big{\|}{\widehat{\mathscr{U}}_{0}}\big{\|}_{H^{1}(y^{\alpha},D)}. We also get the approximation of the trace at $y=0$ for $t\in(\sigma^{\ell+1},\sigma^{\ell})$:

[TABLE]

We then need to estimate the $\mathbb{V}^{\mathcal{X}}_{h,\beta}$-norm. By construction of the cutoff function, we get that $\mathscr{V}(t)\in\mathbb{V}_{h}^{\mathcal{X}}$ and calculate:

[TABLE]

where we used that we can replace the specific lifting $\mathscr{D}_{h}$ with the function $\mathscr{V}$ as it has the “minimum energy” property via Lemma 3.34.

From the stability estimates on each segment $(\sigma^{\ell+1},\sigma^{\ell})$, we get:

[TABLE]

Where we used $\beta<1/2$, a geometric series, and the stability of the lifting of $\widehat{u}_{0}$ and the best approximation. From the approximation property (3.30) we get:

[TABLE]

The estimate $\int_{1}^{\infty}{t^{-2\beta-1}\left\|u_{h,0}\right\|^{2}_{L^{2}(\Omega)}\,dt}\lesssim\left\|u_{h,0}\right\|_{L^{2}(\Omega)}^{2}\lesssim\left\|\widehat{u}_{0}\right\|_{L^{2}(\Omega)}^{2}$ is trivial.

The approximation estimate from (3.29) then follows from the approximation property of $\mathscr{C}_{L}(\widehat{u}_{0})$ from Lemma 3.35 and the best approximation property of $\widehat{u}_{0}$.

∎

Remark 3.37**.**

*For constructing the lifting, Lemma 3.36 relies on the as of yet unpublished work [MKR19]. In the simpler, one dimensional case, the space $\SS^{p,1}(\mathcal{T}_{\Omega})$ coincides with the space $\mathcal{Q}^{p}$ of global polynomials since the reference mesh only consists of a single element. This allows us to replace [MKR19] with results from [BDM07] in this case.  *

We can now give a more constructive characterization of how the triangulation of $\Omega$ must be chosen when working in 1D or 2D to get exponential convergence of the semidiscretization.

Corollary 3.38**.**

Let $\Omega\subset\mathbb{R}^{d}$, $d=1,2$ have analytic boundary. Assume that $u_{0}$ is analytic and $f$ is uniformly analytic in a neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$. For $M\in\mathbb{N}$ and $\sigma\in(0,1)$, use an anisotropic geometric mesh with $M$ layers to discretize in $x$, i.e., $\mathbb{V}_{h}^{\mathcal{X}}:=\SS^{p,1}_{0}(\mathcal{T}^{M}_{\Omega})$. For discretizing in $y$, use $L$ layers and a degree vector $\mathbf{r}$ with linear slope $\mathfrak{s}$, i.e., $\mathbb{V}_{h}^{\mathcal{Y}}:=\SS^{\mathbf{r},1}(\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})})$. Assume that $\sigma^{M}\leq c_{\mathcal{Y}}\mathcal{Y}(\mathfrak{s}L)^{-2}\sigma^{3L/2}$, $\sigma^{L}<p^{2}$ and $u_{h,0}$ is as in Assumption 3.5.

Then there exist constants $b_{1},b_{2}>0$ independent of $L$, $M$, and $p$ such that the following estimate holds:

[TABLE]

Most notably for $M\sim\frac{3}{2}L$ and $p\sim L$, we get exponential convergence:

[TABLE]

Proof.

We choose $t_{0}:=\sigma^{L}$ and $z_{\text{hf}}=\sigma^{L/s}L^{1/s}$ in Theorem 3.26. Assumption 3.5 is met via Lemma 3.36. The assumptions on $\mathbb{V}_{h}^{\mathcal{X}}$ also imply that the necessary scales are resolved, and we get:

[TABLE]

The explicit estimate then follows from the fact that $\dim(\mathbb{V}_{h}^{\mathcal{X}})\sim L^{d+1}$ in this particular construction and $\dim(\mathbb{V}_{h}^{\mathcal{Y}})\sim L^{2}$. We absorb the logarithmic terms $\log(\sigma^{L})\sim L$ into the exponential by slightly reducing the rate $b^{\prime}$. The condition $\sigma^{L}\leq p^{-2}$ is easily verified for such meshes. ∎

For the pointwise and energy errors, the corresponding concrete version reads:

Corollary 3.39**.**

Assume that $u_{0}$ is analytic and $f$, $\dot{f}$ are uniformly analytic in a neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$, and that the meshes and spaces are as in Corollary 3.38. Let $u_{h,0}\in\mathbb{V}^{\mathcal{X}}_{h,\beta}$ be as in Assumption 3.5 for $\beta>0$.

Then there exists a constant $b$, independent of $L$, $M$ and $p$ such that the following estimate holds:

[TABLE]

Or in terms of degrees of freedom, we get

[TABLE]

Proof.

Follows from the fact that using the given parameters, the space $\mathbb{V}_{h}^{\mathcal{X}}$ satisfies the assumptions of Theorem 3.27. The estimate in terms of degrees of freedom follows easily. ∎

Remark 3.40**.**

*In this section, we focused on the case of smooth geometries in 1D and 2D. We would like to point out that we do not see any structural obstacles towards generalizing to the case of curvilinear polygons or smooth 3d geometries. The main ingredient is the necessary generalization of Appendix A.  *

4 Discretization in $t$ – the fully discrete scheme

In this section, we consider the discretization with respect to the time variable $t$. This can be done using mostly standard techniques. We focus on the case of using a discontinuous Galerkin type method. When applied in its $hp$-version, it will allow us to get an exponentially convergent fully discrete scheme, and thus it nicely complements our previous investigations. We follow the presentation in [SS00].

Let $\mathcal{T}_{(0,T)}:=\{(t_{j-1},t_{j})\}_{j=1}^{M}$ be a partition of the time interval $[0,T]$ into subintervals with $0\leq t_{j}<t_{j+1}\leq T$. We set $k_{j}:=t_{j}-t_{j-1}$ and define the one-sided limits

[TABLE]

as well as the jump $[u]_{j}:=u_{j}^{+}-u_{j}^{-}$. We define the DG-bilinear and linear forms:

[TABLE]

Then the DG-approximation is given as the solution to the following problem:

Problem 4.1**.**

Choose $\mathbf{r}_{t}\subseteq\mathbb{N}_{0}$ a polynomial degree distribution, and consider the space $\SS^{\mathbf{r}_{t},0}(\mathcal{T}_{(0,T)})$ of discontinuous piecewise polynomials. Set $\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y},\mathcal{T}}:=\SS^{\mathbf{r}_{t},0}(\mathcal{T}_{(0,T)})\otimes\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$. Find $\mathscr{U}^{h}_{x,y,t}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y},\mathcal{T}}$ such that

[TABLE]

Remark 4.2**.**

*Note that we used the initial condition $u_{0}$ instead of the discrete initial condition $u_{h,0}$. This is due to the fact that we need assumptions on $u_{h,0}$ which make it non-computable in practice. When we talk about “equivalence to time discretization of the semidiscrete problem” we always mean “up to changing the initial condition”, which incurs an additional (but easily treatable) error term.  *

Lemma 4.3**.**

Problem 4.1 is equivalent to solving the “standard” DG-formulation for the semidiscrete semigroup (3.5), i.e., if we define

[TABLE]

Then $u_{h,k}:=\operatorname{tr}(\mathscr{U}^{h}_{x,y,t})\in\SS^{\mathbf{r},0}(\mathcal{T}_{(0,T)})\otimes\mathbb{V}_{h}^{\mathcal{X}}$ solves

[TABLE]

On the other hand, we can recover the extended function by $\mathscr{U}^{h}_{x,y,t}:=\mathscr{D}_{h}u_{h,k}$.

Proof.

We first show that $\mathscr{D}_{h}u_{h,k}$ solves Problem 4.1.

Comparing the two formulations, the only interesting term is $\mathcal{A}(\mathscr{D}_{h}u_{h,k},\mathscr{V})$. We note that we can write:

[TABLE]

where we used that $\mathcal{A}(\mathscr{D}_{h}u_{h,k},\mathscr{W}_{h})=0$ vanishes for functions with $\operatorname{tr}{\mathscr{W}_{h}}=0$ by the definition of the lifting. Thus all the terms in the formulation directly correspond to each other.

We now show the other direction. Let $\mathscr{U}^{h}_{x,y,t}$ be a solution to Problem 4.1. We pick a function $q$, such that $q(t)=0$ outside of a single interval $(t_{j-1},t_{j})$ on which $q(t)$ is a polynomial. We then test (4.1) with functions of the form $\mathscr{V}_{h}(t):=q(t)\mathscr{V}_{0}$, where $\mathscr{V}_{0}\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ satisfies $\operatorname{tr}{\mathscr{V}_{0}}=0$. This means that $\mathscr{V}_{h}(t)\in\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y},\mathcal{T}}$ and we get, since all the terms involving $\operatorname{tr}{\mathscr{V}_{h}}$ vanish:

[TABLE]

Since $\mathscr{U}^{h}_{x,y,t}(t)$ is a polynomial of degree $r_{j}$ in $t$, $\mathcal{A}(\mathscr{U}^{h}_{x,y,t}(t),\mathscr{V}_{h})$ also is such a polynomial. Since the integral vanishes when tested with all similar polynomials, we get that $\mathcal{A}(\mathscr{U}(t),\mathscr{V}_{h})=0$ for all $t\in(t_{j-1},t_{j})$ and all admissible $\mathscr{V}_{0}$. This means we can write $\mathscr{U}^{h}_{x,y,t}=\mathscr{D}_{h}\operatorname{tr}{\mathscr{U}^{h}_{x,y,t}}$ and we can proceed as before to match all the terms in the formulation to their counterpart. ∎

Theorem 4.4** ($h$-version).**

Let $u_{h}$ denote the semidiscrete solution to (3.5). Suppose that Assumption 3.5 is fulfilled with $\beta>0$. Let $\mathbf{r}_{t}=r\equiv\text{const}$ be a fixed parameter. Choose $\mathcal{T}_{(0,T)}$ as a graded mesh with the grading function $h(t):=t^{\beta\left(2r+3\right)}$. Let $N:=\dim(\SS^{\mathbf{r}_{t},0}(\mathcal{T}_{(0,T)}))$.

Assume $u_{0}$ is analytic in $\overline{\Omega}$ and that the right-hand side $f$ satisfies

[TABLE]

with constants $C$ and $d$ independent of $\ell$ and $t$.

Then the following error estimate holds:

[TABLE]

The implied constant depends on $\Omega$, $u_{0}$, $f$, $r$, the terminal time $T$, and the constants from Assumption 3.5.

Proof.

We note that $u_{h,0}\in\mathcal{V}_{h,\beta}^{\mathcal{X}}$ by Assumption and also that the solution to DG-formulation depends continuously on the initial condition. This last statement can be easily seen from the coercivity of $\widetilde{B}$ as shown in [SS00, Lemma 2.7]. Thus, up to an additional error term $C(T)\left\|\Pi_{L^{2}}u_{0}-u_{h,0}\right\|_{L^{2}(\Omega)}^{2}$ we may use $u_{h,0}$ as our initial condition. (This error term is exponentially small by Assumption 3.5).

We want to apply the results from [SS00] and translate our setting into their requirements. They require separable Hilbert spaces $X\subseteq H$ with continuous, dense and compact embedding and a bilinear form $a(\cdot,\cdot):X\times X\to\mathbb{C}$, such that

[TABLE]

for all $u,v\in X$. We set $H:=\left(\mathbb{V}_{h}^{\mathcal{X}},\left\|\cdot\right\|_{L^{2}(\Omega)}\right)$, $X:=\left(\mathbb{V}_{h}^{\mathcal{X}},\left\|\cdot\right\|_{\mathbb{V}_{h}^{\mathcal{X}}}\right)$ and $a(u,v):=\left(\mathcal{L}^{s}_{h}u,v\right)_{L^{2}(\Omega)}$ (extending the real valued bilinear form to a complex one in the canonical way). By Lemma 3.3 this bilinear form satisfies the boundedness and ellipticity conditions. The symmetry follows from the definition and the symmetry of $\mathcal{A}(\cdot,\cdot)$.

The stated result then is a consequence of [SS00, Theorem 5.10]. The main ingredient is the fact that the initial condition is in the interpolation space $\mathcal{V}^{\mathcal{X}}_{h,\beta}$ by Assumption 3.5. Note that [SS00, Theorem 5.10] gives an estimate in the $\mathbb{V}_{h}^{\mathcal{X}}$-norm. In order to get to the more natural $\widetilde{H}^{s}(\Omega)$-norm, we use Lemma 3.3. ∎

Remark 4.5**.**

*For $r:=1$, the scheme in Theorem 4.4 is equivalent to the more common implicit Euler discretization, except that the right hand side is slightly modified. See [Tho06, Page 205] for details.  *

Theorem 4.6** ($hp$-version).**

Let $u_{h}$ denote the semidiscrete solution to (3.5). Consider $\mathcal{T}_{(0,T)}:=\mathcal{T}^{M}_{(\mathbf{0},t_{1})}\cup\mathcal{T}_{(t_{1},T)}$ to be a mesh on $(0,T)$ that is geometrically refined towards [math] and has constant size for larger times $(t_{1},T)$. We choose $\mathbf{r}_{t}$ such that it is linearly increasing on the geometrically refined part and constant afterwards. Let $N:=\dim(\SS^{\mathbf{r}_{t},0}(\mathcal{T}_{(0,T)}))$.

Assume that $u_{0}$ is analytic in $\overline{\Omega}$ and that the right-hand side $f$ satisfies

[TABLE]

with constants $C$ and $d$ independent of $\ell$ and $t$. Suppose that Assumption 3.5 is satisfied.

Then the following error estimate holds:

[TABLE]

The implied constant depends on $\Omega$, $u_{0}$, $f$, $\mu$, the mesh grading and the terminal time $T$ as well as the constants from Assumption 3.5.

Proof.

The proof is analogous to Theorem 4.4, except we now invoke [SS00, Section 5.1.2]. ∎

For the model problem of smooth geometries, we can give explicit bounds for the full discretization.

Corollary 4.7**.**

Assume that we are in the simplified setting of Section 3.3 and let the spaces for $\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ be designed as in Corollary 3.39. Denote the number of layers used in $\mathbb{V}_{h}^{\mathcal{Y}}$ as $M$. Assume that $u_{0}$ is analytic and $f$, $\dot{f}$ are uniformly analytic in a neighborhood $\widetilde{\Omega}\supset\overline{\Omega}$.

Let $\mathcal{T}_{(0,T)}:=\mathcal{T}^{M}_{(\mathbf{0},t_{1})}\cup\mathcal{T}_{(t_{1},T)}$ be a mesh on $(0,T)$ which is geometrically refined towards [math] with $M$ layers and has constant size for larger times $(t_{1},T)$. We chose $\mathbf{r}_{t}$ such that it is linearly increasing on the geometrically refined part and constant afterwards. We take $M\sim L$, where $L$ is the number of levels used for $\mathbb{V}_{h}^{\mathcal{Y}}$.

In addition, assume that the right-hand side $f$ satisfies

[TABLE]

with constants $C$ and $d$ independent of $\ell$ and $t$.

Then there exist constants $C>0$, $b>0$ such that the following error estimate holds:

[TABLE]

The implied constant depends on $u_{0}$, $f$, end time $T$, the domain $\Omega$, $\widetilde{\Omega}$, the mesh grading $\sigma$ as well as on $s$.

Proof.

Follows from Theorem 4.6, Theorem 3.38 and the fact that

[TABLE]

4.1 Practical aspects

In order to efficiently implement the scheme presented, we combine the Schur-form based approach described in [SS00] with the ideas of [BMN*+*18] for dealing with the extended variable.

For each time-inteval, the Schur decomposition in time leads to a sequence of problems of the form

[TABLE]

where $T\in\mathbb{C}^{r\times r}$ is an upper triangular matrix. These problems can be solved using a backward-substitution, where in each step an operator of the form $\frac{k}{\lambda_{j}}\mathcal{L}_{h}^{s}+\mathrm{I}$ has to be inverted. Structurally this is very similar to the operator $\mathcal{G}^{\lambda}$, except that the parameter $\lambda:=\lambda_{j}/k$ is complex valued. Proceeding like in [BMN*+*18] would require simultaneous diagonalization of the matrices

[TABLE]

Since the matrix $A$ is not hermitean if $\operatorname{Im}(\lambda_{j})\neq 0$, it is unclear whether this diagonalization can be done (in practice it appears to be the case). Instead we employ the generalized Schur-form (or QZ-decomposition; see [GVL96, Section 7.72]). It gives unitary matrices $Q$ and $Z$, such that $Q^{H}AZ=:T$ and $Q^{H}BZ=:S$ are both upper triangular. Inserting this decomposition into the definition of $\frac{k}{\lambda_{j}}\mathcal{L}^{s}_{h}+\mathrm{I}$ and using a backward-substitution leads to a sequence of problems of the form

[TABLE]

for $w\in H_{0}^{1}(\Omega)$ with $\kappa_{\ell}\in\mathbb{C}$.

Overall, Problem 4.1 can be solved by solving $\operatorname{dim}(\SS^{\mathbf{r}_{t}}(\mathcal{T}_{(0,T)}))\times\operatorname{dim}(\SS^{\mathfrak{r}}(\mathcal{T}^{M}_{(\mathbf{0},\mathcal{Y})}))$ scalar problems posed on $\Omega$. For the case of the geometric setting of Section 3.3 using the method described in Corollary 4.7, this means that $\mathcal{O}(M^{4})$ problems of size $\mathcal{O}(M^{2})$ need to be solved.

5 Numerical Results

In this section we test the theoretical findings of the previous sections by implementing them using the finite element package NGSolve [Sch14, Sch17] for the discretization in $\Omega$.

5.1 Smooth solution

In order to verify our implementation, we consider an example that has a known exact solution. We work with the simplified model problem of Section 3.3. Namely, working in 1D, we set $\Omega:=(0,1)$, $A=I$ and $c=1$. The initial condition is chosen as $u_{0}(x):=\sin(2\pi\,x)$. As an eigenfunction of the Dirichlet-Laplacian this leads to the exact solution $u(x,t):=e^{-t(2\pi)^{s}}\sin(2\pi\,x)$. We use $s=0.5$ and plot our findings, applying the $hp$-DG method. As seen in Figure 2(a), we get the predicted exponential convergence with respect to the number of refinement layers.

5.2 Singular solution

In order to verify that our method handles startup singularities robustly, we stay in the geometric setting of Section 5.1, but consider the initial condition $u_{0}\equiv 1$ and set $s:=0.75$. We use the trivial right-hand side $f\equiv 0$. Since the initial condition does not satisfy any compatibility condition, we expect startup singularities. As the exact solution is unknown, we precompute a numerical solution with high accuracy using the hp-DG method described in Corollary 4.7 with $M=13$ layers. We integrate up to the terminal time $T=1$. Due to the predicted exponential convergence, we expect a good match of the estimated error to the (unknown) true error.

We compare different time discretization schemes. For the implicit Euler based schemes we chose a fixed polynomial degree for discretizing $x$ and $y$ to be $p=8$. For the $hp-DG$ scheme we chose the same polynomial degree in each variable. As an indicator for comparing the numerical cost, we use the number of systems $N$ we need to solve involving the nonlocal operator $\mathcal{L}_{h}^{s}$. For the implicit Euler, this is proportional to the number of timesteps. For the $hp-DG$ approach it is proportional to the number of layers $M$ squared, i.e. $N\sim M^{2}$. In Figure 5.1 we compare the spacetime $L^{2}$-error to the number of such systems that need solving. We see that, as predicted, the implicit Euler method with a graded stepsize recovers the full convergence rate $\mathcal{O}(N^{-1})$ whereas a uniform approach only yields a reduced rate. It is important to point out that practical considerations may still favor using a uniform grid, as in this case the corresponding matrices can be factorized only once. This yields much faster solution times in each step. Since the reduction of order is small, the uniform approach often outperforms the graded mesh in our experience.

The best performance, as expected, is observed by the $hp-DG$ based method. It provides rapid exponential convergence of order $\mathcal{O}(e^{-b\sqrt{N}})$, confirming Theorem 4.6 and Corollary 4.7.

5.3 A 2d example

Although our theory in the 2D case is restricted to domains with an analytic boundary, we show numerically that the case of polygons can be successfully treated as well. We chose $\Omega:=(0,1)^{2}$, $u_{0}\equiv 1$, $f\equiv 0$, $A:=\operatorname{I}$, $c=0$ and $s:=1/4$. Since no known analytic solution is avaliable, we computed the approximation using $M=10$ levels of refinement in time and used it as our reference solution. All computations were done up to the terminal time $T=1$ and using the $hp$-DG method. For the time discretization and discretization in $y$, we used a geometric grid with $M$ layers. In $\Omega$ we used a geometrically refined grid of $3M/2$ layers in accordance to Corollary 3.26.

In Figure 2(b), we see that also in this case we get the exponential convergence with respect to the number of layers in the $hp$-refinement. This suggests that our methods could also be extended to cover this case.

Appendix A Exponential convergence of $hp$-FEM for singularly perturbed problems with complex coefficients

In this appendix we provide the details for the regularity and approximability by high order FEM on suitably designed meshes for singularly perturbed problems with a complex perturbation parameter. We consider on a domain $\Omega\subset\mathbb{R}^{d}$, $d\in\{1,2\}$ the problem of finding $u_{\varepsilon}\in H_{0}^{1}(\Omega)$ such that

[TABLE]

Concerning the data of this problem, we make the following assumption:

Assumption A.1**.**

The domain $\Omega\subset\mathbb{R}^{d}$, $d\in\{1,2\}$ has an analytic boundary and $\widetilde{\Omega}$ is a domain with $\overline{\Omega}\subset\widetilde{\Omega}$.

The parameter $\varepsilon$ satisfies $\varepsilon>0$, and the matrix valued function $A\in L^{\infty}(\Omega,\mathbb{R}^{d\times d})$ is analytic on $\widetilde{\Omega}$, pointwise and uniformly SPD. The function $f\in L^{2}(\Omega)$ is analytic on $\widetilde{\Omega}$. The function $\zeta\in L^{\infty}(\Omega)$ is such that:

1. (i)

\zeta(x,y)\in\widetilde{\mathcal{S}}:=\big{\{}z\in\mathbb{C}:\left|\pi-\operatorname{Arg}(z)\right|\geq\delta>0\big{\}}, 2. (ii)

$\operatorname{Im}(\zeta)\in\mathbb{R}$* is constant,* 3. (iii)

$\left|\zeta(x,y)\right|\geq\zeta_{0}>0$* in $\Omega$,* 4. (iv)

$\zeta$* is analytic in $\widetilde{\Omega}$.*

Associated with the operator $L_{\varepsilon}$ is the sesquilinear form

[TABLE]

and the energy norm $\left\|\cdot\right\|^{2}_{\varepsilon}:=\varepsilon^{2}\left\|\nabla\cdot\right\|_{L^{2}(\Omega)}^{2}+\left\|\cdot\right\|_{L^{2}(\Omega)}^{2}$.

Lemma A.2**.**

Let Assumption A.1 be valid. Then the bilinear form $a_{\varepsilon}(\cdot,\cdot)$ is bounded and elliptic in the energy norm, i.e., there exists $\theta(\zeta)\in(-\pi,\pi)$ such that

[TABLE]

The implied constants do not depend on $\varepsilon$ or $u$. This implies for the solution $u_{\varepsilon}$ to (A.1):

[TABLE]

Proof.

For $\alpha\in\mathbb{C}$ with $|\alpha|=1$, we compute:

[TABLE]

Thus it remains to show that we can choose $\alpha$ such that $\operatorname{Re}(\alpha)>0$ and $\operatorname{Re}(\alpha\zeta)>0$ uniformly in $\Omega$. If $\operatorname{Im}(\zeta)\geq 0$, we can pick $\alpha:=e^{-i\frac{\pi-\delta}{2}}$, otherwise $\alpha:=e^{i\frac{\pi-\delta}{2}}$ does the trick. The estimate (A.2) follows from the Lax-Milgram lemma. ∎

The previous lemma ensures existence and uniqueness of solutions $u_{\varepsilon}$. In the next one we further prove that $u_{\varepsilon}$ is analytic with explicit bounds on the derivative with respect to the parameter $\varepsilon$.

Lemma A.3**.**

Let Assumption A.1 be valid. Let $u_{\varepsilon}\in H^{1}_{0}(\Omega)$ solve (A.1) Then $u_{\varepsilon}$ is analytic on $\overline{\Omega}$ and satisfies:

[TABLE]

Proof.

The statement is the restriction of [Mel02, Theorem 2.3.1] (for 2D, in 1D the relevant result is Proposition 2.2.1) to the case of smooth domains. While this reference only considers $\zeta>0$, the proof carries over almost verbatim. The only modification needed is the coercivity estimate from Lemma A.2. ∎

While Lemma A.3 will provide exponential convergence in the asymptotic case of sufficiently large polynomial degree, the more practically relevant regime is treated using the following lemma:

Lemma A.4**.**

Let Assumption A.1 be valid. Let $u_{\varepsilon}\in H^{1}_{0}(\Omega)$ solve (A.1) for $\varepsilon\in(0,1]$. Then there exists a smooth cut-off function $\chi$ supported by a tubular neighborhood of $\partial\Omega$ with $\chi\equiv 1$ in a neighborhood of $\partial\Omega$ and constants $C$, $\gamma$, $b>0$ independent of $\varepsilon\in(0,1]$ such that $u_{\varepsilon}$ can be decomposed as

[TABLE]

with the following properties:

1. (i)

The smooth part $w_{\varepsilon}$ is analytic in $\Omega$ and satisfies $\left\|\nabla^{p}w_{\varepsilon}\right\|_{L^{\infty}(\Omega)}\leq C\gamma^{p}p!$ for all $p\in\mathbb{N}_{0}$. 2. (ii)

The remainder $r_{\varepsilon}\in H_{0}^{1}(\Omega)$ satisfies $\left\|r_{\varepsilon}\right\|_{\varepsilon}+\left\|r_{\varepsilon}\right\|_{H^{1}(\Omega)}+\left\|\nabla^{2}r_{\varepsilon}\right\|_{L^{2}(\Omega)}\lesssim Ce^{-b/\varepsilon}$. 3. (iii)

Using boundary fitted coordinates $(\rho,\theta)$, where $\rho=\operatorname{dist}(\cdot,\partial\Omega)$ and $\theta$ is a parametrization of $\partial\Omega$, the boundary layer $u_{\varepsilon}^{\textrm{BL}}$ can be estimated

[TABLE]

Proof.

We focus on the 2D case by adapting [Mel02, Theorem 2.3.4] to the case of smooth geometries and complex data. The 1D case follows from adapting [Mel02, Lemma 7.1.1] instead.

While the somewhat technical proof from [Mel02] only considers real data $\zeta>0$, it can be adapted to our setting in a mostly straight forward way. We make some comments on how to read the proof and how to make the required modifications.

The construction is laid out in [Mel02, Section 7]. The smooth part is constructed inductively:

[TABLE]

The estimate (i) then follows as in [Mel02, Lemma 7.2.1] from Cauchy’s integral theorem. The fact that $\zeta$ is complex does not require modifications, we only need that the function $\frac{1}{\zeta}$ has an analytic extension to a neighborhood of $\Omega$. This is guaranteed by the assumption $\left|\zeta\right|\geq\zeta_{0}>0$.

In order to construct the boundary layer function and prove (iii), one works in boundary adapted coordinates. Proceeding as in [Mel02, Section 7.3.1], $u_{\varepsilon}^{\textrm{BL}}$ is defined via

[TABLE]

where $\widehat{U}_{i}$, using $\lambda:=\sqrt{\frac{\zeta(\theta,0)}{A_{22}(\theta,0)}}$, solves an ODE of the form

[TABLE]

The necessary estimates of [Mel02, Section 7] to conclude (iii) all rely on [Mel02, Lemma 7.3.6] which gives exponential decay for problems of the form (A.4). It is already formulated for complex parameters $\lambda$, we only point out that due to the assumption that $\zeta\in\widetilde{\mathcal{S}}$, we get $\operatorname{Re}(\lambda)>0$ (using the principal branch of the complex square root, satisfying $\operatorname{Re}(z)\geq 0\;\;\forall z\in\mathbb{C}$). The requirement $\operatorname{Re}(\lambda^{2})>0$ made in [Mel02, Lemma 7.3.6] is not satisfied, but inspection of the proof reveals that it is only needed to get unique solvability of (A.4). As seen in Lemma A.2, this is also guaranteed in the current setting.

Finally, (ii) follows form the fact that $r_{\varepsilon}:=u_{\varepsilon}-w_{\varepsilon}-u_{\varepsilon}^{\textrm{BL}}$, solves $L_{\varepsilon}r_{\varepsilon}=f_{\varepsilon}$ where $\left\|f_{\varepsilon}\right\|_{L^{\infty}(\Omega)}$ is exponentially small (see [Mel02, Eqn (7.4.37)] and $r_{\varepsilon}|_{\Gamma}=0$. The stated estimate then again follows from standard a priori estimates, most notably Lemma A.2. ∎

Lemma A.5**.**

Let Assumption A.1 be valid. Let $u_{\varepsilon}$ solve (A.1), let $\mathcal{T}^{L}_{\Omega}$ be an anisotropic geometric mesh refined towards $\partial\Omega$ as in Definition 3.30 and $\mathbb{V}_{h}^{\mathcal{X}}$ be the space of continuous piecewise polynomials of degree $p$ (see (3.27)). Assume that $\sigma^{L}\leq\varepsilon\leq 1$.

Then there exist constants $C$, $b>0$ such that for all $p\in\mathbb{N}$

[TABLE]

Proof.

Analogously to [BMN*+*18, Theorem 7.7], we note that the mesh $\mathcal{T}^{L}_{\Omega}$ contains a so-called admissible boundary layer mesh, i.e. a mesh containing one layer of “needle elements” of size $\lambda p\varepsilon$ (see[MS98] or[Mel02] for the precise definition).

The proof then follows completely analogously to [MS98, Theorem 3.14] (or also [Mel02, Theorem 3.4.8]). The necessary ingredients to generalize to complex parameters $\zeta$ (as described in Assumption A.1) are given by Lemmas A.3 and A.4. ∎

Appendix B Polynomial liftings and interpolation spaces

In this section, we investigate under which conditions we can lift discrete functions from $\mathbb{V}_{h}^{\mathcal{X}}$ to functions in $\mathbb{V}_{h}^{\mathcal{X},\mathcal{Y}}$ in a stable way. This question is deeply related to the theory of interpolation of discrete polynomial spaces. This can be seen in the following proposition:

Proposition B.1** ([Tar07, Lemma 40.1]).**

Let $X_{1}\subseteq X_{1}$ be Banach spaces with continuous embedding. For $\theta\in(0,1)$, denote the interpolation space by $X_{\theta}:=[X_{0},X_{1}]_{\theta,2}$. Then the following statements hold:

1. (i)

If $v$ is a $X_{0}$-valued function such that $v(t)\in X_{1}$ and $\dot{v}(t)\in X_{0}$ for all $t>0$ and $t^{1-\theta}\left\|\dot{v}(t)\right\|_{X_{0}}\in L^{2}(\mathbb{R}_{+},\frac{dt}{t})$, $t^{1-\theta}\left\|v(t)\right\|_{X_{1}}\in L^{2}(\mathbb{R}_{+},\frac{dt}{t})$ then $v(0)\in X_{\theta}$ with

[TABLE] 2. (ii)

If $v_{0}\in X_{\theta}$, there exists a function $v:\mathbb{R}_{+}\to X_{1}$ such that $v(0)=v_{0}$ and

[TABLE]

Proof.

This is just a special case of [Tar07, Lemma 40.1]. We note that in comparison to the statement in the book we changed the roles of $X_{0}$ and $X_{1}$. But since $\displaystyle[X_{1},X_{0}]_{\theta,2}=[X_{0},X_{1}]_{1-\theta,2}$ by [Tar07, Lemma 25.4], the theorem holds in the stated form. ∎

The case of lifting a polynomial on a single element $[0,1]$ to the unit square was addressed in [BDM07]. Namely, the following holds:

Proposition B.2** ([BDM07]).**

For $d=1$, let $\mathcal{Q}^{p}:=\operatorname{span}\{x_{i},\,0\leq i\leq p\}$ denote the space of polynomials. Then the following statements hold:

1. (i)

The interpolation norm coincides with the Sobolev norm, i.e., for all $\theta\in(0,1)$

[TABLE]

with equivalent norms. The implied constant depends only on $\theta$. 2. (ii)

For all $u\in\mathcal{Q}^{p}$, there exists a polynomial \mathscr{U}\in\mathcal{Q}^{p}([0,1]^{2}):=\operatorname{span}_{0\leq i,j\leq p}\big{\{}x^{i}\,y^{j}\big{\}} such that $\operatorname{tr}{\mathscr{U}}=u$, $\mathscr{U}(\cdot,y)\in\mathcal{Q}^{p}$ for all $y\in[0,1]$. For $y>1$, $\mathscr{U}$ can be extended by [math] to $\mathbb{R}_{+}$ such that $\left\|\mathscr{U}\right\|_{H^{1}(y^{\alpha},D)}\lesssim\left\|u\right\|_{{H}^{s}([0,1])}$.

Proof.

See Corollary 4.4 and Corollary 3.3 in Chapter II of [BDM07]. ∎

When working in 2d, it is not sufficient to work only with global polyomials. It is possible to generalize Proposition B.2 to the case of piecewise polynomials on a shape regular mesh. In 2d, this is worked out in [MKR19]:

Proposition B.3** ([MKR19]).**

Let $\Omega\subseteq\mathbb{R}^{2}$ and $\SS^{p,1}(\mathcal{T}_{\Omega})$ denote the space of piecewise polynomials on a shape regular grid of quadrilaterals (see Definition 3.29).

Then the following statements hold:

1. (i)

The interpolation norm coincides with the Sobolev norm, i.e., for all $\theta\in(0,1)$

[TABLE]

with equivalent norms. The implied constant depends only on $\theta$. 2. (ii)

For all $u\in\SS^{p,1}(\mathcal{T}_{\Omega})$, there exists a function \mathscr{U}\in C\big{(}\Omega\times[0,1]\big{)} such that $\operatorname{tr}{\mathscr{U}}=u$, $\mathscr{U}(\cdot,y)\in\SS^{p,1}(\mathcal{T}_{\Omega})$ for all $y\in[0,1]$. For $y>1$, $\mathscr{U}$ can be extended by [math] to $\mathbb{R}_{+}$ such that $\left\|\mathscr{U}\right\|_{H^{1}(y^{\alpha},D)}\lesssim\left\|u\right\|_{{H}^{s}(\Omega)}$.

The previous propositions give a lifting to the space of either polynomials or continuous functions in the extended variable $y$. Since we will be working with piecewise polynomials with a linear degree vector neither is sufficient for our needs. We need the following variation of the previous result:

Lemma B.4**.**

Let $u\in\mathcal{Q}^{p}$ in $1d$ or $u\in\SS^{p,1}(\mathcal{T}_{\Omega})$, where $\mathcal{T}_{\Omega}$ is a shape regular mesh of quadrilaterals as defined in Definition 3.29. Assume that the triangulation $\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})}$ satisfies $\operatorname{diam}(K_{0})\leq h_{x}p^{-2}$, where $K_{0}$ is the element at [math] and $h_{x}:=\min_{K\in\mathcal{T}_{\Omega}}{\operatorname{diam}(K)}$. Then there exists a lifting $\mathscr{U}_{h}\in\SS^{p,1}(\mathcal{T}_{\Omega})\otimes\mathbb{V}_{h}^{\mathcal{Y}}$ such that

[TABLE]

The constant $C$ depends only on $s$ and the mesh grading parameter $\sigma$. The lifting can be chosen to be piecewise linear with respect to $y$.

Proof.

We focus on the 2d case. By Propositions B.1 and B.3 , there exists a lifting $\mathscr{U}\in C(\mathbb{R}_{+},\SS^{p,1}(\mathcal{T}_{\Omega}))$ such that

[TABLE]

Inspecting the proof of Proposition B.1, as given in [Tar07], one can see that the lifting $\mathscr{U}$ is piecewise linear on the grid \big{(}e^{n}\big{)}_{n\in\mathbb{Z}}. By a simple rescaling, we may choose $\mathscr{U}$ as piecewise linear in $y$ on the geometric mesh $\sigma^{n}$ for $n\in\mathbb{Z}$. To get a function which is in the space $\SS^{1,1}(\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})})$, we need to make two modifications: modify $\mathscr{U}$ on the element $K_{0}:=(0,\sigma^{L})$ to also be linear and cut the function off at $\mathcal{Y}$. We define $h_{0}:=\operatorname{diam}(K_{0})=\sigma^{L}$.

We define $\mathscr{U}_{h}(\cdot,t)$ as the linear interpolation between $u=\mathscr{U}(0)$ and $\mathscr{U}(\sigma^{L})$ on $K_{0}$ and $\mathscr{U}_{h}=\mathscr{U}$ otherwise. We need to show:

[TABLE]

We start with the first inequality. Since $\mathscr{U}_{h}$ is the linear interpolant of $\mathscr{U}$, we can write $\partial_{y}\mathscr{U}_{h}=h_{0}^{-1}\int_{0}^{h_{0}}{\partial_{y}\mathscr{U}(\tau)\,d\tau}$. This gives:

[TABLE]

We first investigate the term $I_{1}$. Using the fact that $y\leq h_{0}$ and therefore $h_{0}^{-2}\leq y^{-2}$, we get

[TABLE]

by Hardy’s inequality(see [Gri85, page 28]).

When investigating $I_{2}$, we distinguish $\alpha\geq 0$ and $\alpha\leq 0$. For $\alpha\geq 0$ we have $y^{\alpha}\leq\tau^{\alpha}$ for $y\leq\tau$ and thus after applying Cauchy Schwarz to get the square into the integral:

[TABLE]

For $\alpha\leq 0$, we have $h_{0}^{\alpha}\leq\tau^{\alpha}$ and get:

[TABLE]

which proves (B.1).

We now show (B.2). The proof relies on an inverse estimate and the fact that $\mathscr{U}_{h}$ approximates $\mathscr{U}$. We estimate:

[TABLE]

Since $\mathscr{U}_{h}(\cdot,y)$ and $\mathscr{U}(\cdot,y)$ are (piecewise) polynomials for all fixed $y$, we can use an inverse estimate [Sch98, Theorem 3.91] to get:

[TABLE]

Since $\mathscr{U}_{h}-\mathscr{U}$ vanishes at $y=0$, we can write it as $\mathscr{U}_{h}(y)-\mathscr{U}(y)=\int_{0}^{y}{\partial_{y}{\mathscr{U}}_{h}(\tau)-\partial_{y}{\mathscr{U}}(\tau)\,d\tau}$ and further estimate:

[TABLE]

The term $I_{3}$ is structurally analogous to the term (B.2) and can be estimated using the same techniques. The extra integration in $\tau$ gives an additional power of $h_{0}^{2}$, and we get

[TABLE]

For the term $I_{4}$, we apply Hardy’s inequality and the estimate $h_{0}^{-2}\leq y^{-2}$ to get:

[TABLE]

Overall, since we assumed $h_{0}\leq h_{x}\,p^{-2}$, we get the stability of the modified lifting.

In order to get $\operatorname{supp}{\mathscr{U}_{h}}\subset[0,\mathcal{Y}]$, we pick the cutoff function $\varphi\in\SS^{1,1}(\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})})$ such that $\varphi|_{K_{i}}=1$ on $K_{i}$ for $i=0,\dots,|\mathcal{T}^{L}_{(\mathbf{0},\mathcal{Y})}|-1$ and $\varphi(\mathcal{Y})=0$. We note that the element where $\varphi$ is non-constant has size $\mathcal{O}(1)$, and it can be easily checked that $\varphi\cdot\mathscr{U}_{h}$ is also a stable lifting of $u$. In order to get a function in \SS^{1,1}\big{(}\mathcal{T}^{L}_{\mathcal{Y}},\SS^{p,1}(\mathcal{T}_{\Omega})\big{)} we interpolate the function in the grid points. Since $\mathscr{U}_{h}\cdot\varphi$ is a polynomial of degree at most $2$, interpolating it down to degree 1 is stable in the $L^{2}$ and $H^{1}$ norm (see [BM97, Rem. 13.5 and (13.27)]). Away from [math], the weighted norms are equivalent to the standard ones. This shows that the “cutoff and interpolation”-procedure is stable in $\mathring{H}^{1}(y^{\alpha},\mathcal{C})$. ∎

Appendix C Proof of Proposition 3.12

The following proof consists of condensed and restated results from [Tho06, Chapter 3]. We fix $t_{0}>0$ and consider the discrete backward problem

[TABLE]

For $\tau\in(0,t_{0})$, we get by testing (C.1) with $\theta$ in the $L^{2}$-inner product and using (3.7) and (3.6):

[TABLE]

For $0<\varepsilon<t_{0}$ we get by integrating, since $z_{h}(t_{0})=0$:

[TABLE]

In the limit $\varepsilon\to 0$, the first term converges due to Lemma 3.10 (i) to

[TABLE]

The following stability estimate holds for $z_{h}$ by Lemma 3.10 (iii):

[TABLE]

Combining this estimate with (C.2) completes the proof of (3.8).

Proof of (3.11): For fixed $t>0$, testing the equation (3.7) with $v:=t\,\theta(t)$ and integrating over $\Omega$ gives:

[TABLE]

We integrate in $t$ from $\varepsilon>0$ to $t$ and get:

[TABLE]

We need to bound $\lim_{\varepsilon\to 0}{\varepsilon\left\|\theta(\varepsilon)\right\|^{2}_{L^{2}(\Omega)}}$. Writing $\theta=\Pi_{h}u-u_{h}=\rho+u-u_{h}$, we use the fact that $u_{h}$ and $u$ are bounded by Lemma 3.10 (i). This gives:

[TABLE]

By using Young’s inequality and (3.8), we easily obtain (3.11) from the fact that $\left\|\theta(t)\right\|^{2}_{\widetilde{H}^{s}(\Omega)}\lesssim\left(\mathcal{L}_{h}^{s}\theta(t),\theta(t)\right)_{L^{2}(\Omega)}$ by Lemma 3.3.

Acknowledgments: Financial support by the Austrian Science Fund (FWF) through the research program “Taming complexity in partial differential systems” (grant SFB F65, A.R.).