Latest revision as of 21:13, 27 December 2025

1 Introduction

The modified Helmholtz equation appears in a common physical model, and its numerical solution remains an active line of research. See the discussion in Yaman & zdemir [1] and the recent Yaman et al [2]. In the former, an integral equation method is proposed for the numerical solution. In this work, a Galerkin approach is preferred for the numerical solution of perturbations of the modified Helmholtz equation.

The classical Galerkin approach is the continuous Galerkin, Finite Element Method (FEM). A recent and attractive alternative is the Discontinuous Galerkin (DG) method. Unlike FEM, it is locally conservative and ${\textstyle H^{p}}$-adaptative. On the downside, the application of the Discontinuous Galerkin Method to elliptic problems usually leads to underdetermined linear systems, and penalization or suitable constraints are required. Moreover, Discontinuous Galerkin methods are more expensive because of the need for numerical fluxes at the edges of the mesh elements, thus yielding more coupled unknowns than FEM. See Riviere [3] and Di Pietro & Ern [4] for a thorough exposition.

These shortcomings of the DG method can be addressed in part by its descendant, the Hybrid Discontinuous Galerkin Method (HDG), Bui-Tahn [5]. The numerical flux is hybridized, introducing additional unknowns on the edges of the mesh, which reduces the coupling between elements. The global solution is then obtained by solving small and independent linear systems on each element.

For the elliptic problem under study, we use a simple, physically motivated hybrid numerical flux that yields a well-posed linear system. This is guaranteed by the dominant reaction term in the modified Helmholtz equation.

The outline is as follows. In Section 2, we follow the HDG methodology and introduce a first-order system associated with the perturbed modified Helmholtz equations. Then, we delve into the discretization of the first-order system and show that the resulting finite-dimensional problem is well posed. Numerical tests are presented in Section 3. First, we illustrate its performance as a Poisson solver in a preliminary comparison with a leading method in the literature. Accurate approximations are also shown for some benchmark Helmholtz problems in Coastal Ocean Modeling. In all examples, we highlight approximations on coarse meshes. We close our exposition with conclusions and future work.

2 Numerical solution of the modified Helmholtz equation

In what follows, we build on the theoretical and numerical aspects of the finite element method (FEM) as presented, for instance, in the text [6].

2.1 A first order system for perturbations of the modified Helmholtz equation

For a bounded Lipschitz domain ${\textstyle \Omega }$ in ${\textstyle \mathbb {R} ^{d}}$, let us consider the diffusion-advection-reaction equation for the unknown function ${\textstyle u}$,

The diffusion term is uniformly elliptic, that is ${\textstyle \mathbf {D} (x)}$ is symmetric and there exist constants ${\textstyle 0<\lambda {<\Lambda }}$ such that for almost all ${\textstyle x\in \Omega }$

Also, assume that ${\textstyle \mathbf {v} \in \left(L^{\infty }(\Omega ))\right)^{d}}$, and ${\textstyle c,\,b\in L^{\infty }(\Omega )}$.

For ${\textstyle \mathbf {v} =0}$ and ${\textstyle c>0}$, equation (1) is known as the modified Helmholtz equation. To construct a Discontinuous Galerkin approximation, a first step is to construct a hyperbolic-like first-order system. A small advection term is allowed; the smallness condition is made precise below.

Define the new variable (flux) ${\textstyle \mathbf {z} =-\mathbf {D} \nabla u}$. Since the matrix ${\textstyle \mathbf {D} }$ is invertible, ${\textstyle \nabla u=-\mathbf {D} ^{-1}\mathbf {z} }$, we obtain the following system

Let ${\textstyle m=d+1}$, and ${\textstyle \mathbf {e} ^{i}}$ be the ${\textstyle i-}$th canonical vector. Define

The system becomes

2.2 A Hybrid Discontinuous Galerkin Method

Here we introduce the essentials of HDG, for full details see Bui-Thanh [5].

For any matrix ${\textstyle \mathbf {M} }$, we denote by ${\textstyle \mathbf {M} _{i},\,\mathbf {M} ^{j}}$ its row ${\textstyle i}$ and column ${\textstyle j}$ respectively.

Let us define

and

We apply operators component-wise. For instance, for the divergence operator, we have,

We can write system (3) in the form

Assume we have a valid element partition of the domain ${\textstyle \Omega }$. In DG, functions are approximated locally in each element by polynomials, then coupled with others in adjacent elements by means of a numerical flux. The basic process is as follows.

Let ${\textstyle \tau }$ be an element in the partition. Compute the ${\textstyle (L^{2})^{m}}$ inner product on ${\textstyle \tau }$ of each side of (4) with a test function ${\textstyle \mathbf {v} }$, to obtain

Denoting the inner product in th eboundary of ${\textstyle \tau }$ by ${\textstyle \langle \cdot ,\cdot \rangle _{\partial \tau }}$ we obtain after integrating by parts,

Continuity is not enforced at the boundary of adjacent elements. Therefore, the boundary term ${\textstyle F(\mathbf {u} )}$ is replaced with a boundary numerical flux ${\textstyle F^{*}(\mathbf {u} ^{*})}$ where ${\textstyle \mathbf {u} ^{*}\equiv \mathbf {u} ^{*}(\mathbf {u} ^{-},\mathbf {u} ^{+})}$ solves a Riemann problem with Cauchy data ${\textstyle \mathbf {u} ^{-},\,\mathbf {u} ^{+}}$. As is standard, the – superscript denotes limits from the interior of ${\textstyle e}$, and the ${\textstyle +}$ superscript, limits from the exterior. In this context, element ${\textstyle \tau }$ is denoted by ${\textstyle \tau ^{-}}$ and the outer normal ${\textstyle \mathbf {n} }$ by ${\textstyle \mathbf {n} ^{-}}$.

To hybridize the flux, and break the coupling, ${\textstyle \mathbf {u} ^{*}}$ is regarded as an extra unknown to be solved on the skeleton (the set of element edges ${\textstyle {\mathcal {F}}}$) of the mesh. Renaming ${\textstyle \mathbf {u} ^{*}}$ as ${\textstyle {\hat {\mathbf {u} }}}$ and ${\textstyle F^{*}}$ as ${\textstyle {\hat {F}}}$, the problem is to propose a suitable hybrid numerical flux ${\textstyle {\hat {\mathbf {F} }}(\mathbf {u} )}$.

Let us apply this scheme to equation (1) under the following condition

This condition implies a dominant reaction term, which yields a tamed advection. Consequently, the diffusive flux is assumed to be continuous across any interface. Instead of the Riemann problem solution approach, we set the hybrid flux

Notice the dependence only in the unknown ${\textstyle {\hat {u}}}$.

For each element ${\textstyle \tau }$, the DG local unknown ${\textstyle \mathbf {u} }$ and the extra trace unknown ${\textstyle {\hat {u}}}$ need to satisfy

This is complemented by a weak jump condition in the skeleton. Namely, for each edge ${\textstyle e}$

where ${\textstyle \mathbf {v} }$ in (9), ${\textstyle w}$ in (10) are polynomial test functions defined on elements ${\textstyle \tau }$ and edges ${\textstyle e}$ respectively.

Here, ${\textstyle \lbrack \lbrack \cdot \rbrack \rbrack }$ is the jump operator,

On each element, we are led to solving the equations (9) and (10). Let us show that this finite-dimensional problem is well-posed by showing that under null data, the solution is the trivial one.

Lemma. Assume condition (7) holds. If ${\textstyle f(x)=0}$ and ${\textstyle {\hat {u}}=0}$ on ${\textstyle {\mathcal {F}}}$, then ${\textstyle u=0,\mathbf {z} =\mathbf {0} }$ in ${\textstyle \Omega }$.

Proof.

Let us split the test function in the form ${\textstyle \mathbf {v} =(\mathbf {w} ,w)}$. The equation (9) using the hybrid flux can be written as

By integration by parts in (12)

Let ${\textstyle \mathbf {w} =\mathbf {z} ,w=u}$, and add equations (11) and (13) to obtain

.

Notice that in the boundary term, ${\textstyle u}$ is the inner limit ${\textstyle u^{-}}$.

We are led to

Hence, ${\textstyle u^{-}=0}$ on ${\textstyle \partial \tau }$, and ${\textstyle \mathbf {z} =\mathbf {0} }$ in ${\textstyle \tau }$. Equation (11) becomes

Integrating by parts

Since ${\textstyle \mathbf {w} }$ is arbitrary, ${\textstyle \nabla u=0}$, and consequently ${\textstyle u\equiv 0}$ in each element ${\textstyle u}$.

We conclude that ${\textstyle u=0,\mathbf {z} =\mathbf {0} }$ in ${\textstyle \Omega }$. ${\textstyle \blacksquare }$

Remark. Notice that the result is valid for Poisson problems, ${\textstyle c=0}$, ${\textstyle \mathbf {v} =\mathbf {0} }$.

3 Numerical results

Here, we illustrate some numerical results on a variety of problems, where the perturbed modified Helmholtz equation HDG solver serves as the computational engine.

The full code of the numerical implementation is available on GitHub:

https://github.com/Danalie/Hibrid-Discontinuos-Galerkin-DG2

We will test benchmark problems in rectangular geometries and regular meshes for simplicity.

We use the classical element functions on the reference segment ${\textstyle [-1,1]}$ for 1D applications. We denote by ${\textstyle HDG_{p}}$ a ${\textstyle p}$-order bilinear approximation. In 2D, ${\textstyle HDG_{2},HDG_{4}}$ refer to basis functions on the reference square ${\textstyle [-1,1]\times [-1,1]}$, as products of first-order and second-order one-dimensional 1D basis functions, respectively.

To report the accuracy of the approximation, we use the Root Mean Square Error (RMSE). Sometimes, we normalize by the corresponding norm of the exact solution to obtain the Normalized Root Mean Square Error (NRMSE).

3.1 Poisson problems

For ${\textstyle c=0}$ and ${\textstyle \mathbf {v} =0}$ in (1) we have the pure diffusion equation

To illustrate the HDG method with continuous diffusion numerical flux, as a Poisson solver, we consider ${\textstyle \mathbf {D} }$ as the identity matrix.

First a smooth Dirichlet Problem for Poisson's equation on the domain ${\textstyle \Omega =(0,1)\times (0,1)}$,

The exact solution is ${\textstyle p(x,y)=1+x+y+xy}$. The numerical solution using ${\textstyle HDG_{2}}$ on a ${\textstyle 10\times 10}$ grid yields ${\textstyle RMSE=0}$ to machine precision; no further refinement is necessary.

The DG method is ${\textstyle H^{p}}$-adaptative. Local high order is straightforward, and accurate approximations are possible on coarse meshes. For instance, we use ${\textstyle HDG_{4}}$ for the Dirichlet problem

On a ${\textstyle 10\times 10}$ grid the ${\textstyle RMSE}$ is again zero to machine precision.

A preliminary comparison with the method in [7]. The latter relies on an internal-penalty numerical flux requiring a penalty parameter chosen heuristically, see (34) therein. Accuracy is shown in an exhaustive list of test problems. The simplest of such is the Dirichlet problem

with analytic solution is ${\textstyle q_{e}=\sin(\pi x)\sin(\pi y)}$. Our solution is in Table 1.

By inspection of Figure 7 in [7], it is apparent that our results compare favorably. A full comparison is to be carried out.

3.2 Helmholtz equation

We consider two benchmark problems in Coastal Ocean Modeling (COM) that lead to Helmholtz equations. A comparison is made with the Finite Volume (FVCOM) as presented in Chen et al [8]. Therein, FVCOM is applied for the modeling of tidal simulation in a semienclosed basin with tidal forcing at the open boundary under near-resonance conditions.

Here we apply HDG to illustrate the accurate simulation of the troublesome near resonance case. We remark that for the Helmholtz equation, condition (7) is not valid. Nevertheless, the numerical results are highly satisfactory.

3.2.1 A Rectangular Channel

Consider a fluid layer of uniform density that propagates along a channel aligned with the ${\textstyle x-}$direction. More precisely, a semienclosed narrow channel with length ${\textstyle L}$ and variable depth ${\textstyle H(x)}$ and a closed boundary at ${\textstyle x=L_{1}}$ and an open boundary at ${\textstyle x=L}$.

Neglecting Coriolis force and advection of momentum, the governing equations modeling tidal waves propagation in the semienclosed channel (see Figure 1) are given as

Here, ${\textstyle H(x)}$ is the total water depth, ${\textstyle g}$ is acceleration due to gravity, ${\textstyle u}$ is speed in the ${\textstyle x-}$direction, and ${\textstyle \zeta }$ is sea-level elevation.

Figura 1: Configuration of the semienclosed channel.

Assuming harmonic solutions,

we obtain the one dimensional Helmholtz problem for ${\textstyle \zeta _{0}}$

We specify a periodic tidal forcing with amplitude ${\textstyle A}$ at the mouth of the channel,

and a no-flux boundary condition at the wall,

Let the channel depth decrease linearly toward the end of the channel, so that ${\textstyle H(x)}$ can be written as

It is readily seen that (14) is a Bessel's equation. With the given boundary conditions, the analytic solution is given by

where

${\textstyle J_{0},Y_{0}}$ are the Bessel's functions of degree zero and one respectively.

To compare with the HDG approximation, we solve the system,

Here we have defined ${\textstyle z=-H\zeta _{0}^{\prime }}$ in (14).

The following parameters are considered for a channel very close to resonance.

The HDG method is applied using nodal polynomials of degree 2 (${\textstyle {\mbox{HDG}}_{2}}$). It is compared with the FV method and the analytic solution. For consistency with the FV method, we consider the approximation at the middle point ${\textstyle x_{i}}$ of the element ${\textstyle I_{i}=[x_{i-{\frac {1}{2}}},x_{i+{\frac {1}{2}}}]}$. As illustrated in Table 4, a second-order approximation with HDG in a coarse resolution is of greater quality than FV.

The analytic solution describes a standing wave with a node point near the closed side of the channel. As expected, the reproduction of these features by HDG is accurate. See Figure 2

Figura 2: Solution for the rectangular channel near resonance case, with HDG nodal with polynomials of degree 2 and 80 elements.

Remark. As pointed out in Chen et al [8], regardless of the numerical method, a proper selection of horizontal resolution recover accurately this tidal resonance problem. We argue that the accuracy attained by HDG in coarse meshes yields a better choice.

3.2.2 A Sector Channel

Now consider a flat bottom channel in the form of a semicircular section, which in polar coordinates is defined from ${\textstyle 0}$ to ${\textstyle L}$ in the radial direction and from ${\textstyle -\alpha /2}$ to ${\textstyle \alpha /2}$ in the angular direction.

The semicircular line of radius ${\textstyle L}$ corresponds to an open border, while along the semicircular line of radius ${\textstyle L_{1}}$ and the two sides, they are closed, (Figure 3).

Figura 3: Semienclosed sector channel in the polar region ${\displaystyle L_{1}\leq r\leq L}$, ${\displaystyle -{\frac {\alpha }{2}}\leq \theta \leq {\frac {\alpha }{2}}}$. The sector is open at ${\displaystyle r=L}$ and closed elsewhere.

The following equations govern the nonrotating tidal oscillation,

${\textstyle H_{0}}$ is the constant water depth, ${\textstyle V_{r},V_{\theta }}$ are the radial and angular ${\textstyle r,\theta }$ velocity components, and ${\textstyle \eta }$ is the free surface water elevation.

Assuming a harmonic solution,

we can reduce the equations (15-17) to an elliptic equation

The physical boundary conditions are as follows

1. At the open mouth of the channel, a harmonic tidal forcing is assumed,

${\displaystyle \eta _{0}(r,\theta )={\overline {A}}\cos \left(m\pi {\frac {\theta +{\frac {\alpha }{2}}}{\alpha }}\right),\quad \lbrace L\rbrace \times \left(-{\frac {\alpha }{2}},{\frac {\alpha }{2}}\right),}$

2. On the solid walls, null flux is prescribed,

${\displaystyle -K\nabla \eta _{0}(r,\theta )=0,\quad \lbrace L_{1}\rbrace \times \left({\frac {-\alpha }{2}},{\frac {\alpha }{2}}\right)\bigcup (L_{1},L)\times \left\lbrace -{\frac {\alpha }{2}}\right\rbrace \bigcup (L_{1},L)\times \left\lbrace {\frac {\alpha }{2}}\right\rbrace {.}}$

The analytic solution of this boundary-value problem is:

where

${\textstyle J_{v},Y_{v}}$ are the ${\textstyle v}$ th-order Bessel function of the first and second types.

Let us show the HDG solution in the rectangular domain ${\textstyle (L_{1},L]\times (-\alpha /2,\alpha /2)}$.

Let ${\textstyle \nabla }$ denote the gradient with respect to ${\textstyle (r,\theta )}$ and let

Equation (18) becomes a perturbation of the two-dimensional Helmholtz equation,

We apply the fourth order (${\textstyle {\mbox{HDG}}_{4}}$) scheme developed above for a near-resonance case. The geometric parameter values are:

Figure (4) compares the analytic and numerical solutions. The relative mean square error is shown in Table 3.

Figura 4: Comparison between the analytic solution to the sector problem with the ${\displaystyle {\mbox{HDG}}_{4}}$ solution using ${\displaystyle 30\times 50}$ elements. a) Solution for some ${\displaystyle r}$ values, b) Solution for some ${\displaystyle \theta }$ values.

Remark. Noteworthy, the ${\textstyle 10\times 10}$ ${\textstyle {\mbox{HDG}}_{4}}$ solution is of grater quality that the ${\textstyle 40\times 40}$ FV solution. In regards to execution time, the former is attained in half a second on a personal laptop. A solution of the same quality would require in the order of minutes with the Finite Volume Method in a much finer mesh.

4 Conclussions

We have introduced a first-order system formulation for a regularly perturbed Modified Helmholtz equation. A discretization is derived by means of the Hybrid Discontinuous Galerkin method with an ad-hoc numerical flux. The well-posedness of the finite-dimensional HDG discrete linear system follows from a dominant reaction condition.

The preliminary comparison with a leading Poisson solver shows promise as a competitive alternative. Research on this is ongoing.

The proposed HDG discrete method is also applied to Helmholtz's problems arising from coastal ccean modeling. It outperforms the Finite Volume method commonly used in this setting.

The application to irregular regions with more general meshes is more elaborate but straightforward. Also, a parallel computing implementation is desired. Both tasks are left for future work.

BIBLIOGRAPHY

[1] Yaman, Olha Ivanyshyn and zdemir, Gazi. (2021) "Numerical solution of a generalized boundary value problem for the modified Helmholtz equation in two dimensions", Volume 190. Elsevier. Mathematics and Computers in Simulation 181–191

[2] Ivanyshyn Yaman, Olha and zdemir, Gazi. (2025) "An interior inverse generalized impedance problem for the modified Helmholtz equation in two dimensions", Volume 105. Wiley Online Library. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 1 e202300711

[3] Riviere, Béatrice. (2008) "Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation". SIAM

[4] Di Pietro, Daniele Antonio and Ern, Alexandre. (2011) "Mathematical aspects of discontinuous Galerkin methods", Volume 69. Springer Science & Business Media

[5] Bui-Thanh, Tan. (2015) "From Godunov to a unified hybridized discontinuous Galerkin framework for partial differential equations", Volume 295. Elsevier. Journal of Computational Physics 114–146

[6] Gockenbach, Mark S. (2006) "Understanding and implementing the finite element method". SIAM

[7] Fischer, Nils L and Pfeiffer, Harald P. (2022) "Unified discontinuous Galerkin scheme for a large class of elliptic equations", Volume 105. APS. Physical Review D 2 024034

[8] Chen, Changsheng and Huang, Haosheng and Beardsley, Robert C and Liu, Hedong and Xu, Qichun and Cowles, Geoffrey. (2007) "A finite volume numerical approach for coastal ocean circulation studies: Comparisons with finite difference models", Volume 112. Wiley Online Library. Journal of Geophysical Research: Oceans C3