Staggered discontinuous Galerkin methods for the Helmholtz equations with large wave number

TL;DR

This paper introduces a flexible staggered discontinuous Galerkin method for solving Helmholtz equations with large wave numbers on complex meshes, demonstrating stability, convergence, and effective handling of singular solutions.

Contribution

It presents a novel, flux-free, and mesh-robust staggered discontinuous Galerkin approach for Helmholtz equations with large wave numbers, accommodating rough and distorted grids.

Findings

Method is stable when κh is small.

Error estimates are established in L2 norm.

Numerical experiments confirm theoretical results.

Abstract

In this paper we investigate staggered discontinuous Galerkin method for the Helmholtz equation with large wave number on general quadrilateral and polygonal meshes. The method is highly flexible by allowing rough grids such as the trapezoidal grids and highly distorted grids, and at the same time, is numerical flux free. Furthermore, it allows hanging nodes, which can be simply treated as additional vertices. By exploiting a modified duality argument, the stability and convergence can be proved under the condition that $\kappa h$ is sufficiently small, where $\kappa$ is the wave number and $h$ is the mesh size. Error estimates for both the scalar and vector variables in $L^2$ norm are established. Several numerical experiments are tested to verify our theoretical results and to present the capability of our method for capturing singular solutions.

Tables

Table 1. Table 1: Convergence for ξ = 1 𝜉 1 \xi=1 by using P 1 superscript 𝑃 1 P^{1} and P 2 superscript 𝑃 2 P^{2} approximations.

$P^1$					$P^2$
$h$	${\Vert u-u_h\Vert }_0$	order	${\Vert 𝒑-{𝒑}_h\Vert }_0$	order	${\Vert u-u_h\Vert }_0$	order	${\Vert 𝒑-{𝒑}_h\Vert }_0$	order
$0.5000$	1.72E-1	–	1.65E-1	–	4.19E-2	–	4.78E-2	–
$0.2500$	3.91E-2	2.1358	4.43E-2	1.9018	5.60E-3	2.9069	5.30E-3	3.1641
$0.1250$	1.07E-2	1.8736	1.18E-2	1.9047	7.46E-4	2.9046	6.17E-4	3.1105
$0.0625$	2.70E-3	1.9730	3.00E-3	1.9702	9.49E-5	2.9743	7.42E-5	3.0555
$0.0313$	6.82E-4	1.9933	7.58E-4	1.9920	1.19E-5	2.9937	9.13E-6	3.0242
$0.0156$	1.70E-4	1.9983	1.90E-4	1.9979	1.49E-6	2.9984	1.13E-6	3.0107

Table 2. Table 2: Convergence for ξ = 3 / 2 𝜉 3 2 \xi=3/2 by using P 1 superscript 𝑃 1 P^{1} and P 2 superscript 𝑃 2 P^{2} approximations.

$P^1$					$P^2$
$h$	${\Vert u-u_h\Vert }_0$	order	${\Vert 𝒑-{𝒑}_h\Vert }_0$	order	${\Vert u-u_h\Vert }_0$	order	${\Vert 𝒑-{𝒑}_h\Vert }_0$	order
$0.5000$	1355E-1	–	2.23E-1	–	4.78E-2	–	4.87E-2	–
$0.2500$	3.44E-2	1.9794	5.18E-2	2.1085	6.10E-3	2.9783	7.30E-3	2.7415
$0.1250$	1.04E-2	1.7200	1.59E-2	1.7029	8.86E-4	2.7756	1.10E-3	2.6751
$0.0625$	3.00E-3	1.7801	4.70E-3	1.7572	1.75E-4	2.3379	3.44E-4	1.7298
$0.0313$	9.52E-4	1.6723	1.50E-3	1.6806	5.05E-5	1.7936	1.23E-4	1.4858
$0.0156$	3.22E-4	1.5649	4.89E-4	1.5839	1.70E-5	1.5683	4.44E-5	1.4685

Table 3. Table 3: Convergence for ξ = 2 / 3 𝜉 2 3 \xi=2/3 by using P 1 superscript 𝑃 1 P^{1} and P 2 superscript 𝑃 2 P^{2} approximations.

$P^1$					$P^2$
$h$	${\Vert u-u_h\Vert }_0$	order	${\Vert 𝒑-{𝒑}_h\Vert }_0$	order	${\Vert u-u_h\Vert }_0$	order	${\Vert 𝒑-{𝒑}_h\Vert }_0$	order
$0.5000$	2.33E-1	–	1.92E-1	–	4.68E-2	–	5.79E-2	–
$0.2500$	7.35E-2	1.6675	6.88E-2	1.4832	8.70E-3	2.4330	1.15E-2	2.3287
$0.1250$	3.03E-2	1.2799	2.71E-2	1.3428	2.30E-3	1.8885	5.90E-3	0.9690
$0.0625$	1.59E-2	0.9317	1.64E-2	0.7273	7.36E-4	1.6665	3.60E-3	0.7182
$0.0313$	9.40E-3	0.7498	1.13E-2	0.5355	2.52E-4	1.5477	2.20E-3	0.6817
$0.0156$	5.80E-3	0.6904	7.80E-3	0.5370	1.01E-4	1.3216	1.40E-3	0.6706

Equations

$$-\Delta u-{\kappa }^{2}u$$

$$\nabla u\cdot \mathbf{n}+i\kappa u$$

$$\begin{array}{rl}i\kappa \mathbf{p}& =-\nabla u\text{mbox}in\Omega ,\\ \nabla \cdot \mathbf{p}+i\kappa u& =f\text{mbox}in\Omega ,\\ -\mathbf{p}\cdot \mathbf{n}+u& =g\text{mbox}on\partial \Omega .\end{array}$$

$$\bar{\Omega }=\bigcup _{T\in {\mathcal{T}}_u}T.$$

$$\bar{\Omega }=\bigcup _{\tau \in {\mathcal{T}}_h}\bar{\tau }.$$

$$[w]=w_1-w_2,[\mathbf{v}\cdot \mathbf{n}]={\mathbf{v}}_1\cdot \mathbf{n}-{\mathbf{v}}_2\cdot \mathbf{n},$$

$$S_h:=\{w:w{\mid }_{\tau }\in P^m(\tau )\forall \tau \in {\mathcal{T}}_h;[w]{\mid }_e=0\forall e\in {\mathcal{F}}_u^0\}.$$

$${\mathbf{V}}_h=\{\mathbf{v}:\mathbf{v}{\mid }_{\tau }\in P^m(\tau {)}^2\forall \tau \in {\mathcal{T}}_h;[\mathbf{v}\cdot \mathbf{n}]{\mid }_e=0\forall e\in {\mathcal{F}}_p\}.$$

$$(i\kappa {\mathbf{p}}_h,{\overline{\mathbf{q}}}_h)$$

$$b_h({\mathbf{p}}_h,{\overline{v}}_h)+i\kappa (u_h,{\overline{v}}_h)+⟨u_h,{\overline{v}}_h{⟩}_{\partial \Omega }$$

$$b_h({\mathbf{q}}_h,v_h)$$

$$b_h^{\ast }(v_h,{\mathbf{q}}_h)$$

$$b_h({\mathbf{q}}_h,{\overline{v}}_h)=b_h^{\ast }({\overline{v}}_h,{\mathbf{q}}_h)\forall (v_h,{\mathbf{q}}_h)\in S_h\times {\mathbf{V}}_h.$$

$$A_h(\mathbf{p}-{\mathbf{p}}_h,u-u_h;{\mathbf{q}}_h,v_h)+i\kappa (u-u_h,{\overline{v}}_h)+⟨u-u_h,{\overline{v}}_h{⟩}_{\partial \Omega }=0\forall ({\mathbf{q}}_h,v_h)\in {\mathbf{V}}_h\times S_h.$$

$$⟨I_{h}v-v,\overline{\varphi }{⟩}_e$$

$$(I_{h}v-v,\overline{\varphi }{)}_{\tau }$$

$$⟨(J_h\mathbf{q}-\mathbf{q})\cdot \mathbf{n},\overline{v}{⟩}_e$$

$$(J_h\mathbf{q}-\mathbf{q},\overline{\mathbf{\varphi }}{)}_{\tau }$$

$$b_h(\mathbf{p}-J_h\mathbf{p},\overline{v})$$

$$b_h^{\ast }(u-I_{h}u,\overline{\mathbf{q}})$$

$$\begin{array}{rl}& A_h(J_h\mathbf{p}-{\mathbf{p}}_h,I_{h}u-u_h;{\mathbf{q}}_h,v_h)+i\kappa (I_{h}u-u_h,{\overline{v}}_h)+⟨I_{h}u-u_h,{\overline{v}}_h{⟩}_{\partial \Omega }\\ & =A_h(J_h\mathbf{p}-\mathbf{p},I_{h}u-u;{\mathbf{q}}_h,v_h)+i\kappa (I_{h}u-u,{\overline{v}}_h)+⟨I_{h}u-u,{\overline{v}}_h{⟩}_{\partial \Omega }\\ & =-i\kappa (\overline{J_h\mathbf{p}-\mathbf{p}},{\mathbf{q}}_h)+i\kappa (I_{h}u-u,{\overline{v}}_h)\forall ({\mathbf{q}}_h,v_h)\in {\mathbf{V}}_h\times S_h.\end{array}$$

$$|u_{ϵ}{|}_j$$

$$|u_{\mathcal{A}}{|}_j$$

$$\Vert u-I_{h}u{\Vert }_0+h\Vert \nabla (u-I_{h}u){\Vert }_0$$

$$\Vert \mathbf{p}-J_h\mathbf{p}{\Vert }_0$$

$$\kappa \Vert \mathbf{p}-J_h\mathbf{p}{\Vert }_0$$

$$\Vert u-I_{h}u{\Vert }_0$$

$$\Vert \nabla (u-I_{h}u){\Vert }_0$$

$$i\kappa \mathbf{p}$$

$$\nabla \cdot \mathbf{p}$$

$$-\mathbf{p}\cdot \mathbf{n}+u$$

$$i\kappa ({\mathbf{p}}_h^{+},{\overline{\mathbf{q}}}_h)$$

$$b_h({\mathbf{p}}_h^{+},{\overline{v}}_h)+⟨u_h^{+},{\overline{v}}_h{⟩}_{\partial \Omega }$$

$$B_h({\mathbf{p}}_h^{+},u_h^{+};{\mathbf{q}}_h,v_h)+⟨u_h^{+},{\overline{v}}_h{⟩}_{\partial \Omega }=B_h(\mathbf{p},u;{\mathbf{q}}_h,v_h)+⟨u,{\overline{v}}_h{⟩}_{\partial \Omega }\forall ({\mathbf{q}}_h,v_h)\in {\mathbf{V}}_h\times S_h,$$

$$i\kappa (\mathbf{p}-{\mathbf{p}}_h^{+},{\overline{\mathbf{q}}}_h)$$

$$b_h(\mathbf{p}-{\mathbf{p}}_h^{+},{\overline{v}}_h)+⟨I_{h}u-u_h^{+},{\overline{v}}_h{⟩}_{\partial \Omega }$$

$$B_h({\mathbf{q}}_h,v_h;{\mathbf{p}}_h^{-},u_h^{-})+⟨v_h,\overline{u_h^{-}}{⟩}_{\partial \Omega }=B_h({\mathbf{q}}_h,v_h;\mathbf{p},u)+⟨v_h,\overline{u}{⟩}_{\partial \Omega }\forall ({\mathbf{q}}_h,v_h)\in {\mathbf{V}}_h\times S_h.$$

Full text

Staggered discontinuous Galerkin methods for

the Helmholtz equations with large wave number

Lina Zhao111Department of Mathematics,The Chinese University of Hong Kong, Hong Kong Special Administrative Region. (lzhao@math.cuhk.edu.hk)Department of Computational Science and Engineering, Yonsei University, Seoul 03722, Republic of Korea. (ejpark@yonsei.ac.kr)   Eun-Jae Park11footnotemark: 1   Eric T. Chung333Department of Mathematics,The Chinese University of Hong Kong, Hong Kong Special Administrative Region. (tschung@math.cuhk.edu.hk)

Abstract: In this paper we investigate staggered discontinuous Galerkin method for the Helmholtz equation with large wave number on general quadrilateral and polygonal meshes. The method is highly flexible by allowing rough grids such as the trapezoidal grids and highly distorted grids, and at the same time, is numerical flux free. Furthermore, it allows hanging nodes, which can be simply treated as additional vertices. By exploiting a modified duality argument, the stability and convergence can be proved under the condition that $\kappa h$ is sufficiently small, where $\kappa$ is the wave number and $h$ is the mesh size. Error estimates for both the scalar and vector variables in $L^{2}$ norm are established. Several numerical experiments are tested to verify our theoretical results and to present the capability of our method for capturing singular solutions.

Keywords: Helmholtz problem, Large wave number, Staggered DG method, Duality argument, General quadrilateral and polygonal meshes

1 Introduction

In this paper we develop a staggered discontinuous Galerkin (DG) method for solving the following Helmholtz problem

[TABLE]

where $\kappa>0$ is the wave number and $f\in L^{2}(\Omega)$ represents a harmonic source and $g\in L^{2}(\partial\Omega)$ is a given data function. Here, $\Omega$ is a polygonal or polyhedral domain in $\mathbb{R}^{2}$ and $i=\sqrt{-1}$ is the imaginary unit.

The Helmholtz problem has many practical applications in electrodynamics, especially in optics and in acoustic involving time harmonic wave propagation. The Helmholtz equation with large wave number is indefinite, which makes it difficult to design robust and accurate numerical methods for the Helmholtz problem. For a fixed polynomial order, the pollution effect can be reduced substantially but cannot be avoided in principle [4]. A rigorous analysis for one-dimensional Helmholtz problems and piecewise linear approximation has been given in [2] and [21], respectively, under the condition that $\kappa^{2}h$ is sufficiently small. However, the assumption on $\kappa^{2}h$ is too restrictive and unsatisfactory from a practical point of view. In order to obtain more stable and accurate numerical approximations, a large amount of nonstandard methods have been proposed and analyzed. Among all the methods, we mention in particular the general DG method on regular meshes [34] with mesh condition $\kappa(\kappa h)^{m}\leq C_{0}$, the absolutely stable discontinuous Galerkin (DG) methods [22, 23, 24] without any mesh constraint, continuous interior penalty finite element methods (CIP-FEM) [44] with the mesh condition $\frac{\kappa h}{m}\leq C_{0}(\frac{m}{\kappa})^{\frac{1}{m+1}}$ and a new weak Galerkin (WG) finite element method [38] with the mesh condition $\kappa^{7/2}h^{2}\leq C_{0}$ or $(\kappa h)^{2}+\kappa(\kappa h)^{m+1}\leq C_{0}$, where $m$ is the polynomial order. In addition to the above mentioned approaches, many other numerical methods have also been developed, such as the partition of unity finite element methods [3, 33], the least squares finite element methods [9, 27, 28, 37, 11], the generalized finite element methods [6, 5], the hybridized discontinuous Galerkin methods [26], the interior penalty discontinuous Galerkin methods [1] and the Petrov Galerkin methods [20, 25].

Staggered DG method is initially developed for wave propagation problems [13, 14] on triangular meshes. Since then, it have been successfully applied to a large amount of partial differential equations arising from practical applications, see, e.g., [18, 15, 16, 31, 12, 30, 32, 17, 39, 41] and the references therein. Recently, staggered DG methods have been designed on general quadrilateral and polygonal meshes to solve Darcy law and the Stokes equations [40, 42]. The key features of staggered DG methods can be summarized as follows: First, it can preserve the physical properties, such as the local and global mass conservation, and achieve the superconvergent estimates. Second, it can be flexibly applied to rough grids such as the highly distorted grids and polygonal grids, and at the same time hanging nodes can be simply incorporated into the method. Third, thanks to the staggered continuity property, no numerical flux is needed in the construction of the method. All these distinctive properties make staggered DG method competitive in real applications.

The goal of this paper is to extend staggered DG method on general quadrilateral and polygonal meshes to the Helmholtz problem with large wave number. To the best of our knowledge, very few results are available for the Helmholtz problem on general meshes in the existing literature (cf. [38]). The key idea for staggered DG method is to divide the initial partion (quadrilateral or polygonal meshes) into the union of triangles, then the primal mesh, the dual mesh and the primal simplexes can be constructed. Next, two sets of basis functions for the scalar and vector variables, respectively with staggered continuity for the Helmholtz problem are defined. The primary difficulty of analyzing the Helmholtz problem lies in the strong indefiniteness of the problem which makes it hard to establish the stability for the numerical approximation. To analyze staggered DG method, we exploit a modified duality argument. In addition, the elliptic projections in the spirit of staggered DG method for Darcy problem are defined. The key idea employed here is to use the specially designed elliptic projections in the duality argument to bound the $L^{2}$ errors of the discrete solution by the elliptic projections under the condition that $\kappa h$ is small enough. It is worth mentioning that the superconvergence of the elliptic projections for the scalar variable is the crux to achieve the optimal convergence for both the scalar and vector variables in $L^{2}$ errors with explicit dependence on $\kappa$. In addition, our estimates are comparable to those given in [44] for CIP-FEM. Note that our (undisplayed) analysis shows that if we apply traditional duality argument as proposed in [10], then the stability estimates require that $\kappa^{2}h$ should be sufficiently small, but this mesh condition is too restrictive for large wave number $\kappa$. Thus, we turn to the modified duality argument and improve the mesh condition, namely, we only require $\kappa h$ is small enough. It is worth mentioning that this is the first result on staggered DG method for the Helmholtz problems.

The rest of the paper is organized as follows. In the next section, we present the construction of staggered DG method for the Helmholtz problem. Then in Section 3 we analyze the convergence of staggered DG method, where a modified duality argument is exploited. Finally, some numerical experiments are carried out in Section 4 to confirm the proposed theories.

2 Staggered DG method

The staggered DG method is based on a first order formulation of the Helmholtz problem (1.1)-(1.2), which can be written in mixed form as finding $(\bm{p},u)$ such that

[TABLE]

Next, we will briefly introduce staggered DG method for the Helmholtz problem in mixed form (cf. (2.1)) on general quadrilateral and polygonal meshes, and more details can be referred to [40, 42]. To begin, we construct three meshes: the primal mesh $\mathcal{T}_{u}$, the dual mesh $\mathcal{T}_{p}$, and the primal simplexes $\mathcal{T}_{h}$. For a polygonal domain $\Omega$, consider a general mesh $\mathcal{T}_{u}$ (of $\Omega$) that consists of nonempty connected close disjoint subsets of $\Omega$ (see Figure 1):

[TABLE]

We also let $\mathcal{F}_{u}$ be the set of all primal edges in this partition and $\mathcal{F}_{u}^{0}$ be the subset of all interior edges, that is, the set of edges in $\mathcal{F}_{u}$ that do not lie on $\partial\Omega$. We construct the primal submeshes $\mathcal{T}_{h}$ as a triangular subgrid of the primal grid: for an element $T\in\mathcal{T}_{u}$, elements of $\mathcal{T}_{h}$ are obtained by connecting the interior point $\nu$ to all vertices of $\mathcal{T}_{u}$ (see Figure 1):

[TABLE]

We rename the union of these triangles by $S(\nu)$. Moreover, we will use $\mathcal{F}_{p}$ to denote the set of all the dual edges generated by this subdivision process. For each triangle $\tau\in\mathcal{T}_{h}$, we let $h_{\tau}$ be the diameter of $\tau$ and $h=\max\{h_{\tau},\tau\in\mathcal{T}_{h}\}$. In addition, we define $\mathcal{F}:=\mathcal{F}_{u}\cup\mathcal{F}_{p}$ and $\mathcal{F}^{0}:=\mathcal{F}_{u}^{0}\cup\mathcal{F}_{p}$. The construction for general meshes is illustrated in Figure 1, where the black solid lines are edges in $\mathcal{F}_{u}$ and the red dotted lines are edges in $\mathcal{F}_{p}$.

Finally, we construct the dual mesh. For each interior edge $e\in\mathcal{F}_{u}^{0}$, we use $D(e)$ to denote the dual mesh, which is the union of the two triangles in $\mathcal{T}_{h}$ sharing the edge $e$, and for each boundary edge $e\in\mathcal{F}_{u}\backslash\mathcal{F}_{u}^{0}$, we use $D(e)$ to denote the triangle in $\mathcal{T}_{h}$ having the edge $e$, see Figure 1. We write $\mathcal{T}_{p}$ as the union of all $D(e)$.

For each edge $e$, we define a unit normal vector $n_{e}$ as follows: If $e\in\mathcal{F}\setminus\mathcal{F}^{0}$, then $\bm{n}_{e}$ is the unit normal vector of $e$ pointing towards the outside of $\Omega$. If $e\in\mathcal{F}^{0}$, an interior edge, we then fix $\bm{n}_{e}$ as one of the two possible unit normal vectors on $e$. When there is no ambiguity, we use $\bm{n}$ instead of $\bm{n}_{e}$ to simplify the notation.

Let $D\subset\mathbb{R}^{2}$, we adopt the standard notations for the Sobolev spaces $H^{s}(D)$ and their associated norms $\|\cdot\|_{s,D}$, and semi-norms $|\cdot|_{s,D}$ for $s\geq 0$. In particular, $(\cdot,\cdot)_{D}$ and $\langle\cdot,\cdot\rangle_{\Sigma}$ for $\Sigma\subset D$ denote the $L^{2}$ inner product on complex valued $L^{2}(D)$ and $L^{2}(\Sigma)$, respectively. If $D=\Omega$, the subscript $\Omega$ will be dropped unless otherwise mentioned. In the sequel, we use $C$ to denote a generic positive constant which may have different values at different occurrences.

The following mesh regularity assumptions are also needed throughout the paper (cf. [7, 8]).

Assumption 2.1**.**

We assume there exists a constant $\rho>0$ such that

1. (1)

For every element $S(\nu)\in\mathcal{T}_{u}$ and every edge $e\in\partial S(\nu)$, it satisfies $h_{e}\geq\rho h_{S(\nu)}$, where $h_{e}$ denotes the length of edge $e$ and $h_{S(\nu)}$ denotes the diameter of $S(\nu)$. 2. (2)

Every element $S(\nu)$ in $\mathcal{T}_{u}$ is star-shaped with respect to a ball of radius $\geq\rho h_{S(\nu)}$.

We remark that the above assumptions ensure that the triangulation $\mathcal{T}_{h}$ is shape regular.

Let $m\geq 0$ be the order of approximation. For every $\tau\in\mathcal{T}_{h}$ and $e\in\mathcal{F}$, we define $P^{m}(\tau)$ and $P^{m}(e)$ as the spaces of polynomials of degree less than or equal to $m$ on $\tau$ and $e$, respectively. For $w$ and $\bm{v}$ belonging to the broken Sobolev space, the jump $[v]$ and the jump $[\bm{v}\cdot\bm{n}]$ are defined respectively as

[TABLE]

where $v_{i}=v\mid_{\tau_{i}}$, $\bm{v}_{i}=\bm{v}\mid_{\tau_{i}}$ and $\tau_{1}$, $\tau_{2}$ are the two triangles in $\mathcal{T}_{h}$ having the edge $e\in\mathcal{F}$. In the above definitions, we assume $\bm{n}$ is pointing from $\tau_{1}$ to $\tau_{2}$.

Next, we will introduce some finite dimensional spaces. First, we define the following locally $H^{1}(\Omega)$ conforming space $S_{h}$:

[TABLE]

Notice that, if $w\in S_{h}$, then $w\mid_{D(e)}\in H^{1}(D(e))$ for each edge $e\in\mathcal{F}_{u}$. We next define the following locally $H(\mbox{div};\Omega)-$conforming SDG space $\bm{V}_{h}$:

[TABLE]

Note that if $\bm{v}\in\bm{V}_{h}$, then $\bm{v}\mid_{S(\nu)}\in H(\textnormal{div};S(\nu))$ for each $S(\nu)\in\mathcal{T}_{u}$.

Following [40, 42], we can obtain the staggered DG formulation for (2.1): find $(u_{h},\bm{p}_{h})\in S_{h}\times\bm{V}_{h}$ such that

[TABLE]

where the bilinear forms are defined as

[TABLE]

The following discrete adjoint property can be verified easily by integration by parts:

[TABLE]

Let $A_{h}(\bm{p}_{h},u_{h};\bm{q}_{h},v_{h})=-i\kappa(\overline{\bm{p}}_{h},\bm{q}_{h})-b_{h}^{*}(\overline{u}_{h},\bm{q}_{h})+b_{h}(\bm{p}_{h},\overline{v}_{h})$. Then we have from integration by parts

[TABLE]

3 Convergence analysis

In this section, we aim to derive the convergence estimates and stability of staggered DG method for the Helmholtz problem, our (undisplayed) analysis shows that standard duality argument requires $\kappa^{2}h$ be sufficiently small to achieve the stability. To overcome this issue, a modified duality argument is exploited, where the elliptic projections are the key tools. From which we can get the stability and convergence estimates provided $\kappa h$ is small enough.

To begin, we define the following projection operators, which will be useful for the subsequent analysis. Let $I_{h}:H^{1}(\Omega)\rightarrow S_{h}$ be defined by

[TABLE]

In addition, $J_{h}:H^{1}(\Omega)\rightarrow\bm{V}_{h}$ is defined by

[TABLE]

By the definitions of $I_{h}$ and $J_{h}$, we can get

[TABLE]

In addition, we can also derive from the definition of the bilinear form $A_{h}(\cdot,\cdot;\cdot,\cdot)$ (cf. (2.5))

[TABLE]

The next lemma is found to be useful for the subsequent analysis, one can refer to [36] for proof.

Lemma 3.1**.**

The solution $u$ to the problem (1.1) and (1.2) can be written as $u=u_{\epsilon}+u_{\mathcal{A}}$, and satisfies

[TABLE]

where $M_{f,g}=\|f\|_{0}+\|g\|_{1/2,\partial\Omega}$.

The approximation estimates for $I_{h}$ and $J_{h}$ with high order convergence are given as follows (cf. [19, 14]).

Lemma 3.2**.**

Assume that $u\in H^{m+1}(\Omega)$, then we have

[TABLE]

For our convergence analysis, we also need the following approximation estimates, which only require $H^{2}$ regularity. It uses the decomposition given in Lemma 3.1 and is more subtle than the estimates given in Lemma 3.2. The proof is analogous to that of [35], thus we omit the proof for simplicity.

Lemma 3.3**.**

(approximation properties) Let $(\bm{p},u)$ be the solution to (2.1), then

[TABLE]

The rest of this section is devoted to the error analysis for the Helmholtz problem, to this end we introduce the elliptic projections motivated by those proposed in [44, 43]. For any $u,\bm{p}$, we define its elliptic projections $u_{h}^{+},\bm{p}_{h}^{+}$ as the staggered DG approximations to the first order system

[TABLE]

Then $(\bm{p}_{h}^{+},u_{h}^{+})\in\bm{V}_{h}\times S_{h}$ is the numerical solution to the following discrete formulation

[TABLE]

Let $B_{h}(\bm{p}_{h},u_{h};\bm{q}_{h},v_{h})=i\kappa(\bm{p}_{h},\overline{\bm{q}}_{h})-b_{h}^{*}(u_{h},\overline{\bm{q}}_{h})+b_{h}(\bm{p}_{h},\overline{v}_{h})$, then it follows

[TABLE]

which immediately yields

[TABLE]

In addition, we define $\bm{p}_{h}^{-},u_{h}^{-}$ by

[TABLE]

The following holds

[TABLE]

Indeed, we have

[TABLE]

To be specific,we can obtain from (3.7)

[TABLE]

In other words, $(\bm{p}_{h}^{-},u_{h}^{-})\in\bm{V}_{h}\times S_{h}$ is the staggered DG approximation to the following first order system

[TABLE]

On the other hand we have the following identity by integration by parts

[TABLE]

In addition, if $(\bm{p},u)$ is the solution of (2.1) then

[TABLE]

Now we are ready to prove the next lemma.

Lemma 3.4**.**

Assume $u$ is any function in $H^{2}(\Omega)$ and $\bm{p}$ is any function in $H^{1}(\Omega)^{2}$. Then it holds

[TABLE]

Proof.

The proof for the elliptic projections $(u_{h}^{+},\bm{p}_{h}^{+})$ and $(u_{h}^{-},\bm{p}_{h}^{-})$ are similar, to simplify the presentation, we only give the proof for $(u_{h}^{+},\bm{p}_{h}^{+})$.

Proceeding analogously to (3.3), we can obtain

[TABLE]

Let $\bm{q}_{h}=J_{h}\bm{p}-\bm{p}_{h}^{+},v_{h}=I_{h}u-u_{h}^{+}$, then

[TABLE]

On the other hand, we have from the definition of $A_{h}$ and the discrete adjoint property (2.4) that

[TABLE]

Taking the imaginary part of (3.10), we can obtain

[TABLE]

Thus

[TABLE]

Next, we will estimate $\|I_{h}u-u_{h}^{+}\|_{0}$. Consider the auxiliary problem

[TABLE]

which satisfies the following elliptic regularity estimate

[TABLE]

An application of the Cauchy-Schwarz inequality, the discrete adjoint property (2.4), (3.2), (3.6), (3.11) and the elliptic regularity estimate (3.12) lead to

[TABLE]

Therefore, the proof is complete.

∎

Lemmas 3.3 and 3.4 yield the next lemma.

Lemma 3.5**.**

Let $(\bm{p},u)$ be the solution to the problem (2.1). Then we have

[TABLE]

The $L^{2}$ error estimates for both $\bm{p}_{h}$ and $u_{h}$ are stated in the next lemma.

Lemma 3.6**.**

Let $(u,\bm{p})$ and $(u_{h},\bm{p}_{h})$ denote the solution of (2.1) and (2.2)-(2.3), respectively. If $\kappa h$ is sufficiently small, then the following estimates hold

[TABLE]

Proof.

We first estimate the $L^{2}$ error of $u_{h}$ by introducing the dual problem and exploiting the elliptic projections of the solution to the original continuous problem and of the solution to the dual problem. We consider the following dual problem

[TABLE]

Let $(\bm{p}_{h}^{+},u_{h}^{+})\in\bm{V}_{h}\times S_{h}$ be the elliptic projections defined by (3.4) and let $(\bm{\psi}_{h}^{-},\varphi_{h}^{-})\in\bm{V}_{h}\times S_{h}$ be the elliptic projections defined by (3.7) by replacing $(\bm{p},u)$ by $(\bm{\psi},\varphi)$. It follows from Lemma 3.5 by replacing $u$ by $\varphi$ that

[TABLE]

In addition, Lemma 3.3 yields

[TABLE]

Therefore

[TABLE]

Integration by parts, (3.4), (3.8), (3.9) and (3.14) reveal that

[TABLE]

Now we will estimate each of the above term separately. First, the definition of $B_{h}(\cdot,\cdot;\cdot,\cdot)$ yields

[TABLE]

The Cauchy-Schwarz inequality, Lemma 3.4, (3.15) and (3.16) imply

[TABLE]

We have from the Cauchy-Schwarz inequality, the inverse inequality, Lemma 3.4 and (3.15)

[TABLE]

Similarly, we can get

[TABLE]

An appeal to the Cauchy-Schwarz inequality and (3.17) yields

[TABLE]

The trace inequality, the inverse inequality, Lemma 3.4, (3.16) and (3.17) imply

[TABLE]

Combining the preceding estimates yields (3.13) under the condition that $\kappa h$ is sufficiently small.

Next, we estimate $\|\bm{p}-\bm{p}_{h}\|_{0}$. We have from (3.3)

[TABLE]

Taking the imaginary parts of the above equation implies

[TABLE]

Consequently, we have

[TABLE]

Thus, the proof is complete.

∎

Combining Lemmas 3.3 and 3.6, we can get the following estimates.

Corollary 3.1**.**

If $u\in H^{2}(\Omega)$, then

[TABLE]

Combining Lemmas 3.2 and 3.6, we have

Theorem 3.1**.**

If $(\bm{p},u)\in H^{m}(\Omega)^{2}\times H^{m+1}(\Omega)$ satisfying

[TABLE]

then

[TABLE]

By combining Lemma 3.1 and Corollary 3.1, we can obtain the following stability estimates for staggered DG method.

Corollary 3.2**.**

(stability) Suppose the solution $u\in H^{2}(\Omega)$, then the following estimate holds provided $\kappa h$ is sufficiently small

[TABLE]

4 Numerical experiments

This section presents numerical experiments for the validation of the theoretical results and investigates the accuracy of our staggered DG method. In addition, an singular example is given to test the capability of the proposed method to capture the singularities.

Example 4.1**.**

(Typical smooth solution example)

We first consider a Helmholtz equation defined on the unit square $\Omega=[-0.5,0.5]\times[-0.5,0.5]$. Here, we set $f=\frac{\sin(\kappa r)}{r}$ in (1.1)-(1.2) and $g$ is chosen such that the exact solution is given by

[TABLE]

in polar coordinates, where $J_{\nu}(z)$ are Bessel functions of the first kind.

The exact solution for $\kappa=50$ is displayed in Figure 2, and the numerical solution for $\kappa=50$ by using $P^{1}$ and $P^{3}$ is shown in Figure 3. As expected, higher order polynomial approximation yields better numerical solution.

For the fixed wave number $\kappa$,we first show the dependence of the convergence of $\|u-u_{h}\|_{0}$, $\|\bm{p}-\bm{p}_{h}\|_{0}$ on polynomial order $m$ and mesh size $h$. The convergence history against the number of degrees of freedom for $\kappa=50$ and $\kappa=100$ with different polynomial orders on square grids are reported in Figure 4. It is easy to see that the pollution errors always appear on the coarse meshes, and the optimal convergence can be obtained on fine meshes for different polynomial orders, which confirm the proposed theories.

Next we verify the convergence properties of the SDG method for different wave numbers by piecewise $P^{1},P^{2}$ and $P^{3}$ approximations, respectively. We can see from Figure 5 that in the pre-asymptotic region, the errors always oscillate for different polynomial orders, and optimal convergence can be obtained for fine meshes. In addition, high order polynomial approximation has better performances.

In addition, to test the flexibility of the proposed method on rough grids, we employ the $h$-perturbation grids (cf. [40]) as shown in Figure 6. The $L^{2}$ errors for both the scalar and vector variables by using different polynomial approximations for $\kappa=50$ and $\kappa=10$ are displayed in Figure 7. From which we can see that similar performances as for square grids can be obtained, thus, the proposed method is robust in the sense that it can be flexibly applied to rough grids. Undisplayed numerical experiments on polygonal grids also yield similar results.

Example 4.2**.**

(Singular solution example)

In this example, we consider a square domain $\Omega=(0,1)\times(-0.5,0.5)$. The boundary conditions $g$ and $f$ are chosen such that the exact solution is given by

[TABLE]

for $\xi=1$, $\xi=2/3$ and $\xi=3/2$, respectively, where $J_{\xi}$ denotes the Bessel function of the first kind and order $\xi$. It is well known that $u$ is smooth for $\xi\in\mathbb{N}$, while its derivative has a singularity at origin for $\xi\neq\mathbb{N}$ (cf. [29]). We fix $\kappa=10$, and we compute the numerical solution in the regular case $\xi=1$ and in the singular cases $\xi=1/2$ and $\xi=3/2$. The profiles of the numerical solutions corresponding to these three cases are displayed in Figures 8 and 9. The convergence order against the mesh size $h$ for these three cases by using $P^{1}$ and $P^{2}$ polynomial approximations are reported in Tables 1–3. It is easy to see that when $\xi=1$, optimal convergence can be obtained for both variables in $L^{2}$ norm using linear and quadratic approximations. when $\xi=3/2$, a reduced convergence order can be achieved for both the flux variable and the scalar variable, in addition, the orders are the same for $P^{1}$ and $P^{2}$ polynomial approximations. When $\xi=2/3$, the convergence order gets worse. From the above, we can conclude that a reduced convergence order can be obtained due to the reduced regularity.

Acknowledgements

The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund (Project numbers 14304217 and 14302018) and CUHK Faculty of Science Direct Grant 2018-19.