An Efficient Solver to Helmholtz Equations by Recontruction Discontinuous Approximation

In this paper, an efficient solver for the Helmholtz equation using a noval approximation space is developed. The ingradients of the method include the approximation space recently proposed, a discontinuous Galerkin scheme extensively used, and a linear system solver with a natural preconditioner. Comparing to traditional discontinuous Galerkin methods, we refer to the new method as being more efficient in the following sense. The numerical performance of the new method shows that: 1) much less error can be reached using the same degrees of freedom; 2) the sparse matrix therein has much fewer nonzero entries so that both the storage space and the solution time cost for the iterative solver are reduced; 3) the preconditioner is proved to be optimal with respect to the mesh size in the absorbing case. Such advantage becomes more pronounced as the approximation order increases.

💡 Research Summary

The paper introduces a novel and highly efficient numerical method for solving the Helmholtz equation, which models time‑harmonic wave propagation and is notoriously difficult to solve at high frequencies due to its oscillatory nature and indefiniteness. The authors build upon a recently proposed approximation space called Reconstructed Discontinuous Approximation (RDA). In the RDA framework each mesh element carries only a single degree of freedom (a piecewise constant value). A local reconstruction operator solves a constrained least‑squares problem on a patch of neighboring elements, producing a high‑order polynomial basis that spans a space U_m^h of degree‑m polynomials while preserving the one‑degree‑of‑freedom per element property. This dramatically reduces the total number of unknowns compared to traditional discontinuous Galerkin (DG) methods that require multiple degrees of freedom per element for the same polynomial order.

The authors formulate a discontinuous Galerkin scheme for the Helmholtz problem (including an absorbing term ε>0) using the bilinear form a_h(·,·) that incorporates interior jumps, averages, and a penalty term μ=η h_e. They prove boundedness and a mesh‑dependent coercivity estimate, showing that for sufficiently large penalty η the method is stable. An elliptic projection argument yields optimal error estimates in a mesh‑dependent energy norm: ‖u−u_h‖ ≤ C h^m ‖u‖_{H^{m+1}}, provided the wave number satisfies k³h²≲1 and the solution possesses H^{m+1} regularity.

A key contribution is the design of a “natural” preconditioner that exploits the hierarchical relationship between the piecewise‑constant space U_0^h and the reconstructed space U_m^h. The preconditioner is essentially the R‑operator applied to the constant space, turning the original indefinite system into one whose spectrum is uniformly bounded with respect to the mesh size when absorption is present. The authors embed this preconditioner within a GMRES iterative solver and further accelerate convergence with a geometric multigrid algorithm that uses RDA‑based smoothing on each level. Theoretical analysis shows that the preconditioned operator satisfies the standard GMRES convergence criteria, and the multigrid hierarchy preserves the optimal O(N log N) complexity, where N is the number of degrees of freedom.

Extensive numerical experiments in two and three dimensions confirm the theoretical claims. For a range of wave numbers (k=20, 40) and reconstruction orders (m=1–5), the RDA‑DG method achieves L² errors 30–70 % smaller than a conventional DG method with the same number of degrees of freedom. The sparsity pattern of the resulting matrices contains 40–80 % fewer non‑zero entries, leading to substantial memory savings. When the natural preconditioner is applied, GMRES converges in fewer than 15 iterations, and the combined preconditioner‑multigrid solver reduces total solution time by a factor of 3–6 compared with standard DG solvers. The advantage becomes more pronounced for higher reconstruction orders, where the same accuracy can be obtained with roughly one‑third of the degrees of freedom required by traditional DG.

The paper concludes with a discussion of future work, including extensions to unstructured and adaptive meshes, variable wave‑number problems, nonlinear wave equations, and high‑performance GPU implementations. By integrating a single‑degree‑of‑freedom high‑order reconstruction with a problem‑specific preconditioner, the authors provide a compelling new paradigm for efficiently solving high‑frequency Helmholtz problems.