Superconvergent interpolatory HDG methods for reaction diffusion equations I: An HDG$_{k}$ method

TL;DR

This paper enhances an interpolatory HDG method for reaction diffusion equations by enabling superconvergence through a modified postprocessing step, maintaining computational efficiency and improving accuracy.

Contribution

It introduces a simple modification to the interpolatory HDG method that restores superconvergence for reaction diffusion problems, combining efficiency with higher accuracy.

Findings

Restores superconvergence with a modified postprocessing step.

Maintains computational efficiency of the original interpolatory HDG method.

Numerical results confirm theoretical convergence and performance improvements.

Abstract

In our earlier work [8], we approximated solutions of a general class of scalar parabolic semilinear PDEs by an interpolatory hybridizable discontinuous Galerkin (Interpolatory HDG) method. This method reduces the computational cost compared to standard HDG since the HDG matrices are assembled once before the time integration. Interpolatory HDG also achieves optimal convergence rates; however, we did not observe superconvergence after an element-by-element postprocessing. In this work, we revisit the Interpolatory HDG method for reaction diffusion problems, and use the postprocessed approximate solution to evaluate the nonlinear term. We prove this simple change restores the superconvergence and keeps the computational advantages of the Interpolatory HDG method. We present numerical results to illustrate the convergence theory and the performance of the method.

Equations

$$\begin{array}{rl}{\partial }_{t}u-\Delta u+F(u)& =f\text{mbox}in\Omega \times (0,T],\\ u& =0\text{mbox}on\partial \Omega \times (0,T],\\ u(\cdot ,0)& =u_0\text{mbox}in\Omega .\end{array}$$

$$(w,v{)}_{{\mathcal{T}}_h}:=\sum _{K\in {\mathcal{T}}_h}(w,v{)}_K,{⟨\zeta ,\rho ⟩}_{\partial {\mathcal{T}}_h}:=\sum _{K\in {\mathcal{T}}_h}{⟨\zeta ,\rho ⟩}_{\partial K},$$

$${\mathbf{V}}_h$$

$$W_h$$

$$Z_h$$

$$M_h$$

$$({\mathbf{q}}_h,{\mathbf{r}}_h{)}_{{\mathcal{T}}_h}-(u_h,\nabla \cdot {\mathbf{r}}_h{)}_{{\mathcal{T}}_h}+{⟨{\hat{u}}_h,{\mathbf{r}}_h\cdot \mathbf{n}⟩}_{\partial {\mathcal{T}}_h}$$

$$({\partial }_t{u}_h,v_h{)}_{{\mathcal{T}}_h}-({\mathbf{q}}_h,\nabla v_h{)}_{{\mathcal{T}}_h}+{⟨{\hat{\mathbf{q}}}_h\cdot \mathbf{n},v_h⟩}_{\partial {\mathcal{T}}_h}+({\mathcal{I}}_{h}F(u_h^{\star }),v_h{)}_{{\mathcal{T}}_h}$$

$${⟨{\hat{\mathbf{q}}}_h\cdot \mathbf{n},{\hat{v}}_h⟩}_{\partial {\mathcal{T}}_h\backslash {\epsilon }_h^{\partial }}$$

$$u_h(0)$$

$${\hat{\mathbf{q}}}_h\cdot \mathbf{n}$$

$$(\nabla {\mathfrak{q}}_h^{k+1}({\mathbf{q}}_h,u_h),\nabla z_h{)}_K$$

$$({\mathfrak{q}}_h^{k+1}({\mathbf{q}}_h,u_h),w_h{)}_K$$

$$[{\mathcal{P}}^{k+1}(K){]}^{\perp }=\{z\in {\mathcal{P}}^{k+1}(K):(z,w_0{)}_K=0\text{for all }{w}_0\in {\mathcal{P}}^0(K)\}.$$

$${\partial }_t^{+}{u}_h^n=\frac{u_h^n-u_h^{n-1}}{\Delta t}.$$

$$({\mathbf{q}}_h^n,\mathbf{r}{)}_{{\mathcal{T}}_h}-(u_h^n,\nabla \cdot \mathbf{r}{)}_{{\mathcal{T}}_h}+{⟨{\hat{u}}_h^n,\mathbf{r}\cdot \mathbf{n}⟩}_{\partial {\mathcal{T}}_h}$$

$$({\partial }_t^{+}{u}_h^n,w{)}_{{\mathcal{T}}_h}-({\mathbf{q}}_h^n,\nabla w{)}_{{\mathcal{T}}_h}+{⟨{\hat{\mathbf{q}}}_h^n\cdot \mathbf{n},w⟩}_{\partial {\mathcal{T}}_h}+({\mathcal{I}}_{h}F(u_h^{n\star }),w{)}_{{\mathcal{T}}_h}$$

$${⟨{\hat{\mathbf{q}}}_h^n\cdot \mathbf{n},\mu ⟩}_{\partial {\mathcal{T}}_h\backslash {\epsilon }_h^{\partial }}$$

$$u_h^0$$

$${\hat{\mathbf{q}}}_h^n\cdot \mathbf{n}={\mathbf{q}}_h^n\cdot \mathbf{n}+\tau (u_h^n-{\hat{u}}_h^n),$$

$$(\nabla u_h^{n\star },\nabla z_h{)}_K$$

$$(u_h^{n\star },w_h{)}_K$$

$$(\nabla u_h^{n\star },\nabla z_h{)}_K+({\eta }_h^n,z_h{)}_K$$

$$(u_h^{n\star },w_h{)}_K$$

$$\begin{array}{r}{\mathbf{q}}_h^n=\sum _{j=1}^{N_1}{\alpha }_j^n{\mathbf{\phi }}_j,u_h^n=\sum _{j=1}^{N_2}{\beta }_j^n{\varphi }_j,\\ u_h^{n\star }=\sum _{j=1}^{N_3}{\gamma }_j^n{\chi }_j,{\hat{u}}_h^n=\sum _{j=1}^{N_4}{\zeta }_j^n{\psi }_j.\end{array}$$

$$A_1$$

$$A_4$$

$$A_7$$

$$M$$

$$b_1$$

$$\displaystyle{\begin{bmatrix}A_{1}^{k}&(b_{1}^{k})^{T}\\ b_{1}^{k}&0\end{bmatrix}}{\left[{\begin{array}[]{*{20}{c}}\bm{\gamma}^{n}_{k}\\ \bm{\eta}^{n}_{k}\\ \end{array}}\right]}={\begin{bmatrix}-A_{2}^{k}&0\\ 0&b_{2}^{k}\end{bmatrix}}{\left[{\begin{array}[]{*{20}{c}}\bm{\alpha}^{n}_{k}\\ \bm{\beta}^{n}_{k}\end{array}}\right]},$$

$$\displaystyle{\left[{\begin{array}[]{*{20}{c}}\bm{\gamma}^{n}_{k}\\ \bm{\eta}^{n}_{k}\\ \end{array}}\right]}$$

$$\displaystyle={\begin{bmatrix}B_{11}^{k}&B_{12}^{k}\\ B_{21}^{k}&B_{22}^{k}\end{bmatrix}}{\left[{\begin{array}[]{*{20}{c}}\bm{\alpha}^{n}_{k}\\ \bm{\beta}^{n}_{k}\end{array}}\right]},$$

$${\mathbf{\gamma }}_k^n=B_{11}^k{\mathbf{\alpha }}_k^n+B_{12}^k{\mathbf{\beta }}_k^n.$$

$$[({\mathcal{I}}_{h}F(u_h^{n\star }),{\varphi }_i{)}_{{\mathcal{T}}_h}]$$

$$=A_9\mathcal{F}(B_{11}{\mathbf{\alpha }}^n+B_{12}{\mathbf{\beta }}^n),$$

$$\mathcal{F}({\mathbf{\gamma }}^n)$$

$$\displaystyle\underbrace{\begin{bmatrix}A_{3}&-A_{4}&A_{5}\\ A_{4}^{T}&\Delta t^{-1}M+A_{6}&-A_{7}\\ A_{5}^{T}&A_{7}^{T}&-A_{8}\end{bmatrix}}_{K}\underbrace{\left[{\begin{array}[]{*{20}{c}}\bm{\alpha}^{n}\\ \bm{\beta}^{n}\\ \bm{\zeta}^{n}\\ \end{array}}\right]}_{\bm{x}_{n}}+\underbrace{\left[{\begin{array}[]{*{20}{c}}0\\ A_{9}\mathcal{F}(B_{11}\bm{\alpha}^{n}+B_{12}\bm{\beta}^{n})\\ 0\end{array}}\right]}_{\mathscr{F}(\bm{x}_{n})}=\underbrace{\left[{\begin{array}[]{*{20}{c}}0\\ b_{3}^{n}+{\Delta t}^{-1}M\bm{\beta}^{n-1}\\ 0\end{array}}\right]}_{\bm{b}_{n}},$$

$$K{\mathbf{x}}_n+\mathcal{F}({\mathbf{x}}_n)={\mathbf{b}}_n.$$

Full text

Keywords Interpolatory hybridizable discontinuous Galerkin method, superconvergence

1 Introduction

In our earlier work [8], we introduced an interpolatory hybridizable discontinuous Galerkin (Interpolatory HDG) method to approximate the solution of semilinear parabolic PDEs. In contrast to standard HDG, the Interpolatory HDG method uses an elementwise interpolation procedure to approximate the nonlinear term; therefore, all quadrature for the nonlinear term can be performed once before the time integration, which yields a significant computational cost reduction. The Interpolatory HDG method still converged at optimal rates, but superconvergence using element-by-element postprocessing was lost.

The superconvergence is an excellent feature of HDG methods, and therefore in this work we modify the Interpolatory HDG method from [8] and restore the superconvergence for reaction diffusion PDEs.

Specifically, we consider the following class of scalar reaction diffusion PDEs on a Lipschitz polyhedral domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, with boundary $\partial\Omega$:

[TABLE]

In Section 2, we provide background on HDG methods and describe the new Interpolatory HDG approach in detail. We use the HDGk method to approximate the linear terms in the equation; i.e., $k$th order discontinuous polynomials are used to approximate the flux $\bm{q}=-\nabla u$, the scalar variable $u$, and its trace, and the stabilization function is chosen as $O(1)$ piecewise constant. For the nonlinear term, we again use an elementwise Lagrange interpolation operator, as in [8], but now we also approximate $u$ using a postprocessing approach. This modified approximate nonlinearity restores the superconvergence and, as in [8], we have a simple explicit expressions for the nonlinear term and Jacobian matrix, which leads to an efficient and unified implementation.

We analyze the semidiscrete Interpolatory HDGk method in Section 3. We first assume the nonlinearity satisfies a global Lipschitz condition and prove the superconvergence. Next, we establish the superconvergence under a local Lipschitz condition, assuming the mesh is quasi-uniform.

In Section 4, we illustrate the convergence theory with numerical experiments and also demonstrate the performance of the Interpolatory HDGk method on a reaction diffusion PDE system.

We note that interpolatory finite element methods for nonlinear PDEs are well-known to have computational advantages and have a long history. The approach has been given many different names, including finite element methods with interpolated coefficients, product approximation, and the group finite element method. For more information, see [18, 7, 20, 19, 24, 32, 35, 39, 38, 4, 6, 23, 22, 40, 34, 41, 36, 17] and the references therein.

2 Interpolatory HDGk formulation and implementation

Hybridizable discontinuous Galerkin (HDG) methods were proposed by Cockburn et al. in [13]. HDG methods work with the mixed formulation of the PDE, and on each element the approximate solution and flux are expressed in terms of the approximate solution trace on the element boundary. The approximate trace is uniquely determined by requiring the normal component of the numerical trace of the flux to be continuous across element boundaries. This allows the approximate solution and approximate flux variables to be eliminated locally on each element; the result is a global system of equations for the approximate solution trace only. Therefore, the number of globally coupled degrees of freedom for HDG methods is significantly lower than for standard DG methods. HDG methods have been successfully applied to linear PDEs [13, 14, 15, 27] and nonlinear PDEs [29, 31, 26, 25, 2, 30, 21, 28, 16].

To describe the Interpolatory HDGk method, we introduce notation below. We mostly follow the notation used in [13], where HDG methods were considered for linear, steady-state diffusion.

Let $\mathcal{T}_{h}$ be a collection of disjoint simplexes $K$ that partition $\Omega$. Let $\partial\mathcal{T}_{h}$ denote the set $\{\partial K:K\in\mathcal{T}_{h}\}$. For an element $K$ in the collection $\mathcal{T}_{h}$, let $e=\partial K\cap\Gamma$ denote the boundary face of $K$ if the $d-1$ Lebesgue measure of $e$ is nonzero. For two elements $K^{+}$ and $K^{-}$ of the collection $\mathcal{T}_{h}$, let $e=\partial K^{+}\cap\partial K^{-}$ denote the interior face between $K^{+}$ and $K^{-}$ if the $d-1$ Lebesgue measure of $e$ is nonzero. Let $\varepsilon_{h}^{o}$ and $\varepsilon_{h}^{\partial}$ denote the sets of interior and boundary faces, respectively, and let $\varepsilon_{h}$ denote the union of $\varepsilon_{h}^{o}$ and $\varepsilon_{h}^{\partial}$. We use the mesh-dependent inner products

[TABLE]

where $(\cdot,\cdot)_{D}$ denotes the $L^{2}(D)$ inner product for a set $D\subset\mathbb{R}^{d}$ and $\langle\cdot,\cdot\rangle_{\Gamma}$ denotes the $L^{2}(\Gamma)$ inner product for a set $\Gamma\subset\mathbb{R}^{d-1}$.

Let $\mathcal{P}^{k}(D)$ denote the set of polynomials of degree at most $k$ on a domain $D$. We consider the discontinuous finite element spaces

[TABLE]

All spatial derivatives of functions in these spaces should be understood piecewise on each element $K\in\mathcal{T}_{h}$.

We consider the HDG method that approximates the scalar variable $u$, flux $\bm{q}=-\nabla u$, and boundary trace $\widehat{u}$ using the spaces $W_{h}$, $\bm{V}_{h}$, and $M_{h}$, respectively; i.e., polynomials of degree $k$ are used for all variables. We call this specific method HDGk to distinguish it from the wide variety of other available HDG methods, see, e.g., [12, 10, 11, 9]. The space $Z_{h}$ is used for postprocessing.

For the Interpolatory HDGk method, we use an elementwise interpolatory procedure along with postprocessing to approximate the nonlinear term. Let $\mathcal{I}_{h}$ be the elementwise interpolation operator with respect to the finite element nodes for the postprocessing space $Z_{h}$. Therefore, for any function $g$ that is continuous on each element we have $\mathcal{I}_{h}g\in Z_{h}$.

The Interpolatory HDGk formulation reads: find $(\bm{q}_{h},u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ such that, for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$, we have

[TABLE]

where $\Pi$ is a projection mapping into $W_{h}$ and the numerical trace for the flux is defined by

[TABLE]

Here, the stabilization function $\tau$ is nonnegative, constant on each element, and $O(1)$. Furthermore, the postprocessed scalar variable $u_{h}^{\star}=\mathfrak{q}_{h}^{k+1}(\bm{q}_{h},u_{h})\in Z_{h}$ is determined on each element $K$ by

[TABLE]

for all $(z_{h},w_{h})\in[\mathcal{P}^{k+1}(K)]^{\perp}\times\mathcal{P}^{0}(K)$, where

[TABLE]

Remark 2.1*.*

In our original Interpolatory HDG work [8], we used $\mathcal{I}_{h}^{k}F(u_{h})$ to approximate the nonlinear term, where $\mathcal{I}_{h}^{k}$ is the elementwise interpolation operator mapping into $W_{h}$. We proved optimal convergence rates for all variables, but we did not observe superconvergence after an element-by-element postprocessing. In this work, we approximate the nonlinearity using $\mathcal{I}_{h}$ and postprocessing, i.e., $\mathcal{I}_{h}F(u_{h}^{\star})$. Note that this approximate nonlinearity is in $Z_{h}$ instead of $W_{h}$ as in our first work. This simple change yields the superconvergence and keeps all the advantages of the original Interpolatory HDGk method proposed in [8]. We provide details on the computational advantages of this approach in Section 2.1.

2.1 Implementation

In our original work [8] on Interpolatory HDG, we provided details of the implementation for the method. Since we changed the discretization of the nonlinear term in this work, the implementation is different; therefore, we provide details for the implementation of this new formulation and show how all matrices need only be assembled once before the time integration. As in our earlier work [8], we describe the implementation using a simple time discretization approach: backward Euler with a Newton iteration to solve the nonlinear system at each time step. Using Interpolatory HDG with other time discretization approaches is also possible.

Let $N$ be a positive integer and define the time step $\Delta t=T/N$. We denote the approximation of $(\bm{q}_{h}(t),u_{h}(t),\widehat{u}_{h}(t))$ by $(\bm{q}^{n}_{h},u^{n}_{h},\widehat{u}^{n}_{h})$ at the discrete time $t_{n}=n\Delta t$, for $n=0,1,2,\ldots,N$. We replace the time derivative $\partial_{t}u_{h}$ in (2.5) by the backward Euler difference quotient

[TABLE]

This gives the following fully discrete method: find $(\bm{q}^{n}_{h},u^{n}_{h},\widehat{u}^{n}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ satisfying

[TABLE]

for all $(\bm{r},w,\mu)\in\bm{V}_{h}\times W_{h}\times M_{h}$ and $n=1,2,\ldots,N$. In (2.10), $f^{n}=f(t_{n},\cdot)$, the numerical trace for the flux on $\partial\mathcal{T}_{h}$ is defined by

[TABLE]

and the postprocessed approximate solution $u_{h}^{n\star}$ is determined on each element $K$ by solving

[TABLE]

for all $(z_{h},w_{h})\in[\mathcal{P}^{k+1}(K)]^{\perp}\times\mathcal{P}^{0}(K)$.

As is discussed below, the Interpolatory HDGk method takes great advantage of nodal basis functions; however, the postprocessing (2.12) uses an orthogonal complement space, which complicates the implementation. To avoid this, on each element $K$, we introduce a Lagrange multiplier $\eta_{h}^{n}\in\mathcal{P}^{0}(K)$ such that

[TABLE]

holds for all $(z_{h},w_{h})\in\mathcal{P}^{k+1}(K)\times\mathcal{P}^{0}(K)$.

Remark 2.2*.*

In this work, we used $w_{h}\in\mathcal{P}^{\ell}(K)$ with $\ell=0$ in (2.13). Actually, $\ell=0,1,\ldots,k-1$ works both in the analysis and the numerical experiments. In Part II of this work, we use the same postprocessing (2.13) with $\ell=k$.

Assume $\bm{V}_{h}=\mbox{span}\{\bm{\varphi}_{i}\}_{i=1}^{N_{1}}$, $W_{h}=\mbox{span}\{\phi_{i}\}_{i=1}^{N_{2}}$, $Z_{h}=\mbox{span}\{\chi_{i}\}_{i=1}^{N_{3}}$, and $M_{h}=\mbox{span}\{\psi_{i}\}_{i=1}^{N_{4}}$. Then

[TABLE]

Also, define the following matrices

[TABLE]

and vectors

[TABLE]

Since $\bm{V}_{h}$, $W_{h}$, and $Z_{h}$ are discontinuous finite element spaces, many of the matrices are block diagonal with small blocks.

Substitute (2.14) into the postprocessing equation (2.13) and use the corresponding test functions to test (2.13) on each element $K\in\mathcal{T}_{h}$. This gives the following local postprocessing equation

[TABLE]

were $A_{1}^{k}$ is the $k$th block of the matrix $A_{1}$, and $A_{2}^{k},b_{1}^{k}$, and $b_{2}^{k}$ are defined similarly. That is,

[TABLE]

i.e.,

[TABLE]

Let $B_{11}$ and $B_{12}$ be the block diagonal matrices with $k$th blocks $B_{11}^{k}$ and $B_{12}^{k}$, respectively.

As in [8], once we test (2.10b) using $w=\phi_{i}$ we can express the Interpolatory HDGk nonlinear term by the matrix-vector product

[TABLE]

where $\mathcal{F}$ is defined by

[TABLE]

Then the system (2.10) can be rewritten as

[TABLE]

i.e.,

[TABLE]

To apply Newton’s method to solve the nonlinear equations (2.34), define $G:\mathbb{R}^{N_{1}+N_{2}+N_{4}}\to\mathbb{R}^{N_{1}+N_{2}+N_{4}}$ by

[TABLE]

At each time step $t_{n}$ for $1\leq n\leq N$, given an initial guess $\bm{x}_{n}^{(0)}$, Newton’s method yields

[TABLE]

where the Jacobian matrix $G^{\prime}(\bm{x}_{n}^{(m-1)})$ is given by

[TABLE]

Similar to our earlier work [8] on Interpolatory HDG, the term $\mathscr{F}^{\prime}(\bm{x}_{n}^{(m-1)})$ is easily computed by

[TABLE]

where $A_{10}^{n,(m)}$ and $A_{11}^{n,(m)}$ can be efficiently computed using sparse matrix operations by

[TABLE]

Therefore, equation (2.36) can be rewritten as

[TABLE]

where

[TABLE]

This equation can be solved by locally eliminating the unknowns $\bm{\alpha}^{n,(m)}$ and $\bm{\beta}^{n,(m)}$; see [8] for details.

Remark 2.3*.*

In this new Interpolatory HDGk formulation, we only need to assemble the HDG matrices and the HDG postprocessing matrices $B_{11}$ and $B_{12}$ once before the time integration. Hence, we keep all the advantages from our earlier work [8]: the new approach eliminates the computational cost of matrix reassembly and gives simple explicit expressions for the nonlinear term and Jacobian matrix, which leads to a simple unified implementation for a variety of nonlinear PDEs.

3 Error analysis

In this section, we give a rigorous error analysis for the semidiscrete Interpolatory HDGk method. Below, we state our assumptions and briefly outline the main results. Then we provide an overview of the projections required for the analysis in Section 3.1. The proofs of the main results follow. We first assume in Section 3.2 that the nonlinearity satisfies a global Lipschitz condition. Finally, in Section 3.3 we extend the results to locally Lipschitz nonlinearities; however, we assume the mesh is quasi-uniform and $h$ is sufficiently small for this case.

We use the standard notation $W^{m,p}(\Omega)$ for Sobolev spaces on $\Omega$ with norm $\|\cdot\|_{m,p,\Omega}$ and seminorm $|\cdot|_{m,p,\Omega}$. We also write $H^{m}(\Omega)$ instead of $W^{m,2}(\Omega)$, and we omit the index $p$ in the corresponding norms and seminorms.

Throughout, we assume the solution of the PDE (1.1) exists and is unique for $t\in[0,T]$, the function $F$, the problem data, and the solution of the PDE are smooth enough, and the semidiscrete Interpolatory HDGk equations (2.5) have a unique solution on $[0,T]$. Furthermore, we assume the mesh is uniformly shape regular, $h\leq 1$, and the projection $\Pi$ used for the initial condition in (2.5d) is $\Pi=\Pi_{W}$, where $\Pi_{W}$ is defined below in Section 3.1.

We also make the following regularity assumption on the dual problem: there exists a constant $C$ such that for any $\Theta\in L^{2}(\Omega)$, the solution $(\bm{\Phi},\Psi)$ of the dual problem

[TABLE]

satisfies $(\bm{\Phi},\Psi)\in[H^{1}(\Omega)]^{d}\times H^{2}(\Omega)$ and

[TABLE]

This assumption is satisfied if $\Omega$ is convex.

We show for all $0\leq t\leq T$ the solution $(\bm{q}_{h},u_{h},u_{h}^{\star})$ of the semidiscrete Interpolatory HDGk equations (2.5) satisfies

[TABLE]

In our error estimates, the constants $C$ can vary from line to line and may depend on the exact solution and the final time $T$. As in the linear case [3], superconvergence is only obtained for $k\geq 1$.

Remark 3.1*.*

In [3], the $L^{\infty}(L^{2})$ error for $u-u_{h}^{*}$ superconverges at a rate of $\sqrt{\log\kappa}\,\,h^{k+2}$, where $\kappa$ depends on the mesh and the term $\sqrt{\log\kappa}$ grows very slowly as $h$ tends to zero. The term $\sqrt{\log\kappa}$ results from the parabolic duality argument used in [3]. It appears this parabolic duality argument is not applicable to Interpolatory HDG. Therefore, in this work we use a duality argument based on Wheeler’s work [37] and avoid the term $\sqrt{\log\kappa}$ in our error estimates; however, we require the solution has higher regularity than the regularity needed in [3] for the linear case.

3.1 Projections and basic estimates

We first introduce the HDGk projection operator $\Pi_{h}(\bm{q},u):=(\bm{\Pi}_{V}\bm{q},\Pi_{W}u)$ defined in [15], where $\bm{\Pi}_{V}\bm{q}$ and $\Pi_{W}u$ denote components of the projection of $\bm{q}$ and $u$ into $\bm{V}_{h}$ and $W_{h}$, respectively. For each element $K\in\mathcal{T}_{h}$, the projection is determined by the equations

[TABLE]

for all faces $e$ of the simplex $K$. The approximation properties of the HDGk projection (3.3) are given in the following result from [15]:

Lemma 3.2**.**

Suppose $k\geq 0$, $\tau|_{\partial K}$ is nonnegative and $\tau_{K}^{\max}:=\max\tau|_{\partial K}>0$. Then the system (3.3) is uniquely solvable for $\bm{\Pi}_{V}\bm{q}$ and $\Pi_{W}u$. Furthermore, there is a constant $C$ independent of $K$ and $\tau$ such that

[TABLE]

for $\ell_{\bm{q}},\ell_{u}$ in $[0,k]$. Here $\tau_{K}^{*}:=\max\tau|_{{\partial K}\backslash e^{*}}$, where $e^{*}$ is a face of $K$ at which $\tau|_{\partial K}$ is maximum.

Next, for each simplex $K$ in $\mathcal{T}_{h}$ and each boundary face $e$ of $K$, let $\Pi_{\ell}$ (for any $\ell\geq 0$) and $P_{M}$ denote the standard $L^{2}$ orthogonal projection operators $\Pi_{\ell}:L^{2}(K)\to\mathcal{P}^{\ell}(K)$ and $P_{M}:L^{2}(e)\to\mathcal{P}^{k}(e)$ satisfying

[TABLE]

The following error estimates for the $L^{2}$ projections and the elementwise interpolation operator $\mathcal{I}_{h}$ from Section 2 are standard and can be found in [1]:

Lemma 3.3**.**

Suppose $k,\ell\geq 0$. There exists a constant $C$ independent of $K\in\mathcal{T}_{h}$ such that

[TABLE]

3.2 Error analysis under a global Lipschitz condition

In this section, we assume the nonlinearity $F$ is globally Lipschitz:

Assumption 3.4**.**

There is a constant $L>0$ such that

[TABLE]

for all $u,v\in\mathbb{R}$.

We remove this restriction in the next section. Our proof relies on techniques used in [8, 3]. We split the proof of the main result into several steps.

To begin, we first rewrite the semidiscrete interpolatory HDG equations (2.5). First, subtract (2.5c) from (2.5b) and integrate by parts to give the following formulation:

Lemma 3.5**.**

The Interpolatory HDGk method finds $(\bm{q}_{h},u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ satisfying

[TABLE]

for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

We also define the HDGk operator $\mathscr{B}$:

[TABLE]

This allows us to rewrite the semidiscrete Interpolatory HDGk formulation (3.7) as follows: find $(\bm{q}_{h},u_{h},\widehat{u}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ such that

[TABLE]

for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

3.2.1 Step 1: Equations for the projection of the errors

Lemma 3.6**.**

For $\varepsilon_{h}^{\bm{q}}=\bm{\Pi}_{V}\bm{q}-\bm{q}_{h}$, $\varepsilon_{h}^{u}=\Pi_{W}u-u_{h}$, and $\varepsilon_{h}^{\widehat{u}}=P_{M}u-\widehat{u}_{h}$, we have

[TABLE]

for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

Proof.

By the definition of the operator $\mathscr{B}$ in (3.8), we have

[TABLE]

where we used the HDGk projection (3.3) and the $L^{2}$ projection $P_{M}$ (3.5b). Use (3.9) and subtract to obtain the result. ∎

3.2.2 Step 2: Estimate of $\varepsilon_{h}^{u}$ in $L^{\infty}(L^{2})$ by an energy argument

Lemma 3.7**.**

For any $t\in[0,T]$, we have

[TABLE]

where $\delta_{k0}$ denotes the Kronecker delta symbol so that $\delta_{k0}=1$ for $k=0$ and $\delta_{k0}=0$ for $k\geq 1$.

Proof.

We begin the proof with the case $k\geq 1$. The proof is very similar to a proof in [33], but we include it for completeness. By the properties of $\Pi_{W}$ and $\Pi_{k+1}$, we obtain

[TABLE]

Hence, for all $w_{0}\in\mathcal{P}^{0}(K)$, we have

[TABLE]

Let $e_{h}=u_{h}^{\star}-u_{h}+\Pi_{W}u-\Pi_{k+1}u$. Using the postprocessing equation (2.7), $\bm{q}=-\nabla u$, and an inverse inequality gives

[TABLE]

Since $(e_{h},1)_{K}=0$, apply the Poincaré inequality and the above estimate (3.11) to give

[TABLE]

Next, estimate the last term using an inverse inequality:

[TABLE]

This implies

[TABLE]

Hence, we have

[TABLE]

This completes the proof for the case $k\geq 1$.

For the case $k=0$, we follow the above steps in the proof of the $k\geq 1$ case but we replace the projection $\Pi_{W}$ with $\Pi_{k}$ to obtain

[TABLE]

The optimality of the $L^{2}$ projection gives $\|\Pi_{k}u-u\|_{\mathcal{T}_{h}}\leq\|\Pi_{W}u-u\|_{\mathcal{T}_{h}}$, and this completes the proof. ∎

To bound the error in the nonlinear term, we split $F(u)-\mathcal{I}_{h}F(u_{h}^{\star})$ as

[TABLE]

A bound for the first term $R_{1}$ follows directly from the standard FE interpolation error estimate (3.6a) in Lemma 3.3 due to the smoothness assumption for the function $F$. Error bounds for $R_{2}$ and $R_{3}$ are given in the following result:

Lemma 3.8**.**

We have

[TABLE]

The proofs of the estimates in Lemma 3.8 are similar to proofs in [8] and also use $\|\Pi_{k+1}u-u\|_{\mathcal{T}_{h}}\leq\|u-\mathcal{I}_{h}u\|_{\mathcal{T}_{h}}$. We omit the details.

Lemma 3.9**.**

We have the estimate

[TABLE]

where

[TABLE]

Proof.

Take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(\varepsilon_{h}^{\bm{q}},\varepsilon_{h}^{u},\varepsilon_{h}^{\widehat{u}})$ in the error equation (3.10) to give

[TABLE]

Apply the Cauchy-Schwarz inequality to each term of the right-hand side of the above identity and use $h\leq 1$, Lemma 3.7, and Lemma 3.8 to get

[TABLE]

Gronwall’s inequality, $\varepsilon_{h}^{u}(0)=0$, and $e^{Ct}\leq e^{CT}$ give the result. ∎

3.2.3 Step 3: Estimate of $\varepsilon_{h}^{q}$ in $L^{\infty}(L^{2})$ by an energy argument

Lemma 3.10**.**

We have

[TABLE]

where

[TABLE]

Proof.

Take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(\bm{r}_{h},0,0)$ in the error equation (3.10) and differentiate the result with respect to time. Also take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(0,v_{h},\widehat{v}_{h})$ in (3.10) to get

[TABLE]

for all $(\bm{r}_{h},w_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

Next, take $\bm{r}_{h}=\varepsilon_{h}^{\bm{q}}$ in (3.15a), $v_{h}=\partial_{t}\varepsilon_{h}^{u}$ in (3.15b), and $\widehat{v}_{h}=\partial_{t}\varepsilon_{h}^{\widehat{u}}$ in (3.15b) to obtain

[TABLE]

Integrating in time gives

[TABLE]

Use the Cauchy-Schwarz inequality, $h\leq 1$, Lemma 3.7, and Lemma 3.8 to get

[TABLE]

Apply the integral Gronwall’s inequality to obtain

[TABLE]

Next, let $t=0$ in (3.13) and use $\varepsilon_{h}^{u}(0)=0$ to get

[TABLE]

Therefore,

[TABLE]

and the estimate for $\|\varepsilon_{h}^{u}\|^{2}_{\mathcal{T}_{h}}$ in Lemma 3.9 completes the proof. ∎

3.2.4 Step 4: Superconvergent estimate for $\varepsilon_{h}^{u}$ in $L^{\infty}(L^{2})$ by a duality argument

To get a superconvergent rate for $\|\varepsilon_{h}^{u}\|_{\mathcal{T}_{h}}$, we adopt a duality argument from Wheeler [37]. In that work, an elliptic projection is used and it commutes with the time derivative. It is easy to check that the operator $\Pi_{W}$ defined in (3.3) commutes with the time derivative, i.e, $\partial_{t}\Pi_{W}u=\Pi_{W}u_{t}$.

For any $t\in[0,T]$, let $(\overline{\bm{q}}_{h},\overline{u}_{h},\widehat{\overline{u}}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ be the solution of the following steady state problem

[TABLE]

for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

The following estimates are proved in Section 6.

Lemma 3.11**.**

For any $t\in[0,T]$, we have

[TABLE]

Lemma 3.12**.**

Let $e_{h}^{\bm{q}}=\bm{q}_{h}-\overline{\bm{q}}_{h}$, $e_{h}^{u}=u_{h}-\overline{u}_{h}$, and $e_{h}^{\widehat{u}}=\widehat{u}_{h}-{\widehat{\overline{u}}}_{h}$. Then for any $t\in[0,T]$ we have

[TABLE]

where

[TABLE]

Proof.

By the definition of the operator $\mathscr{B}$ in (3.8), we have

[TABLE]

Take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(e_{h}^{\bm{q}},e_{h}^{u},e_{h}^{\widehat{u}})$ and use Lemma 3.7, Lemma 3.8, Lemma 3.10, the bound

[TABLE]

and also Lemma 3.11 to give

[TABLE]

Integration from [math] to $t$ gives

[TABLE]

By Gronwall’s inequality and $e_{h}^{u}(0)=u_{h}(0)-\overline{u}_{h}(0)=\Pi_{W}u_{0}-\overline{u}_{h}(0)$, we have

[TABLE]

∎

A combination of Lemma 3.11 and Lemma 3.12 gives the following lemma:

Lemma 3.13**.**

For any $t\in[0,T]$, we have

[TABLE]

where

[TABLE]

Using $u-u_{h}=(u-\Pi_{W}u)+\varepsilon_{h}^{u}$, $\bm{q}-\bm{q}_{h}=(\bm{q}-\Pi_{V}\bm{q})+\varepsilon_{h}^{\bm{q}}$, $u-u_{h}^{\star}=(u-\Pi_{k+1}u)+(\Pi_{k+1}u-u_{h}^{\star})$, and the triangle inequality gives the main result:

Theorem 3.14**.**

If the nonlinearity $F$ satisfies the global Lipschitz condition in 3.4 and the assumptions at the beginning of Section 3 hold, then for all $0\leq t\leq T$ the solution $(\bm{q}_{h},u_{h},u_{h}^{\star})$ of the semidiscrete Interpolatory HDGk equations satisfy

[TABLE]

where $\bm{\mathscr{G}}$, $\bm{\mathscr{H}}$ and $\bm{\mathscr{L}}$ are defined in (3.14), (3.12) and (3.19), respectively.

By Lemma 3.2, Lemma 3.3, and Theorem 3.14, we obtain convergence rates for smooth solutions.

Corollary 3.15**.**

If, in addition, $u$, $\bm{q}$, and $F(u)$ are sufficiently smooth for $t\in[0,T]$, then for all $0\leq t\leq T$ the solution $(\bm{q}_{h},u_{h},u_{h}^{\star})$ of the semidiscrete Interpolatory HDGk equations satisfy

[TABLE]

3.3 Error analysis under a local Lipschitz condition

In applications, the nonlinearity $F$ might not satisfy the global Lipschitz condition, see 3.4. Instead, let

[TABLE]

In this section, we assume the mesh is quasi-uniform, the polynomial degree satisfies $k\geq 1$, and the nonlinearity $F$ satisfies the following local Lipschitz condition:

Assumption 3.16**.**

There is a constant $L(M)>0$ such that

[TABLE]

for all $u,v\in[-M,M]$.

Our proof relies on techniques used in [36]. Below, we use the notation $\rho_{h}^{u^{\star}}=\Pi_{k+1}u-u_{h}^{\star}$.

Lemma 3.17**.**

*If $h$ is small enough and $k\geq 1$, then there exists $t^{*}_{h}\in(0,T]$ such that Lemma 3.10 and Lemma 3.13 hold for all $t\in[0,t^{*}_{h}]$. *

Proof.

Take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(\varepsilon_{h}^{\bm{q}},\varepsilon_{h}^{u},\varepsilon_{h}^{\widehat{u}})$ in (3.10) to give

[TABLE]

Take $t=0$ and use the fact $\varepsilon_{h}^{u}(0)=0$ to obtain

[TABLE]

By Lemma 3.7, we have

[TABLE]

The inverse inequality gives

[TABLE]

Since the exact solution is smooth at $t=0$, we can choose $h$ small enough so that $\|\rho_{h}^{u^{\star}}(0)\|_{L^{\infty}(\Omega)}<1/2$. Also, since the error equation (3.10) is continuous with respect to the time $t$, again using an inverse inequality shows that there exists $t^{*}_{h}\in(0,T]$ such that for all $h$ small enough,

[TABLE]

For all $h$ sufficiently small we have

[TABLE]

This implies for all $h$ small enough and all $t\in[0,t^{*}_{h}]$,

[TABLE]

Therefore, $u$, $\Pi_{k+1}u$, and $u_{h}^{\star}$ are located in the interval $[-M,M]$, where $M$ is defined in (3.20). This allows us to take advantage of the local Lipschitz condition of $F(u)$ for all $t\in[0,t^{*}_{h}]$. Hence, Lemma 3.10 and Lemma 3.13 hold for all $t\in[0,t^{*}_{h}]$. ∎

Lemma 3.18**.**

*For $h$ small enough and $k\geq 1$, the conclusions of Lemma 3.10 and Lemma 3.13 are true on the whole time interval $[0,T]$. *

Proof.

Fix $h^{*}>0$ so that (3.22), (3.23), and Lemma 3.17 are true for all $h\leq h^{*}$, and assume $t^{*}_{h}$ is the largest value for which (3.22) is true for all $h\leq h^{*}$. Define the set $\mathcal{A}=\{h\in[0,h^{*}]:t^{*}_{h}\neq T\}$. If the result is not true, then $\mathcal{A}$ is nonempty, $\inf\{h:h\in\mathcal{A}\}=0$, and also

[TABLE]

However, by the inverse inequality, $k\geq 1$, and since Lemma 3.17 is true, we have

[TABLE]

Since $C$ does not depend on $h$, there exists $h^{*}_{1}\leq h^{*}$ such that $\|\rho_{h}^{u^{\star}}(t^{*}_{h})\|_{L^{\infty}(\Omega)}<1/2$ for all $h\in\mathcal{A}$ such that $h\leq h^{*}_{1}$. This contradicts (3.24), and therefore $t^{*}_{h}=T$ for all $h$ small enough. ∎

Theorem 3.19**.**

If the nonlinearity $F$ satisfies the local Lipschitz condition in 3.16, the mesh is quasi-uniform, $k\geq 1$, and the assumptions at the beginning of Section 3 hold, then for all $h$ small enough the conclusions of Theorem 3.14 and Corollary 3.15 are true for all $0\leq t\leq T$.

4 Numerical Results

In this section, we present two examples to demonstrate the performance of the Interpolatory HDGk method.

Example 4.1** (The Allen-Cahn or Chaffee-Infante equation).**

We begin with an example with an exact solution in order to illustrate the convergence theory. The domain is the unit square $\Omega=[0,1]\times[0,1]\subset\mathbb{R}^{2}$, the nonlinear term is $F(u)=u^{3}-u$, and the source term $f$ is chosen so that the exact solution is $u=\sin(t)\sin(\pi x)\sin(\pi y)$. Backward Euler and Crank-Nicolson are applied for the time discretization when $k=0$ and $k=1$, respectively, where $k$ is the degree of the polynomial. The time step is chosen as $\Delta t=h$ when $k=0$ and $\Delta t=h^{2}$ when $k=1$. We report the errors at the final time $T=1$ for polynomial degrees $k=0$ and $k=1$ in Table 1. The observed convergence rates match the theory.

Example 4.2** (The Schnakenberg model).**

Next, we consider a more complicated example of a reaction diffusion PDE system with zero Neumann boundary conditions that does not satisfy the assumptions for the convergence theory established here. We consider such an example to demonstrate the applicability of the Interpolatory HDGk method to more general problems.

Specifically, we consider the Schnakenberg model, which has been used to model the spatial distribution of a morphogen; see [42] for more details. The Schnakenberg system has the form

[TABLE]

with initial conditions

[TABLE]

and homogeneous Neumann boundary conditions. The parameter values are $\kappa=100$, $a=0.1305$, $b=0.7695$, $D_{1}=0.05$, and $D_{2}=1$. We choose polynomial degree $k=1$ and apply Crank-Nicolson for the time discretization with time step $\Delta t=0.001$.

We vary the spatial domain, but keep all of parameters in the model unchanged. The first domain is the unit square $\Omega=[0,1]\times[0,1]$, and the domain is partitioned into $2048$ elements. The second domain is the circle $\Omega=\{(x,y):(x-0.5)^{2}+(y-0.5)^{2}<0.5^{2}\}$ and we use $7168$ elements.

Numerical results are shown in Figure 1–Figure 2. Spot patterns form on the square and circular domains. Our numerical results are very similar to results reported in [42].

5 Conclusion

In our earlier work [8], we considered an Interpolatory HDGk methods for semilinear parabolic PDEs with a general nonlinearity of the form $F(\nabla u,u)$. The interpolatory approach achieves optimal convergence rates and reduces the computational cost compared to standard HDG since all of the HDG matrices are assembled once before the time stepping procedure. However, the method does not have superconvergence by postprocessing.

In this work, we proposed a new superconvergent Interpolatory HDGk method for approximating the solution of reaction diffusion PDEs. Unlike our earlier Interpolatory HDGk work [8], the new method uses a postprocessing $u_{h}^{\star}$ to evaluate the nonlinear term. This change provides the superconvergence, and the new method also keeps all of the computational advantages of using an interpolatory approach for the nonlinear term. We proved the superconvergence under a global Lipschitz condition for the nonlinearity, and then extended the superconvergence results to a local Lipschitz condition assuming the mesh is quasi-uniform.

In the second part of this work [5], we again consider reaction diffusion equations and extend the ideas here to derive other superconvergent interpolatory HDG methods inspired by hybrid high-order methods [9]. However, it is currently not clear whether the present approach can be used to obtain the superconvergence for semilinear PDEs with a general nonlinearity $F(\nabla u,u)$. We are currently exploring this issue.

6 Appendix A

Recall the steady state problem (3.16) from Section 3.2.4, which we repeat here for convenience: let $(\overline{\bm{q}}_{h},\overline{u}_{h},\widehat{\overline{u}}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ be the solution of

[TABLE]

for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$. Since $\Pi_{W}$ commutes with the time derivative, taking the partial derivative of (6.1) with respect to $t$ shows $(\partial_{t}\overline{\bm{q}}_{h},\partial_{t}\overline{u}_{h},\partial_{t}\widehat{\overline{u}}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$ is the solution of

[TABLE]

for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

The proof of the following lemma is very similar to a proof in [15], hence we omit it here.

Lemma 6.1**.**

For $\varepsilon_{h}^{\overline{\bm{q}}}=\bm{\Pi}_{V}\bm{q}-\overline{\bm{q}}_{h}$, $\varepsilon_{h}^{\overline{u}}=\Pi_{W}{u}-\overline{u}_{h}$, and $\varepsilon_{h}^{\widehat{\overline{u}}}=P_{M}u-\widehat{\overline{u}}_{h}$, we have

[TABLE]

for all $(\bm{r}_{h},v_{h},\widehat{v}_{h})\in\bm{V}_{h}\times W_{h}\times M_{h}$.

The next step is the consideration of the dual problem (3.1), which we again repeat for convenience: Let

[TABLE]

By the assumption at the beginning of Section 3, this boundary value problem admits the regularity estimate

[TABLE]

for all $\Theta\in L^{2}(\Omega)$.

Lemma 6.2**.**

We have

[TABLE]

Proof.

Let $\Theta=\varepsilon_{h}^{\overline{u}}$ in the dual problem (6.4), and take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(-\bm{\Pi}_{V}\bm{\Phi},\Pi_{W}\Psi,P_{M}\Psi)$ in the definition of $\mathscr{B}$ (3.8) to get

[TABLE]

On the other hand, take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(-\bm{\Pi}_{V}\bm{\Phi},\Pi_{W}\Psi,P_{M}\Psi)$ in (6.3) to get

[TABLE]

Comparing the above two equalities gives

[TABLE]

Hence, by the regularity of the dual PDE (6.5), we have

[TABLE]

Next, take $(\bm{r}_{h},v_{h},\widehat{v}_{h})=(\varepsilon_{h}^{\overline{\bm{q}}},\varepsilon_{h}^{\overline{u}},\varepsilon_{h}^{\widehat{\overline{u}}})$ in (6.3) to obtain

[TABLE]

This implies

[TABLE]

Next, use $h\leq 1$ and substitute (6.8) into (6.7) to yield the result. ∎

Following the same steps, we obtain the following result:

Lemma 6.3**.**

We have

[TABLE]

Acknowledgements

G. Chen thanks Missouri University of Science and Technology for hosting him as a visiting scholar; some of this work was completed during his research visit. J. Singler and Y. Zhang were supported in part by National Science Foundation grant DMS-1217122. J. Singler and Y. Zhang thank the IMA for funding research visits, during which some of this work was completed.