# A Mixed Discontinuous Galerkin Method for the Helmholtz Equation

Authors: Qingjie Hu; Yinnian He; Tingting Li; Jing Wen
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we introduce and analyze a mixed discontinuous Galerkin method for the Helmholtz equation. The mixed discontinuous Galerkin method is designed by using a discontinuousPp+1−1−Pp−1 finite element pair for the flux variable and the scattered field with p≥0. We can get optimal order convergence for the flux variable in both Hdiv-like norm and L2 norm and the scattered field in L2 norm numerically. Moreover, we conduct the numerical experiments on the Helmholtz equation with perturbation and the rectangular waveguide, which also demonstrate the good performance of the mixed discontinuous Galerkin method.


---

## Body

## 1. Introduction

We consider the following nonhomogeneous Helmholtz equation with the Robin boundary condition:(1)−Δu−k2u=f,inΩ∈Rd,(2)∂u∂nΩ−iku=−g,onΓ=∂Ω,where d=2,3, i=−1 denotes the imaginary unit, k∈R+ is a given positive number and known as the wave number, Ω is an open and bounded domain, and f∈L2Ω represents a harmonic source. The Robin boundary condition (2) is the lowest-order absorbing boundary condition [1]. The applications of the Helmholtz equation are extensive in many practical applications, such as geophysics and radar detecting, simulation of ground penetrating, biomedical imaging, acoustic noise control, and seismic wave propagation. The numerical solution of the Helmholtz equation is fundamental to the simulation of time harmonic wave phenomena in acoustics, electromagnetics, and elasticity. However, it remains a challenge to design robust and efficient numerical algorithms for the Helmholtz equation, especially with the large wave number or highly oscillatory solution.In recent years, many numerical methods have been developed to solve and analyze the Helmholtz equation, for instance, finite difference method [2–5], conforming finite element method [6, 7], boundary element method [8], weak Galerkin finite element method [9–12], spectral method [13], and adaptive finite element method [14] . It is also well known that the discontinuous Galerkin methods are flexible and highly parallelizable, and hence discontinuous Galerkin methods are widely used to solve the Helmholtz equation numerically, such as interior penalty discontinuous Galerkin method [15], hybridizable discontinuous Galerkin method [16, 17], local discontinuous Galerkin method [18], and the references therein. However, the local discontinuous Galerkin method [19–21] is known to be more physical and flexible on designing discontinuous Galerkin schemes. Two local discontinuous Galerkin methods are studied in [18], where the P1−1−P1−1 finite element pair was used to approximate the flux variable and the scattered field. They obtain the suboptimal convergence for the flux variable. In the numerical fluxes, the authors choose the penalty parameter β>0 in theoretical analysis and numerical experiments.In this paper, we set the penalty parameterβ=0 in numerical fluxes and use the discontinuous Pp+1−1−Pp−1p≥0 finite element pair to approximate the flux variable and the scattered field. In addition, we obtain the optimal convergence for the scattered field and the flux variable in the L2 norm and the flux variable in the Hdiv-like norm numerically.The paper is organized as follows. In Section2, we give the definition of the numerical fluxes and introduce the mixed discontinuous Galerkin method. A series of numerical experiments are presented to validate the effectiveness of the mixed discontinuous Galerkin method in Section 3. The paper concludes with conclusions in Section 4.


## 2. Mixed Local Discontinuous Galerkin Method

### 2.1. Meshes and Notations

In order to establish the mixed discontinuous Galerkin method, we first introduce some notations. The standard space, norm, and inner product notations are adopted in this paper. In particular, for a bounded domainD, ⋅,⋅D and ·,·∂D denote the L2- inner product on complex-valued spaces L2D and L2∂D, respectively.LetTh be a shape regular triangulation of the domain Ω with mesh size h, hK be the diameter of any K∈Th and h=maxK∈ThhK, ℰh be the set of all edges in Th, ℰhI be the collection of all interior edges, and ℰhB=ℰh/ℰhI be the set of boundary edges. For an interior edge e shared by two elements K1 and K2, n1 and n2 are the unit outward normal vectors of them, respectively. For any scalar function v and vector-valued function τ, let vi=v∂Kiandτi=τ∂Kii=1,2, and the average · and jump · are defined as follows:(3)v=12v1+v2,τ=12τ1+τ2,e∈ℰhI,v=v1n1+v2n2,τ=τ1·n1+τ2·n2,e∈ℰhI.On a boundary edgee∈ℰhB, we set(4)v=v,τ=τ,v=vn,τ=τ⋅n,where n is the outward unit normal vector on ∂Ω.


### 2.2. Mixed Discontinuous Galerkin Scheme

In this section, we introduce the mixed discontinuous Galerkin method for the nonhomogeneous Helmholtz equations (1) and (2). Firstly, we introduce an auxiliary variable σ=−Δu and rewrite equation (1) into(5)σ+Δu=0,inΩ,(6)divσ−k2u=f,inΩ.Next, we introduce the mixed discontinuous Galerkin method to find the solution pairu,σ of (5) and (6) numerically. Multiplying (5) and (6) by test functions τ¯ and v¯, respectively, and integrating both equations over an element K∈Th yield to(7)∫Kσ⋅τ¯dx−∫Kudivτ¯dx+∫∂KunK⋅τ¯ds=0,(8)−∫Kσ⋅∇v¯dx−k2∫Kuv¯dx+∫∂Kσ⋅nKv¯ds=∫Kfv¯dx,where nK denotes the unit outward normal vector to ∂K and the overbar denotes complex conjugation. The mixed discontinuous Galerkin spaces Vh and Qh are defined as(9)Vh=v∈L2Ω;vK∈PrK,∀K∈Th,Qh=τ∈L2Ωd;τK∈PlKd,∀K∈Th,where PrKr≥1 stands for the set of all polynomials of degree less than or equal to r on K. Based on the weak formulations (7) and (8), the solution uh,σh of the mixed discontinuous Galerkin method can be defined by(10)∫Kσh·τ¯hdx−∫Kuhdivτ¯hdx+∫∂Ku^KnK·τ¯hds=0−∫Kσh·∇v¯hdx−k2∫Kuhv¯hdx+∫∂Kσ^K·nKv¯hds=∫Kfv¯hdx,for any vh,τh∈Vh×Qh, where u^K and σ^K are the numerical fluxes on the boundary ∂K. The numerical fluxes have to be suitably defined in order to ensure the stability and accuracy of the mixed discontinuous Galerkin method. Summing the above equations over all elements, K∈Th, and using the following integration by parts identity,(11)∑K∈Th∫∂Kvτ¯·nds=∑e∈ℰh∫evτ¯ds+∑e∈ℰhI∫evτ¯ds,we get(12)σh,τhΩ−uh,divhτhΩ+∑e∈ℰhu^K,τhe+∑e∈ℰhIu^K,τhe=0,(13)divhσh,vhΩ−k2uh,vhΩ+∑e∈ℰhIσ^K−σhvhe+∑e∈ℰhσ^K−σh,vhe=f,vhΩ.Here, for any piecewise smooth scalar functionv and vector-valued function τ, let divhτK=divτK on any element K∈Th. In this paper, we choose the numerical fluxes in (12) and (13) as follows:(14)σ^h=σh,u^h=uh+iηehe−1σh,one∈ℰhI,(15)σ^h=σh,u^h=ikσh⋅n−ikg,one∈ℰhB,where ηe is positive constant. Substituting the numerical fluxes (14) and (15) into (12) and (13), we get the mixed discontinuous Galerkin formulation of (1) and (2) as follows: finding uh,σh∈Vh×Qh such that(16)σh,τhΩ−uh,divhτhΩ+∑e∈ℰhIuh,τhe+∑e∈ℰhIiηehe−1σh,τhe+ikσh·n,τh·nΓ=ikg,τh·nΓdivhσh,vhΩ−k2uh,vhΩ−∑e∈ℰhIσh,vhe=f,vhΩ,for any vh,τh∈Vh×Qh.


## 2.1. Meshes and Notations

In order to establish the mixed discontinuous Galerkin method, we first introduce some notations. The standard space, norm, and inner product notations are adopted in this paper. In particular, for a bounded domainD, ⋅,⋅D and ·,·∂D denote the L2- inner product on complex-valued spaces L2D and L2∂D, respectively.LetTh be a shape regular triangulation of the domain Ω with mesh size h, hK be the diameter of any K∈Th and h=maxK∈ThhK, ℰh be the set of all edges in Th, ℰhI be the collection of all interior edges, and ℰhB=ℰh/ℰhI be the set of boundary edges. For an interior edge e shared by two elements K1 and K2, n1 and n2 are the unit outward normal vectors of them, respectively. For any scalar function v and vector-valued function τ, let vi=v∂Kiandτi=τ∂Kii=1,2, and the average · and jump · are defined as follows:(3)v=12v1+v2,τ=12τ1+τ2,e∈ℰhI,v=v1n1+v2n2,τ=τ1·n1+τ2·n2,e∈ℰhI.On a boundary edgee∈ℰhB, we set(4)v=v,τ=τ,v=vn,τ=τ⋅n,where n is the outward unit normal vector on ∂Ω.


## 2.2. Mixed Discontinuous Galerkin Scheme

In this section, we introduce the mixed discontinuous Galerkin method for the nonhomogeneous Helmholtz equations (1) and (2). Firstly, we introduce an auxiliary variable σ=−Δu and rewrite equation (1) into(5)σ+Δu=0,inΩ,(6)divσ−k2u=f,inΩ.Next, we introduce the mixed discontinuous Galerkin method to find the solution pairu,σ of (5) and (6) numerically. Multiplying (5) and (6) by test functions τ¯ and v¯, respectively, and integrating both equations over an element K∈Th yield to(7)∫Kσ⋅τ¯dx−∫Kudivτ¯dx+∫∂KunK⋅τ¯ds=0,(8)−∫Kσ⋅∇v¯dx−k2∫Kuv¯dx+∫∂Kσ⋅nKv¯ds=∫Kfv¯dx,where nK denotes the unit outward normal vector to ∂K and the overbar denotes complex conjugation. The mixed discontinuous Galerkin spaces Vh and Qh are defined as(9)Vh=v∈L2Ω;vK∈PrK,∀K∈Th,Qh=τ∈L2Ωd;τK∈PlKd,∀K∈Th,where PrKr≥1 stands for the set of all polynomials of degree less than or equal to r on K. Based on the weak formulations (7) and (8), the solution uh,σh of the mixed discontinuous Galerkin method can be defined by(10)∫Kσh·τ¯hdx−∫Kuhdivτ¯hdx+∫∂Ku^KnK·τ¯hds=0−∫Kσh·∇v¯hdx−k2∫Kuhv¯hdx+∫∂Kσ^K·nKv¯hds=∫Kfv¯hdx,for any vh,τh∈Vh×Qh, where u^K and σ^K are the numerical fluxes on the boundary ∂K. The numerical fluxes have to be suitably defined in order to ensure the stability and accuracy of the mixed discontinuous Galerkin method. Summing the above equations over all elements, K∈Th, and using the following integration by parts identity,(11)∑K∈Th∫∂Kvτ¯·nds=∑e∈ℰh∫evτ¯ds+∑e∈ℰhI∫evτ¯ds,we get(12)σh,τhΩ−uh,divhτhΩ+∑e∈ℰhu^K,τhe+∑e∈ℰhIu^K,τhe=0,(13)divhσh,vhΩ−k2uh,vhΩ+∑e∈ℰhIσ^K−σhvhe+∑e∈ℰhσ^K−σh,vhe=f,vhΩ.Here, for any piecewise smooth scalar functionv and vector-valued function τ, let divhτK=divτK on any element K∈Th. In this paper, we choose the numerical fluxes in (12) and (13) as follows:(14)σ^h=σh,u^h=uh+iηehe−1σh,one∈ℰhI,(15)σ^h=σh,u^h=ikσh⋅n−ikg,one∈ℰhB,where ηe is positive constant. Substituting the numerical fluxes (14) and (15) into (12) and (13), we get the mixed discontinuous Galerkin formulation of (1) and (2) as follows: finding uh,σh∈Vh×Qh such that(16)σh,τhΩ−uh,divhτhΩ+∑e∈ℰhIuh,τhe+∑e∈ℰhIiηehe−1σh,τhe+ikσh·n,τh·nΓ=ikg,τh·nΓdivhσh,vhΩ−k2uh,vhΩ−∑e∈ℰhIσh,vhe=f,vhΩ,for any vh,τh∈Vh×Qh.


## 3. Numerical Experiments

In this section, we present some numerical examples of the mixed discontinuous Galerkin method. In the following experiments, we set penalty parameterηe=10.


### 3.1. Example 1: Helmholtz Equation in Convex Domain

We first consider a Helmholtz problem defined on the domainΩ=−0.5,0.5×−0.5,0.5, and the exact solution is given by(17)u=coskrk−cosk+isinkkJ0k+iJ1kJ0kr,f=sinkrr,where r=x2+y2 and Jnz are the Bessel functions of the first kind and order n. The function g is chosen by the exact solution. We test uniformly triangular meshes as shown in Figure 1(a). The numerical errors and convergence rate for k=5 and k=10 with the discontinuous Galerkin pair Pp+1−1−Pp−1p=0,1 are given in Tables 1–4. We can observe that the convergence rate is optimal for the flux variable in both Hdiv-like norm and L2 norm and the scattered field in L2 norm.Figure 1
The meshes for example 3.1 (a) and example 3.2 (b).
(a)
(b)Table 1
Numerical errors and convergence rate withk=5 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
1.0505e – 01
—
4.0575e − 01
—
2.6311e + 00
—
1/4
5.2449e – 02
1.0020
1.3732e − 01
1.5631
1.3206e + 00
0.9945
1/8
2.7299e – 02
0.9421
3.8612e − 02
1.8304
6.8713e − 01
0.9425
1/16
1.3851e − 02
0.9788
9.9766e − 03
1.9524
3.4859e − 01
0.9790
1/32
6.9536e − 03
0.9942
2.5156e − 03
1.9877
1.7499e − 01
0.9942
1/64
3.4804e − 03
0.9985
6.3026e − 04
1.9969
8.7586e − 02
0.9985Table 2
Numerical errors and convergence rates withk=10 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
5.4580e − 02
—
9.9906e − 01
—
8.9505e + 00
—
1/4
2.1681e − 02
1.3319
5.0097e − 01
0.9958
3.2653e + 00
1.4547
1/8
1.8045e − 02
0.2648
2.1682e − 01
1.2082
1.5504e + 00
1.0746
1/16
1.2249e − 02
0.5590
6.3883e − 02
1.7630
1.1230e + 00
0.4653
1/32
6.6122e − 03
0.8894
1.6672e − 02
1.9381
6.1967e − 01
0.8578
1/64
3.3701e − 03
0.9723
4.2134e − 03
1.9843
3.1747e − 01
0.9649Table 3
Numerical errors and convergence rates withk=5 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
2.1579e − 02
—
4.3274e − 02
—
5.7183e − 01
—
1/4
6.6959e − 03
1.6883
7.0819e − 03
2.6113
1.8744e − 01
1.6092
1/8
1.7359e − 03
1.9476
8.8260e − 04
3.0043
4.8484e − 02
1.9508
1/16
4.3740e – 04
1.9886
1.1042e − 04
2.9987
1.2215e − 02
1.9889
1/32
1.0956e − 04
1.9973
1.3799e − 05
3.0003
3.0595e − 03
1.9973
1/64
2.7402e − 05
1.9993
1.5553e − 06
3.1493
7.6522e − 04
1.9993Table 4
Numerical errors and convergence rates withk=10 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
4.3837e − 02
—
5.3549e − 01
—
4.7389e + 00
—
1/4
1.4489e − 02
1.5972
7.2364e − 02
2.8875
1.2945e + 00
1.8722
1/8
4.1745e − 03
1.7953
6.4364e − 03
3.4909
3.9645e − 01
1.7072
1/16
1.0708e − 03
1.9629
8.4948e − 04
2.9216
1.0206e − 01
1.9578
1/32
2.6887e − 04
1.9937
1.1549e − 04
2.8788
2.5641e − 02
1.9928
1/64
6.7280e − 05
1.9986
1.3821e − 05
3.0629
6.4172e − 03
1.9984


### 3.2. Example 2: Helmholtz Equation in Nonconvex Domain

In this numerical test, we solve the Helmholtz problem defined on a nonconvex domain as shown in Figure1(b). We take f=0 and choose the exact solution as follows:(18)u=Jξkx+12+y2cosξarctanyx+1.The boundary conditiong is given by the exact solution u for ξ=1,ξ=3/2, and ξ=2/3. We can check that u is smooth for ξ∈ℕ, while its derivative has a singularity at −1,0 for ξ∉ℕ. For ξ=1,3/2,and2/3 and wave number k=4, Figures 2 and 3 show the errors of the mixed discontinuous Galerkin method with the finite element pair P1−1−P0−1 and P2−1−P1−1, respectively. It shows that convergent rates are optimal in three errors for cases ξ=1. We lose the optimal convergent in the flux variable σ−σh0,Ω, while the scalar variable u−uh0,Ω and the flux variable divσ−σh0,Ω still converge at the optimal order, for the case ξ=2/3 and ξ=3/2.Figure 2
Convergence rates withk=4andξ=1,2/3,and3/2 for the P1−1−P0−1 element. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)Figure 3
Convergence rates withk=4andξ=1,2/3,and3/2 for the P2−1−P1−1 element. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)


### 3.3. Example 3: Helmholtz Equation with Large Wave Number

This test is devoted to studying the numerical test for the same setting as example 1 with a large wave number. We solve the Helmholtz problem with different mesh sizes for four wave numbersk=5,10,50,and100. The errors in u−uh0,Ω, σ−σh0,Ω, and divσ−σh0,Ω are shown in Figure 4 with p=0 and Figure 5 with p=1. It indicates that the mixed discontinuous Galerkin method is convergent for the cases k=5 and k=10, and the errors begin to reduce as h becomes to be quite small for large wave numbers k=50 and 100. The surface plots of the mixed discontinuous Galerkin solution and the exact solution are shown in Figure 6. It shows that the mixed discontinuous Galerkin solution has the correct shape.Figure 4
The errors for four wave numbersk=5,10,50,and100andP1−1−P0−1. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)Figure 5
The errors for four wave numbersk=5,10,50,and100andP2−1−P1−1. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)Figure 6
The surface plot fork=100,h=1/64,andP1−1−P0−1. (a) Numerical solution. (b) Exact solution.
(a)
(b)


### 3.4. Example 4: Helmholtz Equation with Perturbation

We consider the Helmholtz equation with perturbation [14] in the following form:(19)−Δu−k21+qxu=f,inΩ∈R2,∂u∂nΩ+iku=g,onΓ=∂Ω,where f=k2qxuix,g=−∂ui/∂nΩ−ikui, uix=expikx is the incident field and the perturbation function qx has compact support in the domain Ω. Note that the solution is propagating for qx>−1, while it is evanescent for qx<−1. The numerical fluxes for the Helmholtz equation with perturbation are the same as (14) and (15). In this test, let Ω=−1.5,1.5×−1.5,1.5. Due to the area of compact support of perturbation qx which is little relative to the entire area Ω, with the aim of better observing the experimental phenomena of the scattered field, we limit the interval to −1,1×−1,1 as shown in Figures 7 and 8.(i)
Take the perturbation functionqx=1.5exp−160x2 with wave number k=40. The graphics of the real part of the scattered field and perturbation function is shown in Figure 7. We can find that the wave propagates through the origin of the positive perturbation function.(ii)
In Figure8, we present the geometry of the real part of the scattered field and perturbation function with the perturbation function qx=−1.5exp−160x2 and wave number k=40. The results show that due to the scattering of the plane wave impinged on the negative perturbation, the wave gradually disappears near the origin.Figure 7
Scattered field and perturbation function withk=40,h=1/256,andP1−1−P0−1. (a) Real part of the scattered field. (b) Perturbation function.
(a)
(b)Figure 8
Scattered field and perturbation function withk=40,h=1/256,andP1−1−P0−1. (a) Real part of he scattered field. (b) Perturbation function.
(a)
(b)


### 3.5. Example 5: Rectangular Waveguide

In this example, we consider the domainΩ is Ω=0,1×0,1, and the 2D rectangular waveguide model problem satisfies the Helmholtz equation (1) with the following boundary conditions:(20)∂u∂n=0,onΓ1=0,1×0andΓ3=0,1×1,(21)∂u∂n=g2,onΓ4=0×0,1,(22)∂u∂n=−iku,onΓ4=1×0,1.From [3, 22], we set in (1) and g2=e−1/2y−1/22 in (21). Similarly, we define the numerical fluxes for the rectangular waveguide (1) and (20) and (22) as the following:(23)σ^h=σh,u^h=uh+iηehe−1σh,one∈εhI,σ^h=σh,u^h=−ikσh⋅n,one∈Γ2,σ^h=−g2n,u^h=uh,one∈Γ4,σ^h=0,u^h=uh,one∈Γ1∪Γ3.The real part of the waveguide solutions fork=16π,h=1/64 and k=64π,h=1/128 is shown in Figure 9, which is consistent with the results in [3, 22].Figure 9
Real part of the waveguide solution withP1−1−P0−1. (a) k=16π,h=1/64. (b) k=64π,h=1/128.
(a)
(b)


## 3.1. Example 1: Helmholtz Equation in Convex Domain

We first consider a Helmholtz problem defined on the domainΩ=−0.5,0.5×−0.5,0.5, and the exact solution is given by(17)u=coskrk−cosk+isinkkJ0k+iJ1kJ0kr,f=sinkrr,where r=x2+y2 and Jnz are the Bessel functions of the first kind and order n. The function g is chosen by the exact solution. We test uniformly triangular meshes as shown in Figure 1(a). The numerical errors and convergence rate for k=5 and k=10 with the discontinuous Galerkin pair Pp+1−1−Pp−1p=0,1 are given in Tables 1–4. We can observe that the convergence rate is optimal for the flux variable in both Hdiv-like norm and L2 norm and the scattered field in L2 norm.Figure 1
The meshes for example 3.1 (a) and example 3.2 (b).
(a)
(b)Table 1
Numerical errors and convergence rate withk=5 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
1.0505e – 01
—
4.0575e − 01
—
2.6311e + 00
—
1/4
5.2449e – 02
1.0020
1.3732e − 01
1.5631
1.3206e + 00
0.9945
1/8
2.7299e – 02
0.9421
3.8612e − 02
1.8304
6.8713e − 01
0.9425
1/16
1.3851e − 02
0.9788
9.9766e − 03
1.9524
3.4859e − 01
0.9790
1/32
6.9536e − 03
0.9942
2.5156e − 03
1.9877
1.7499e − 01
0.9942
1/64
3.4804e − 03
0.9985
6.3026e − 04
1.9969
8.7586e − 02
0.9985Table 2
Numerical errors and convergence rates withk=10 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
5.4580e − 02
—
9.9906e − 01
—
8.9505e + 00
—
1/4
2.1681e − 02
1.3319
5.0097e − 01
0.9958
3.2653e + 00
1.4547
1/8
1.8045e − 02
0.2648
2.1682e − 01
1.2082
1.5504e + 00
1.0746
1/16
1.2249e − 02
0.5590
6.3883e − 02
1.7630
1.1230e + 00
0.4653
1/32
6.6122e − 03
0.8894
1.6672e − 02
1.9381
6.1967e − 01
0.8578
1/64
3.3701e − 03
0.9723
4.2134e − 03
1.9843
3.1747e − 01
0.9649Table 3
Numerical errors and convergence rates withk=5 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
2.1579e − 02
—
4.3274e − 02
—
5.7183e − 01
—
1/4
6.6959e − 03
1.6883
7.0819e − 03
2.6113
1.8744e − 01
1.6092
1/8
1.7359e − 03
1.9476
8.8260e − 04
3.0043
4.8484e − 02
1.9508
1/16
4.3740e – 04
1.9886
1.1042e − 04
2.9987
1.2215e − 02
1.9889
1/32
1.0956e − 04
1.9973
1.3799e − 05
3.0003
3.0595e − 03
1.9973
1/64
2.7402e − 05
1.9993
1.5553e − 06
3.1493
7.6522e − 04
1.9993Table 4
Numerical errors and convergence rates withk=10 for the P1−1−P0−1 element.
h
u
−
u
h
0
Rate
σ
−
σ
h
0
Rate
div
σ
−
σ
h
0
Rate
1/2
4.3837e − 02
—
5.3549e − 01
—
4.7389e + 00
—
1/4
1.4489e − 02
1.5972
7.2364e − 02
2.8875
1.2945e + 00
1.8722
1/8
4.1745e − 03
1.7953
6.4364e − 03
3.4909
3.9645e − 01
1.7072
1/16
1.0708e − 03
1.9629
8.4948e − 04
2.9216
1.0206e − 01
1.9578
1/32
2.6887e − 04
1.9937
1.1549e − 04
2.8788
2.5641e − 02
1.9928
1/64
6.7280e − 05
1.9986
1.3821e − 05
3.0629
6.4172e − 03
1.9984


## 3.2. Example 2: Helmholtz Equation in Nonconvex Domain

In this numerical test, we solve the Helmholtz problem defined on a nonconvex domain as shown in Figure1(b). We take f=0 and choose the exact solution as follows:(18)u=Jξkx+12+y2cosξarctanyx+1.The boundary conditiong is given by the exact solution u for ξ=1,ξ=3/2, and ξ=2/3. We can check that u is smooth for ξ∈ℕ, while its derivative has a singularity at −1,0 for ξ∉ℕ. For ξ=1,3/2,and2/3 and wave number k=4, Figures 2 and 3 show the errors of the mixed discontinuous Galerkin method with the finite element pair P1−1−P0−1 and P2−1−P1−1, respectively. It shows that convergent rates are optimal in three errors for cases ξ=1. We lose the optimal convergent in the flux variable σ−σh0,Ω, while the scalar variable u−uh0,Ω and the flux variable divσ−σh0,Ω still converge at the optimal order, for the case ξ=2/3 and ξ=3/2.Figure 2
Convergence rates withk=4andξ=1,2/3,and3/2 for the P1−1−P0−1 element. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)Figure 3
Convergence rates withk=4andξ=1,2/3,and3/2 for the P2−1−P1−1 element. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)


## 3.3. Example 3: Helmholtz Equation with Large Wave Number

This test is devoted to studying the numerical test for the same setting as example 1 with a large wave number. We solve the Helmholtz problem with different mesh sizes for four wave numbersk=5,10,50,and100. The errors in u−uh0,Ω, σ−σh0,Ω, and divσ−σh0,Ω are shown in Figure 4 with p=0 and Figure 5 with p=1. It indicates that the mixed discontinuous Galerkin method is convergent for the cases k=5 and k=10, and the errors begin to reduce as h becomes to be quite small for large wave numbers k=50 and 100. The surface plots of the mixed discontinuous Galerkin solution and the exact solution are shown in Figure 6. It shows that the mixed discontinuous Galerkin solution has the correct shape.Figure 4
The errors for four wave numbersk=5,10,50,and100andP1−1−P0−1. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)Figure 5
The errors for four wave numbersk=5,10,50,and100andP2−1−P1−1. (a) u−uh0. (b) σ−σh0. (c) divσ−σh0.
(a)
(b)
(c)Figure 6
The surface plot fork=100,h=1/64,andP1−1−P0−1. (a) Numerical solution. (b) Exact solution.
(a)
(b)


## 3.4. Example 4: Helmholtz Equation with Perturbation

We consider the Helmholtz equation with perturbation [14] in the following form:(19)−Δu−k21+qxu=f,inΩ∈R2,∂u∂nΩ+iku=g,onΓ=∂Ω,where f=k2qxuix,g=−∂ui/∂nΩ−ikui, uix=expikx is the incident field and the perturbation function qx has compact support in the domain Ω. Note that the solution is propagating for qx>−1, while it is evanescent for qx<−1. The numerical fluxes for the Helmholtz equation with perturbation are the same as (14) and (15). In this test, let Ω=−1.5,1.5×−1.5,1.5. Due to the area of compact support of perturbation qx which is little relative to the entire area Ω, with the aim of better observing the experimental phenomena of the scattered field, we limit the interval to −1,1×−1,1 as shown in Figures 7 and 8.(i)
Take the perturbation functionqx=1.5exp−160x2 with wave number k=40. The graphics of the real part of the scattered field and perturbation function is shown in Figure 7. We can find that the wave propagates through the origin of the positive perturbation function.(ii)
In Figure8, we present the geometry of the real part of the scattered field and perturbation function with the perturbation function qx=−1.5exp−160x2 and wave number k=40. The results show that due to the scattering of the plane wave impinged on the negative perturbation, the wave gradually disappears near the origin.Figure 7
Scattered field and perturbation function withk=40,h=1/256,andP1−1−P0−1. (a) Real part of the scattered field. (b) Perturbation function.
(a)
(b)Figure 8
Scattered field and perturbation function withk=40,h=1/256,andP1−1−P0−1. (a) Real part of he scattered field. (b) Perturbation function.
(a)
(b)


## 3.5. Example 5: Rectangular Waveguide

In this example, we consider the domainΩ is Ω=0,1×0,1, and the 2D rectangular waveguide model problem satisfies the Helmholtz equation (1) with the following boundary conditions:(20)∂u∂n=0,onΓ1=0,1×0andΓ3=0,1×1,(21)∂u∂n=g2,onΓ4=0×0,1,(22)∂u∂n=−iku,onΓ4=1×0,1.From [3, 22], we set in (1) and g2=e−1/2y−1/22 in (21). Similarly, we define the numerical fluxes for the rectangular waveguide (1) and (20) and (22) as the following:(23)σ^h=σh,u^h=uh+iηehe−1σh,one∈εhI,σ^h=σh,u^h=−ikσh⋅n,one∈Γ2,σ^h=−g2n,u^h=uh,one∈Γ4,σ^h=0,u^h=uh,one∈Γ1∪Γ3.The real part of the waveguide solutions fork=16π,h=1/64 and k=64π,h=1/128 is shown in Figure 9, which is consistent with the results in [3, 22].Figure 9
Real part of the waveguide solution withP1−1−P0−1. (a) k=16π,h=1/64. (b) k=64π,h=1/128.
(a)
(b)


## 4. Conclusion

The paper developed the mixed discontinuous Galerkin method for the Helmholtz equation. We define the numerical fluxes and introduce the mixed Galerkin method. Compared with [18], we set penalty parameter β=0 in numerical fluxes and use the discontinuous PP+1−1−Pp−1p≥0 finite element pair to approximate the flux variable and the scattered field. We present several numerical experiments to validate the effectiveness of the mixed discontinuous Galerkin method. The optimal convergence for the scattered field and the flux variable in the L2 norm and the flux variable in the Hdiv-like norm is obtained numerically. Furthermore, we test the Helmholtz equation with a large wave number and the Helmholtz equation with perturbation function and the rectangular waveguide, which indicate the well performance of the mixed discontinuous Galerkin method.



---

*Source: 9582583-2020-05-04.xml*