¹

1 Introduction↩︎

Helmholtz equation and Variational Principles↩︎

The Helmholtz equation with an impedance boundary condition \[\label{Helm95PDE} -\Delta u - k^2 u = f\quad \text{ in } \varOmega \qquad \partial_{\mathbf{n}} u = i k u \quad \text{ on } \partial\varOmega\tag{1}\] is the time-harmonic analogue of the wave equation with absorbing boundary conditions of first order. Here, the function \(u\) is complex-valued, \(u:\varOmega\to\mathbb{C}\), \(\varOmega\subset\mathbb{R}^{\nu }\) is a bounded Lipschitz domain, \(f\in L^2(\varOmega)\) is a complex-valued function and \(k>0\) is the wavenumber. We use the compact notation \({\mathscr{L}}u:=-\Delta u - k^2 u\). The use of various absorbing boundary conditions is motivated by the need to reduce wave propagation problems posed on large or unbounded domains to problems defined on computationally tractable bounded domains. While this model is classical and broadly applicable, the design of numerical algorithms capable of accurately approximating it at high wavenumbers remains a significant challenge.

In this work, we investigate whether variational principles can be associated with the Helmholtz equation and in particular with 1 . Specifically, we ask whether there exist energies motivated from first principles whose stationary points, the zeros of their first variation, coincide with solutions of 1 . To the best of our knowledge, such principles have not been systematically developed for impedance Helmholtz problems such as 1 .

Starting from Hamilton’s principle for the wave equation, we derive time-harmonic energies that are, in general, indefinite. These energies apply naturally to the Helmholtz equation posed on unrestricted domains. When the domain is truncated and impedance boundary conditions are imposed, suitable modifications become necessary. In unrestricted domains, the corresponding Lagrangian functionals may be chosen to be real-valued. However, to account for the influence of the artificial boundary introduced by the absorbing condition, we employ complex-valued Lagrangians. Nevertheless, the variational structure is preserved, since the stationarity condition has a unique critical point corresponding to the solution of 1 .

We then construct strongly coercive augmentations of these energies. These regularised Lagrangians are carefully designed modifications of the original indefinite functionals and provide stability bounds in the case of star-shaped domains.

For all formulations considered, we address existence and uniqueness questions for the corresponding variational problems by employing the complex version of the Lax-Milgram theorem. Throughout, our derivations and proposed formulations are motivated by the goal of remaining as close as possible to the underlying physical properties of the problem.

Apart from its theoretical appeal, the derivation of new variational principles has significant practical implications. Modern methods that utilise neural-network-based discrete spaces typically rely on discrete minimisation problems, where the loss functional is constructed either by approximating a continuous energy associated with the problem or by using the \(L^2\) residual of the differential operator Karniadakis:pinn:orig:2019?. While the residual-based approach is straightforward, it is desirable to employ alternative methods, especially for PDEs with highly oscillatory behaviour such as 1 . The hope is that, by employing a variational principle with strong coercivity, one obtains a more robust foundation for efficient and reliable methods. In fact, we show that the proposed variational principle may lead to an effective neural-network formulation for 1 . We also present prototype numerical evidence that the same principle leads to efficient finite element formulations and outline how its properties motivate discretisations with improved characteristics. Strongly coercive bilinear forms for the Helmholtz equation were previously developed in Moiola_Spence2014?, leading to finite element methods that, in general, produce non-symmetric matrices. While we employ similar analytical tools, namely the Rellich and Morawetz identities, which are effective in the analysis of 1 , our motivation and design strategy differ substantially, see Remark [remark95MS] and Section ¿sec:sec:design?.

Bibliography and Main Contributions↩︎

The Helmholtz equation has been extensively investigated from multiple perspectives, including the treatment of various boundary conditions and alternative formulations; see, for example MorawetzLudwig1968?, FixMarin1978?, Goldstein1982?, BGT1985?, Aziz_et_al_1988?, KellerGivoli1989?, Givoli1992?, Douglas_et_al_1993?, perthame1999morrey?, Yafaev2000?, SingerTurkel2004?, CummingsFeng2006?, Mitsoudis2007?, Athanassoulis_et_al_2008?, chandler2008wave?, MitsoudisPlexousakis2009?, MMP2012?, SauterTorres2018?, GROTE_control_H_2019?, GrahamSauter2020? and references therein. Frequency-explicit stability estimates are important because, among other reasons, when combined with finite element error bounds, they highlight the constraints imposed by moderate to high wavenumbers and mesh resolution. These analyses reveal subtle correlations between the wavenumber and the mesh size, an aspect first systematically investigated within the research programme initiated by I. Babuška and his collaborators (see, e.g., IhlenburgBabuska1995?, Ihlenburg1998?, BabuskaSauter2000?, and the references therein). Stability analyses based on energy methods employing nonstandard test functions were developed in Makridakis_et_al_1996? and Melenk1995?. These test functions belong to a broader class of so-called Morawetz multipliers, introduced by Morawetz and Ludwig MorawetzLudwig1968? for the justification of geometric optics for the reduced wave equations. Such multipliers are instrumental in deriving Rellich-type identities, as used, for instance, in CummingsFeng2006?, Spence_Chandler-Wilde_integral_2011?, Moiola_Spence2014? and related works. The multipliers \(\mathcal{M}\) typically involve \(x\cdot\nabla u\) together with possible lower-order terms; the identities obtained by testing the PDE with \(\mathcal{M}u\) are known as Morawetz identities. These identities play a crucial role in establishing stability estimates for a wide range of problems in both linear and nonlinear wave propagation, Morawetz_PIrish_1972?, Lesky2003?, MetcalfeEtal2005?, Tao_2006?, DAncona_2012?, Hintz_Zworski_2017?, Tao_MI_notices?, cossetti2024virial?, Bigniardi_Moiola_wave_2025?. Such identities will also be instrumental in the design and analysis of the variational principles introduced in this paper.

The main contributions of this paper are as follows. We derive variational principles for the Helmholtz equation 1 from first principles, starting from Hamilton’s principle (least action) and a Lagrangian formulation of the wave equation to motivate the associated energy functionals. As expected, these functionals are generally indefinite. Incorporating impedance boundary conditions presents an additional challenge, which we address by employing complex-valued Lagrangians. We then construct regularised versions of these functionals that are strongly coercive. Our main idea is to perturb the physical Lagrangian by adding a term of the form \(\gamma \| {\mathscr{L}}u-f \|^2\) and a corresponding boundary term. This regularisation enhances the stability of the variational principle, though in a nontrivial way. It turns out that, by choosing \(\gamma\) suitably, one can obtain a coercive functional with respect to an appropriately strengthened norm. The resulting coercivity arises from a combination of the original physical Lagrangian and suitable Morawetz identities, which are employed to derive lower bounds first for the bulk energy terms and subsequently for the boundary contributions. The resulting energies form the basis for well-posed variational principles and minimisation problems whose unique minimisers approximate the solutions of the Helmholtz equation.

Finally, we show how these variational principles yield discrete approximation schemes. We present computational methods based on classical conforming finite elements and neural-network discretisations. The primary focus is the mathematical analysis underpinning the design of the new variational principles, but we also include prototype formulations for both finite elements and neural networks to illustrate their potential for computing approximations to 1 . It is worth emphasising that neural network approximations depend on appropriate minimisation principles, and, to our knowledge, this is a new coercive augmented variational formulation for impedance Helmholtz problems that can be used with neural trial spaces. Preliminary results suggest several desirable properties, a detailed numerical investigation is deferred to future work. Regarding the extensive literature on finite element schemes for approximating the Helmholtz equation, as mentioned above, analyses combining frequency-explicit stability bounds for weak solutions provide both theoretical and computational insights into the behaviour of such methods and the pollution effects that become particularly significant at high wavenumbers BabuskaSauter2000?. Classical error estimates and quasi-optimality bounds based on inf-sup conditions trace back to ideas introduced by Aziz et al. Aziz_et_al_1988?, Makridakis et al. Makridakis_et_al_1996? and Melenk Melenk1995?, which build upon the approach of Schatz Schatz1974? for indefinite problems. Further detailed studies of the pollution effect and the improved performance of higher-order methods can be found in Sauter_refined?, MelenkSauter2010?. Least squares methods were considered in First-orderLS_Helm_Lee_et_al_2000?, methods employing basis functions that solve local problems and thus achieve enhanced performance are discussed in HipMP_serveyT_2016?. Discontinuous Galerkin methods were investigated in Feng_dg_2011?, while specially designed algorithms that avoid the pollution effect were proposed in Peterseim_polloutionHelm_2017?. Extensions to problems with variable coefficients have also been considered, for example, in GrahamSauter2020?.

The paper is organised as follows. Section 2 recalls the physical derivation of the wave equation via Hamilton’s principle and motivates natural Lagrangian functionals for the Helmholtz equation. Section 3 introduces notation and functional-analytic preliminaries and defines the variational formulations and the regularised Lagrangians, and Section 4 proves the main results. Section 5 discusses discrete approximations based on neural networks and conforming finite elements and illustrates their formulation within the variational framework.

2 Energy principles↩︎

The basic model: Transverse vibrations of a string↩︎

In this section, we explain how the wave and the Helmholtz equations can be derived from Hamilton’s principle. Such a derivation provides us with a natural Lagrangian, which in turn arises in the standard variational reformulation of certain boundary value problems for the Helmholtz equation. For simplicity, we restrict the discussion to the case \(\nu = 1\), and we consider the transverse motion of a string occupying the interval \(\Omega=(\alpha,\beta)\) along the \(x\)-axis. Analogous results can be derived for \(\nu = 2\) or \(\nu = 3\), corresponding to the transverse motion of a membrane or to sound waves in an ideal compressible fluid, respectively. We assume that the mass density (per unit length) of the string is \(\rho(x)\) and that the string is subjected to tension by opposite forces of strength \(S\) at its ends. Then, the transverse deflection \(U(x,t)\) satisfies the wave equation \[\label{eq:wave-nh} \frac{1}{c^2(x)}\frac{\partial^2 U}{\partial t^2}=\frac{\partial^2 U}{\partial x^2} {- \frac{F(x,t)}{S}}, \;\; t>0\;, \;\; \alpha < x< \beta,\tag{2}\] where \(c(x)= \left(S/\rho(x)\right)^{1/2}\) is the phase velocity of the waves propagating along the string, and \(F(x,t)\) denotes an external transverse load per unit length (the sign convention in 2 is immaterial for what follows). We assume that \(c(x)\) and \(F(x,t)\) are sufficiently smooth functions.

At the end points \(x=\alpha,\beta\), it is necessary to impose boundary conditions, usually in the form of linear combination of \(U(x,t)\) and \(\partial U/\partial x (x,t)\), for all \(t>0\), and also we must prescribe initial data for \(U(x,0)\) and \(\partial U/\partial t(x,0)\), for all \(\alpha < x< \beta\), in order to be able to construct a unique solution of the problem.

In many applications, as well as for fundamental theoretical reasons in scattering theory, it is necessary to study the corresponding time-harmonic problem, and, in particular, the associated Green’s function. To this end, we consider that the motion of the string starts from rest at time \(t = 0\), and it subsequently evolves under the action of a transverse force distribution which varies periodically in time with frequency \(\omega\). Such a force is represented by the distribution \(F(x,t)=\Re\bigl(e^{-i\omega t} f(x)\bigr) H(t)\), where, typically, \(f(x)\) is a localised function and \(H(t)\) is the Heaviside function. If the Limiting Amplitude Principle applies, then as \(t \to \infty\) the motion becomes time-harmonic of the form \(U(x,t) = \Re\bigl(u(x) e^{-i\omega t}\bigr)\), where the complex amplitude \(u(x)\) satisfies the Helmholtz equation \[\label{eq:helm-nh} u''(x) + k^2(x) u(x) = \frac{f(x)}{S}, \qquad k(x) =\frac{\omega}{c(x)}, \quad \alpha < x< \beta \;,\tag{3}\] and certain boundary conditions dictated by \(U(x,t)\) and \(\partial U/\partial x (x,t)\) at \(x=\alpha,\beta\). Whenever \(\alpha=-\infty\) and/or \(\beta=+\infty\), the boundary conditions must be replaced by outgoing (Sommerfeld-type) radiation conditions to ensure uniqueness of the solution.

Variational derivation of the wave equation↩︎

Let \[\mathcal{T}(x,t)= \frac{1}{2} \rho(x) \Bigl(\frac{\partial U}{\partial t}\Bigr)^2\] be the kinetic energy density (that is, energy per unit length), and \[\mathcal{V}(x,t)= {\frac{1}{2} S \Bigl(\frac{\partial U}{\partial x}\Bigr)^2}\] be the potential energy density of the string ; the external loading enters through the work term in the Lagrangian rather than as stored potential energy. For any time interval \(0<t<T\), \(T<\infty\), the action of the system is given by the space-time integral \[\label{eq:dvar} \mathcal{I}(U)= \frac{1}{T}\int_{0}^{T}\int_{\alpha}^{\beta} \mathcal{L}\Bigl(x, {U},\frac{\partial U}{\partial t},\frac{\partial U}{\partial x}\Bigr) {\mathrm d}x {\mathrm d}t \;,\tag{4}\] of the Lagrangian density \[\label{eq:lagrdens} \begin{align} \mathcal{L}\Bigl(x, {U},\frac{\partial U}{\partial t},\frac{\partial U}{\partial x}\Bigr) &= \mathcal{T}(x,t)-\mathcal{V}(x,t) { - F(x,t) U(x,t)} \\ &= \frac{1}{2} \rho(x) \Bigl(\frac{\partial U}{\partial t}\Bigr)^2 - \frac{1}{2} S \Bigl(\frac{\partial U}{\partial x}\Bigr)^2 { - F(x,t) U(x,t)} \;. \end{align}\tag{5}\] According to Hamilton’s variational principle Lanczos_VP_1949?, Gelfand_Fomin_CV?, the actual motion of the string follows from the stationarity of the action, i.e., the derivative of the action with respect to the field \(U\) in the direction of any appropriate admissible function \(\eta=\eta(x,t)\) Evans?, \[\label{eq:deul1} {D_{U}\, \mathcal{I}{(U)}[\eta]=\frac{d}{d\epsilon}I(U+\epsilon \eta)|_{\epsilon=0}=0} \;\tag{6}\]

Equation 6 implies \[\label{eq:lagrange} {D_{U}\, \mathcal{I}(U)}[\eta]=\frac{1}{T} \, \int_{0}^{T}\int_{\alpha}^{\beta} \frac{\delta {\mathcal{L}}}{\delta U}\, \eta \, dx \, dt + \frac{1}{T}\,\int_{0}^{T}\int_{\alpha}^{\beta} \, \frac{\partial}{\partial t} \Biggl(\frac{\partial {\mathcal{L}}}{\partial \Bigl(\frac{\partial U}{\partial t}\Bigr)} \, \eta \Biggr)+ \frac{\partial}{\partial x} \Biggl(\frac{\partial {\mathcal{L}}}{\partial \Bigl(\frac{\partial U}{\partial x}\Bigr)} \, \eta \Biggr) \, dx \, dt =0\;,\tag{7}\] where\(\frac{\delta {\mathcal{L}}}{\delta U}\) is the variational derivative \[\label{eq:varder} \frac{\delta {\mathcal{L}}}{\delta U}= \, \frac{\partial {\mathcal{L}}}{\partial U} -\frac{\partial }{\partial t}\Biggl( \frac{\partial {\mathcal{L}}}{\partial \Bigl( \frac{\partial U}{\partial t}\Bigr)} \Biggr) - \frac{\partial }{\partial x}\Biggl( \frac{\partial {\mathcal{L}}}{\partial \Bigl(\frac{\partial U}{\partial x} \Bigr)} \Biggr) \;,\tag{8}\] For the Lagrangian density 5 , equation 7 gives the wave equation \[\label{eq:wave} \frac{1}{c^2(x)} \frac{\partial^2 U}{\partial t^2}-\frac{\partial^2 U}{\partial x^2}= {- \frac{F}{S}} \;,\tag{9}\] assuming that the boundary term \[\label{eq:deul2} {\mathcal{R}}=\frac{1}{T}\,\int_{0}^{T}\int_{\alpha}^{\beta} \, \Biggl\{ \frac{\partial}{\partial t} \Biggl(\frac{\partial {\mathcal{L}}}{\partial \Bigl(\frac{\partial U}{\partial t}\Bigr)} \, \eta\Biggr)+ \frac{\partial}{\partial x} \Biggl(\frac{\partial {\mathcal{L}}}{\partial \Bigl(\frac{\partial U}{\partial x}\Bigr)} \, \eta \Biggr)\Biggr\}\, dx\, dt\tag{10}\] vanishes.

By integrating by parts, the boundary term 10 is written as follows \[\begin{align} \label{eq:remainder1} { \mathcal{R}}=\frac{1}{T}\, \int_{\alpha}^{\beta}\rho \Biggl( U_t(x,T)\eta(x,T)-U_t(x,0) \eta(x,0)\Biggr) \, dx \nonumber \\ -S\, \frac{1}{T}\,\int_{0}^{T}\Biggl(U_x(\beta, \tau)\eta(\beta, \tau)-U_x(\alpha, \tau)\eta(\alpha, \tau)\Biggr)\, d\tau \;. \end{align}\tag{11}\]

If we choose \(\eta\) to have compact support in \(G=(\alpha, \beta)\times(0,T)\), then \(\mathcal{R}\equiv 0\), and 9 holds in the sense of distributions on \(G\). Formally, since in this case we do not vary \(U\) on \(\partial\, G\), we can handle inhomogeneous initial data, as well as impedance boundary conditions of the form \(U_{t}=\pm \mu \, U_{x}\) for some \(\mu >0\). Moreover, if we assume that the admissible function \(\eta\) satisfies \(\eta(x,0)=\eta(x,T)=0\) for all \(\alpha \le x \le \beta\), that is we do not vary \(U(x,t)\) at the initial and final time, then we must impose at \(x=\alpha \;, \beta\), the natural Neumann (zero slope) boundary condition \(\frac{\partial U}{\partial x}=0\), to ensure that \({\mathcal{R}}=0\).

Variational derivation of the Helmholtz equation↩︎

The Helmholtz equation can be derived from the time-harmonic version of Hamilton’s principle. We choose \(T=\frac{2\pi}{\omega}\), the period of the time-harmonic wave, and define the stationary density \(\ell\) as the time average of the time-dependent Lagrangian density \(\mathcal{L}\), \[\label{eq:meanlang} \ell\Bigl(x,u,\tfrac{\partial u}{\partial x}\Bigr) := \frac{2}{T}\int_{0}^{T} \mathcal{L}\Bigl(x, {U},\frac{\partial U}{\partial t},\frac{\partial U}{\partial x}\Bigr) {\mathrm d}t \;,\tag{12}\] where we substitute the time-harmonic ansatz \(U(x,t)=\Re\bigl(u(x) e^{-i\omega t}\bigr)\), and \(F(x,t)=\Re\bigl(f(x) e^{-i\omega t}\bigr)\).

After time integration in 12 , we derive \[\label{eq:hlagrdens} \ell\Bigl(x,u,\tfrac{\partial u}{\partial x}\Bigr) = {\frac{1}{2} \, S\, \Bigl(k(x)^2 |u|^2- \bigl|\tfrac{\partial u}{\partial x}\bigr|^2\Bigr) - \Re\bigl(f(x) \overline{u(x)}\bigr),} \qquad k(x)=\frac{\omega}{c(x)}.\tag{13}\] (Any positive constant multiple of \(\ell\) produces the same stationary points, but the coefficients in 13 match the time average of 5 .)

The Helmholtz equation is derived from the stationarity of the action \[\mathcal{J} (u) := \int_{\alpha}^{\beta} \ell\Bigl(x,u,\tfrac{\partial u}{\partial x}\Bigr) {\mathrm d}x .\] This action arises naturally by rewriting the action \(\mathcal{I}\) in 4 in terms of the time average 12 , \[\begin{align} \label{eq:dvarh} \mathcal{I} = \frac{1}{T}\int_{0}^{T}\int_{\alpha}^{\beta} \mathcal{L}\Bigl(x, {U},\frac{\partial U}{\partial t},\frac{\partial U}{\partial x}\Bigr) {\mathrm d}x {\mathrm d}t = \frac{1}{2} \,\int_{\alpha}^{\beta} \ell\Bigl(x,u,\tfrac{\partial u}{\partial x}\Bigr) {\mathrm d}x . \end{align}\tag{14}\]

The use of the action \(\mathcal{J}\) for the variational derivation of the Helmholtz equation is not standard, and it has been inspired by the variational derivation of the Schrödinger equation from a real Lagrangian Mo_book?, since both the Schrödinger wave function and the wave function \(u\) in the Helmholtz equation are complex.

The time-harmonic version of Hamilton’s principle dictates that \[\label{eq:svar951} {D_{u}\, \mathcal{J}{(u)}[\theta]=\frac{d}{d\epsilon}\mathcal{J}{(u)}(u+\epsilon \theta)|_{\epsilon=0}=0} \;.\tag{15}\] For an admissible function \(\theta=\theta(x)\) with compact support in \((\alpha,\beta)\), this condition implies that the Helmholtz equation 3 holds in the distributional sense. Similarly to the variational derivation of the wave equation, for other admissible functions \(\theta(x)\) we can accommodate various boundary conditions.

With variations taken with respect to the real inner product on \(H^1((\alpha,\beta);\mathbb{C})\), the Euler-Lagrange equation for 15 is precisely the Helmholtz equation in the sense of distributions. See Section 3 for the detailed computation

Impedance boundary conditions↩︎

The total energy density (per unit length) of the string at time \(t\) is \[\label{eq:endens} \mathcal{E}(x,t)= \mathcal{T}(x,t)+\mathcal{V}(x,t) .\tag{16}\] For \(F\equiv 0\), the total energy \(E(t)=\int_{\alpha}^{\beta} \mathcal{E}(x,t) {\mathrm d}x\) is constant for suitable choices of initial and boundary conditions. In fact, by direct calculation and using 2 , \[\label{eq:dere} \frac{{\mathrm d}E(t)}{{\mathrm d}t} = \int_{\alpha}^{\beta}\Biggl( \frac{1}{2} \rho(x) \frac{\partial}{\partial t}\Bigl(\frac{\partial U}{\partial t}\Bigr)^2 + \frac{1}{2} S \frac{\partial}{\partial t}\Bigl(\frac{\partial U}{\partial x}\Bigr)^2 \Biggr) {\mathrm d}x = S \Biggl[\frac{\partial U}{\partial x} \frac{\partial U}{\partial t}\Biggr]_{x=\alpha}^{x=\beta} .\tag{17}\] If, for example, either the deflection \(U\) or the slope \(\partial U/\partial x\) at the endpoints is zero for all \(t>0\), then the last term in 17 vanishes. For an infinite string, where \(\alpha=-\infty\) and \(\beta=+\infty\), the derivative \(\partial U/\partial t\) vanishes at both ends due to the finite speed of propagation.

However, the boundary conditions \[\label{eq:time-imp} \frac{\partial U}{\partial x}=-\mu \frac{\partial U}{\partial t}\quad \text{at } x=\beta, \qquad \frac{\partial U}{\partial x}=+\mu \frac{\partial U}{\partial t}\quad \text{at } x=\alpha,\tag{18}\] where \(\mu>0\), are dissipative, since by substituting 18 into 17 we get \[{\frac{{\mathrm d}E(t)}{{\mathrm d}t} = -\mu\, S\bigl((\partial_t U)^2(\beta,t)+(\partial_t U)^2(\alpha,t)\bigr)< 0,}\] unless \(\partial_t U(\alpha,t)=\partial_t U(\beta,t)=0\). In the sequel we will take the string tension force \(S=1\).

For the time-harmonic problem, the boundary conditions 18 imply the following dissipative (impedance) boundary conditions for the Helmholtz equation: \[\label{eq:time-harm-imp} \frac{{\mathrm d}u}{{\mathrm d}x}= i \sigma_{\beta} k_{\beta} u \quad \text{at } x=\beta, \qquad \frac{{\mathrm d}u}{{\mathrm d}x}= - i \sigma_{\alpha} k_{\alpha} u \quad \text{at } x=\alpha,\tag{19}\] where \(k_{\alpha}=\omega/c(\alpha)\), \(k_{\beta}=\omega/c(\beta)\), and \(\sigma_{\alpha}=\mu c(\alpha)\), \(\sigma_{\beta}=\mu c(\beta)\).

In higher dimensions the Helmholtz equation is stated in a bounded domain \(\Omega\subset\mathbb{R}^{\nu }\), where \(u\) represents the deflection of a membrane for \({\nu }=2\), and, e.g., the pressure of a barotropic fluid for \({\nu }=3\). In these cases the impedance condition on the boundary has the form \[{\partial_{\mathbf{n}} u} = i k \sigma u \quad \text{on } \partial\Omega,\] with \(\mathbf{n}\) the outward unit normal and \(\sigma>0\). The canonical choice \(\sigma\equiv 1\) yields the boundary condition used in 1 . The impedance boundary conditions cannot be derived from Hamilton’s variational principle associated with the real Lagrangian 13 , since real Lagrangians generate conservative systems and therefore do not capture dissipative boundary effects.

However, this obstruction can be remedied formally by adding a purely imaginary boundary term to the real action \(\mathcal{J}=\int_{\Omega} \ell {\mathrm d}x\) thereby obtaining a complex action. This additional term models the boundary \(\partial\Omega\) as an energy-dissipating membrane, see FR2015?, where boundary conditions are interpreted as physical processes occurring within thin boundary layers.

To proceed, we first make the following observation. We write the Lagrangian \(\ell\) in the equivalent form \[\label{eq:mdhlagrdens} \ell\Big(x, u, \overline{u}, \nabla u, \nabla \overline{u}\Big) := \frac{1}{2} \Big( \nabla u \cdot \nabla \overline{u}- k^2 u\overline{u}\Big) - \Re(f\overline{u}),\tag{20}\] and consider the action \(\mathcal{J}=\int_{\Omega} \ell {\mathrm d}x\) as a functional of the complex-valued functions \(u\) and \(\overline{u}\) (equivalently, of \(\Re u\) and \(\Im u\)). Following the variational formalism used in the derivation of the Schrödinger equation from a real Lagrangian Mo_book?, we vary \(u\) and \(\overline{u}\) independently to derive simultaneously the Helmholtz equation and its complex conjugate.

Then we construct the complex action \[\mathcal{S}(u\;, \overline{u}) := \int_{\Omega} \ell\Bigl(x, u, \overline{u}, \nabla u, \nabla \overline{u}\Bigr) {\mathrm d}x - \frac{ik}{2}\int_{\partial \Omega}u \overline{u}{\mathrm d}S \;,\] and compute the Wirtinger derivative with respect to \(\overline{u}\) in the direction \(\overline{v}\), \[\label{eq:wirtubar} D_{\overline{u}}\mathcal{S}(u\;, \overline{u}) := \frac{d}{d\epsilon}\mathcal{S}(u, \overline{u}+\epsilon \overline{v})\Big|_{\epsilon=0} \;.\tag{21}\] We obtain \[\label{Section2:complexVF} D_{\overline{u}}\mathcal{S}(u\;, \overline{u})= \frac{1}{2} \int_\Omega \Bigl(\nabla u\cdot \nabla \overline{v}- k^2 \, u\overline{v}\Bigr)\, {\mathrm d}x -\frac{1}{2} \int_\Omega f \overline{v}\, {\mathrm d}x -\frac{ i\, k}{2}\int_{\partial \Omega} u \overline{v}\, {\mathrm d}S \, .\tag{22}\] By imposing the stationarity condition \(D_{\overline{u}}\mathcal{S}(u\;, \overline{u})=0\) and integrating by parts in the first integral, we recover both the Helmholtz equation for \(u\) in \(\Omega\) and the impedance boundary condition on \(\partial \Omega\). This stationarity condition is equivalent to the weak formulation of the impedance boundary value problem, see 25 with \(\zeta_\Omega =f\) and \(\eta_{\partial\Omega} =0\) in the next section. It is important to note that, when we impose the stationarity condition \(D_{u}\mathcal{S}(u\;, \overline{u})=0\), we recover the Helmholtz equation for \(\overline{u}\) in \(\Omega\) together with the conjugate impedance boundary condition \({\partial_{\mathbf{n}} \overline{u}} = -i k \overline{u}\) on \(\partial \Omega\). This distinction arises from the fact that the problem is non-selfadjoint.

3 Variational principles for the Helmholtz equation↩︎

Notation and preliminaries↩︎

We introduce some standard notation. Throughout, \(\Omega\subset\mathbb{R}^{\nu }\) is a bounded Lipschitz domain with outward unit normal \(\mathbf{n}\) Grisvard_book?. For a domain \(D\), we write \((\cdot,\cdot)_D\) for the \(L^2\) inner product on \(D\), \(\|\cdot\|_D\) for the corresponding \(L^2\) norm, and \(\|\cdot\|_{m,D}\) for the Sobolev norm on \(H^m(D)\). When \(D=\Omega,\) we omit the subscript \(D.\) On the boundary we use the trace spaces \(H^s(\partial\Omega)\), \(s\in\{\tfrac32,\tfrac12,-\tfrac12,-\tfrac32\}\), and the trace operator \(\tau :H^1(\Omega)\to H^{1/2}(\partial\Omega)\) Lions_Magenes_I?. When no confusion may arise we will use just \(u\) to denote \(\tau u\, .\)

It will be useful to consider a generalised impedance problem on \(\Omega\): seek a complex-valued function \(w\) such that \[\label{eq:gen-Helm} - \Delta w - k^2 w = \zeta_\Omega \quad \text{in } \Omega, \qquad \partial_{\mathbf{n}} w - i k w = \eta_{\partial\Omega} \quad \text{on } \partial\Omega,\tag{23}\] where the source terms \(\zeta_\Omega\) and \(\eta_{\partial\Omega}\) are given. For the minimal weak setting we take \(\zeta_\Omega\in H^{-1}(\Omega)\) and \(\eta_{\partial\Omega}\in H^{-1/2}(\partial\Omega)\). In several places below we also work in an \(L^2\)-based setting, writing \(\zeta_\Omega\in L^2(\Omega)\) and \(\eta_{\partial\Omega}\in L^2(\partial\Omega)\) when stronger regularity is required by the functionals.

We denote by \({\mathscr H}\) the space \(H^1(\Omega)\) equipped with the wavenumber-dependent norm \[\|v\|_{H^1_k(\Omega)} := \Bigl(\|\nabla v\|_\Omega^2 + k^2 \|v\|_\Omega^2\Bigr)^{1/2}\, .\]

For energies that involve the bulk residual \({\mathscr{L}}v:=-\Delta v-k^2 v\) in \(L^2(\Omega)\) we use the space \[{\mathscr V}:= \Bigl\{ v\in H^1(\Omega) : {\mathscr{L}}v \in L^2(\Omega),\;v\in H^1(\partial \Omega), \text{and}\;\partial_{\mathbf{n}} v - i k v \in L^2(\partial \Omega ) \Bigr\},\] endowed with the norm \[\|v\|_{{\mathscr V}}^2 := \|\nabla v\|_\Omega^2 + k^2 \|v\|_\Omega^2 + \|{\mathscr{L}}v\|_\Omega^2 + k ^2 \|v \|_{\partial \Omega} ^2 + \| \nabla v \|_{\partial \Omega} ^2 + \| \partial_{\mathbf{n}} v - i k v \|_{\partial \Omega} ^2 .\] The homogeneous-impedance subspace is \[{\mathscr V}_{BC}:= \Bigl\{ v\in {\mathscr V}: \partial_{\mathbf{n}} v - i k v = 0 \;\text{in } L^2(\partial\Omega) \Bigr\},\] and in the sequel we shall use the notation \[\|v\|_{\mathscr {V}_{BC}}^2:=\|\nabla v\|_\Omega^2 + k^2 \|v\|_\Omega^2 + \|{\mathscr{L}}v\|_\Omega^2 + k ^2 \|v \|_{\partial \Omega} ^2 + \| \nabla v \|_{\partial \Omega} ^2 .\]

The regularity and other properties of 23 can be formulated with the aid of its variational formulation. To this end, consider the sesquilinear form \({\mathcal{B}}\) defined by \[\label{eq:bf} {\mathcal{B}}(u,v) := \int_\Omega \nabla u\cdot \nabla \overline{v}\, {\mathrm d}x - k^2 \int_\Omega u\overline{v}\, {\mathrm d}x - i\, k\int_{\partial \Omega} u \overline{v}\, {\mathrm d}S \, .\tag{24}\] Given \(\zeta_\Omega\) and \(\eta_{\partial\Omega}\) as above, the weak problem is:

Find \(w\in {\mathscr H}\) such that \[\label{eq:wf-gen} {\mathcal{B}}(w,v) = (\zeta_\Omega,v)_\Omega + {\langle \eta_{\partial\Omega}, v \rangle_{H^{-1/2}, H^{1/2}}} \qquad \text{for all } v\in {\mathscr H}.\tag{25}\] In the \(L^2\)-based setting the duality pairing reduces to the \(L^2(\partial\Omega)\) inner product.

By construction, \({\mathscr V}\) is the class of functions for which the data-to-solution relation in 25 is meaningful with \(\zeta_\Omega\in L^2(\Omega)\) and \(\eta_{\partial\Omega}\in L^2(\partial\Omega)\), and \({\mathscr V}_{BC}\subset{\mathscr V}\) corresponds to the case \(\eta_{\partial\Omega}=0\). The solution \(u\) of 1 satisfies 25 with \(\zeta_\Omega=f\) and \(\eta_{\partial\Omega}=0\), hence \(u\in{\mathscr V}_{BC}\), see Moiola_Spence2014? and Remark [Rem3.4].

Variational Principles↩︎

In Section 2, we show that the physical energy associated with the Helmholtz equation in unrestricted domains is \[\label{Energy95H95main95s2} \mathcal{E}_P(v) = \int_{\Omega } \frac{1}{2} \Bigl( |\nabla v |^2 - k^2|v|^2 \Bigr) {\mathrm d}x - \Re \int_{\Omega } f \overline{v} {\mathrm d}x .\tag{26}\] Then the map \(\mathcal{E}_P : {\mathscr H}\to {\mathbb{R}}\) is clearly differentiable. We take variations with respect to the real inner product on \(H^1(\Omega;\mathbb{C})\), i.e. \(\langle \phi,\psi\rangle_{\mathbb{R}}:=\Re \int_{\Omega}\phi \overline{\psi} {\mathrm d}x\), so that \(D\mathcal{E}_P(u)[v]\in\mathbb{R}\) for all \(u,v\). To find its derivative we first notice \[\label{Energy95H95main95Re} \mathcal{E}_P(v) = \Re \Bigl\{ \int_{\Omega } \frac{1}{2} \bigl( |\nabla v |^2 - k^2|v|^2 \bigr) {\mathrm d}x - \int_{\Omega } f \overline{v} {\mathrm d}x \Bigr\}.\tag{27}\] Then, \[\label{Energy95H95main95Re2}\begin{align} {\langle } D\mathcal{E}_P(u) , v {\rangle } =& \Re \Bigl\{ \int_{\Omega } \tfrac12 \Bigl( \nabla u\cdot \nabla \overline{v} + \nabla v\cdot \nabla \overline{u}- k^2( u \overline{v}+ v \overline{u}) \Bigr) {\mathrm d}x - \int_{\Omega } f \overline{v} {\mathrm d}x \Bigr\}\\[2pt] = & \Re \Bigl\{ \int_{\Omega } \bigl( \nabla u\cdot \nabla \overline{v}- k^2 u \overline{v}\bigr) {\mathrm d}x - \int_{\Omega } f \overline{v} {\mathrm d}x \Bigr\}. \end{align}\tag{28}\] Consider \(u= u_R + i u_I\), \(v= v_R + i v_I\), and \(f = f_R + i f_I\). Then stationary points of \(\mathcal{E}_P,\) i.e., \(u\in {\mathscr H}\) satisfying, \[\label{Energy95H95main95stat95p}\begin{align} {\langle } D\mathcal{E}_P(u) , v {\rangle } = 0, \quad \text{for all } \;v \in {\mathscr H}, \end{align}\tag{29}\] are such that, for any \(v\in {\mathscr H},\) \[\label{Energy95H95main95Re95system} \begin{align} & \int_{\Omega } \Bigl( \nabla u _R \cdot \nabla v_R \;- k^2 u_R v_R \Bigr) {\mathrm d}x - \int_{\Omega }{ f_R} v_R {\mathrm d}x = 0 , \\ & \int_{\Omega } \Bigl( \nabla u _I \cdot \nabla v_I \;- k^2 u_I v_I \Bigr) {\mathrm d}x - { \int_{\Omega } f_I v_I {\mathrm d}x } = 0 . \\ \end{align}\tag{30}\] Since \(C_c^\infty(\Omega)\subset {\mathscr H}\) is dense, testing 29 with \(v\in C_c^\infty(\Omega)\) yields the Euler-Lagrange equation in distributional form. Considering test functions \(v\in {\mathscr H}\cap C ^\infty _ c (\Omega)\) we conclude that stationary points of \(\mathcal{E}_P\) satisfy \(- \Delta u - k^2 u = f\) in \(\Omega,\) in the sense of distributions.

As discussed in Section 2, we cannot incorporate the absorbing boundary conditions into the energy when we restrict the computational domain. One way to view this is to split the domain as \(\Omega = \Omega_{\text{Restr}} \cup \Omega_2\). The corresponding energies can then be written as follows, assuming that \(f\) has compact support in \(\Omega_{\text{Restr}}\): \[\label{Energy95H95main95s2952} \begin{align} \mathcal{E}_P(v) = &\int_{ \Omega _{\text{Restr}} } \frac{1}{2} \Bigl( |\nabla v |^2 - k^2|v|^2 \Bigr) {\mathrm d}x - \Re \int_{\Omega _{\text{Restr}} } f \overline{v} {\mathrm d}x + \int_{\Omega _{2} } \frac{1}{2} \Bigl( |\nabla v |^2 - k^2|v|^2 \Bigr) {\mathrm d}x \\ =&: \mathcal{E}_{P, \Omega _{\text{Restr}} }(v) + \mathcal{E}_{P, \Omega _{2}}(v) \end{align}\tag{31}\] An energetically consistent approach requires representing the contribution of the complementary domain, \(\mathcal{E}_{P, \Omega _{2}}(v)\), through a surface integral incorporating appropriate boundary conditions. As discussed in Section 2, it remains unclear whether impedance conditions can be used to represent this energy contribution. To preserve the variational structure associated with the boundary conditions in (1.1), we therefore leverage complex Lagrangians.

To this end, we define \[\label{Energy95H95complex951} \begin{align} \mathcal{E}_{P, {\mathbb{C}}} (v) = \mathcal{E}_{P, {\mathbb{C}}} (v, \overline{v}) = &\int_{ \Omega } \frac{1}{2} \Bigl( |\nabla v |^2 - k^2|v|^2 \Bigr) {\mathrm d}x - \Re \int_{\Omega } f \overline{v} {\mathrm d}x - \frac{ik}{2}\int_{\partial \Omega}|v|^2 {\mathrm d}x \\ = &\int_{ \Omega } \frac{1}{2} \Bigl( \nabla v\cdot \nabla \overline{v}- k^2v \overline{v}\Bigr) {\mathrm d}x - \Re \int_{\Omega } f \overline{v} {\mathrm d}x - \frac{ik}{2}\int_{\partial \Omega}v \overline{v}{\mathrm d}x \end{align}\tag{32}\]

The calculations leading to 22 show the following

Regularised Lagrangians for Helmholtz↩︎

The Lagrangians derived above provide a variational characterisation of the Helmholtz equation from first principles. However, they are indefinite and lack coercivity. To address this, we construct regularised functionals that penalise the residual of the Helmholtz operator in a least-squares sense. In doing so, we aim to gain a deeper understanding of the analytical properties of the reduced wave equation and, in turn, develop tools for the design of numerical approximations with potentially improved properties.

Specifically, the introduction of the following real Lagrangian will be instrumental \[\begin{align} \mathcal{F}_{\gamma} (v) &= \int_{\Omega } \frac{1}{2} \Bigl( |\nabla v |^2 - k^2|v|^2 \Bigr) {\mathrm d}x - {\Re \int_{\Omega } f \overline{v} {\mathrm d}x} \\ &\quad + \gamma_1 \int_{\Omega} \bigl| {\mathscr{L}}v - f \bigr|^2 {\mathrm d}x + \gamma_2 \int_{\partial \Omega} \Bigl| \partial_{\mathbf{n}} v - i k v \Bigr|^2 {\mathrm d}S. \end{align}\] Following the above reasoning, we introduce the corresponding complex Lagrangian to incorporate the influence of the impedance boundary conditions into the principal part \[\label{Complex-Energy95H95gamma95wbc} \begin{align} \mathcal{F}_{\gamma,\, {\mathbb{C}}\, } (v) &= \int_{\Omega } \frac{1}{2} \Bigl( |\nabla v |^2 - k^2|v|^2 \Bigr) {\mathrm d}x - {\Re \int_{\Omega } f \overline{v} {\mathrm d}x} \\ &\quad + \gamma_1 \int_{\Omega} \bigl| {\mathscr{L}}v - f \bigr|^2 {\mathrm d}x + \gamma_2 \int_{\partial \Omega} \Bigl| \partial_{\mathbf{n}} v - i k v \Bigr|^2 {\mathrm d}S - \frac{ik}{2}\int_{\partial \Omega}|v|^2 {\mathrm d}x . \end{align}\tag{34}\] To emphasise the role of \(v\) and \(\overline{v}\) as independent variables when differentiating complex Lagrangians, and using the obvious notation, we denote \[\label{Complex-Energy95v95vbar} \begin{align} \mathcal{F}_{\gamma,\, {\mathbb{C}}\, } (v, \overline{v}) = \, &\mathcal{F}_{\gamma,\, {\mathbb{C}}\, } (v)\, . \end{align}\tag{35}\] The next section is devoted to the proof of our key result, stating that the homogeneous part of \(F_\gamma(u),\) \(F_\gamma(u):=\mathcal{F}_\gamma(u)\big|_{f=0}\) is coercive in \({\mathscr V}.\) In particular, we shall assume the following geometric assumption on \(\Omega\), compare to Moiola_Spence2014?; see Remark [Rem3.6].

The following coercivity result is the key analytical input. Its proof is given in Section 4.

As a consequence of Theorem 1, we have the following result.

Statement (i) is precisely Theorem 1. Assuming this coercivity estimate, statement (ii) follows from the complex version of the Lax–Milgram theorem Dautray_Lions_v2_1988?. To this end, consider the sesquilinear form \(\mathcal{A}_{\gamma,\, {\mathbb{C}}\, }\) defined by \[\label{Aform95reg95wbc} \begin{align} \mathcal{A}_{\gamma,\, {\mathbb{C}}\, }(u,v) &:= \int_\Omega \nabla u\cdot \nabla \overline{v} {\mathrm d}x - k^2 \int_\Omega u \overline{v} {\mathrm d}x \\ &\quad + 2\gamma_1 \int_{\Omega} ({\mathscr{L}}u) \overline{{\mathscr{L}}v} {\mathrm d}x + 2\gamma_2 \int_{\partial\Omega} \Bigl(\partial_{\mathbf{n}}u - i k u \Bigr) \overline{\Bigl(\partial_{\mathbf{n}}v - i k v \Bigr)} {\mathrm d}S \\ &\quad - {ik} \int_{\partial \Omega}u \, \overline{v}{\mathrm d}x . \end{align}\tag{36}\] This form is linear in its first argument and antilinear in its second, and it arises as the first variation of \(\mathcal{F}_{\gamma,\, {\mathbb{C}}\, }\) with respect to \(\overline{v}.\) We conclude that \[\label{Energy95H95gamma95Re95295wbc} \begin{align} \langle D_{\overline{v}\, } \mathcal{F}_{\gamma,\, {\mathbb{C}}\, }(u), v \rangle = \Bigl\{ \mathcal{A}_{\gamma,\, {\mathbb{C}}\, }(u,v) - \int_{\Omega} f \overline{\bigl(v + 2\gamma_1 {\mathscr{L}}v\bigr)} {\mathrm d}x \Bigr\}. \end{align}\tag{37}\] Thus, the stationary points of \(\mathcal{F}_\gamma\) satisfy \[\label{Energy95H95gamma95Re95395wbc} \mathcal{A}_{\gamma,\, {\mathbb{C}}\, } (u,v) = \int_{\Omega} f \overline{\bigl(v + 2\gamma_1 {\mathscr{L}}v\bigr)} {\mathrm d}x , \qquad \text{for all } v \in {\mathscr V}.\tag{38}\] Theorem 1 implies \[\Re \; \frac{1}{2} \mathcal{A}_{\gamma,\, {\mathbb{C}}\, }(v,v) = \Re \; \mathcal{F}_{\gamma,\, {\mathbb{C}}\, } (v)\Big|_{f=0} \ge \tilde{c} \|v\|_{{\mathscr V}}^{2}, \qquad \text{for all } v\in {\mathscr V}.\] Thus \(\mathcal{A}_{\gamma,\, {\mathbb{C}}\, }\) is bounded on \({\mathscr V}\times{\mathscr V}\) and \(\mathcal{A}_{\gamma,\, {\mathbb{C}}\, }\) is \({\mathscr V}\)-elliptic, while the right-hand side \(v\mapsto \int_\Omega f \overline{(v+2\gamma {\mathscr{L}}v)}\) defines a bounded antilinear functional on \({\mathscr V}\). Therefore, by the Lax-Milgram theorem Dautray_Lions_v2_1988? the variational problem 38 has a unique solution. To complete the proof, consider \(u\) being the unique weak solution of 1 . Then \(u\) solves \[{\mathcal{B}}(u,v) = \int_{\Omega } f \overline{v} {\mathrm d}x, \quad \text{for all } v\in {\mathscr H}\, .\] Furthermore, since \(\Omega\) is a Lipschitz domain, \(u \in {\mathscr V},\) see Remark [Rem3.4] . In addition, \({\mathscr{L}}u =f\) in \(L^2 (\Omega)\) and the impedance boundary conditions are satisfied in the \(L^2 ( \partial \Omega)\) sense. Since \({\mathscr V}\subset {\mathscr H}\) we conclude that \(u\) satisfies 38 and it is exactly the unique stationary point of \(\mathcal{F}_{\gamma,\, {\mathbb{C}}\, }\) with respect to \(\overline{v}\, .\)

Related minimisation problems↩︎

Our analysis may lead to alternative approaches for approximating 1 . In particular, it provides a mathematical foundation for iterative methods that incorporate boundary conditions either as Lagrange multipliers or through fixed-point iterations. To fix ideas, consider the following problem: assume that an approximation of \(u,\) denoted by \(\tilde{w}\), is given, and we seek to update it by solving the following problem: Find \(w \in {\mathscr V}\) such that \[\label{fixed95p951} \begin{align} \int_\Omega \nabla w\cdot \nabla \overline{v} {\mathrm d}x &- k^2 \int_\Omega w \overline{v} {\mathrm d}x + 2\gamma_1 \int_{\Omega} ({\mathscr{L}}w) \overline{{\mathscr{L}}v} {\mathrm d}x + 2\gamma_2 \int_{\partial\Omega} \Bigl(\partial_{\mathbf{n}}w - i k w \Bigr) \overline{\Bigl(\partial_{\mathbf{n}}v - i k v \Bigr)} {\mathrm d}S \\ &\quad = {ik} \int_{\partial \Omega}\widetilde{w} \, \overline{v}\, {\mathrm d}S + \int_{\Omega} f \overline{\bigl(v + 2\gamma_1 {\mathscr{L}}v\bigr)} {\mathrm d}x , \qquad \text{for all } v \in {\mathscr V}. \end{align}\tag{39}\] Consider this problem with all data fixed except for \(\widetilde{w}.\) We shall show that it admits a unique solution, and denote by \(T\) the solution operator that maps \(\widetilde{w} \mapsto w :\) \[T : {\mathscr V}\to {\mathscr V}\, , \qquad w = T ( \widetilde{w} )\, .\] We consider the sesquilinear form \(\mathcal{A}_{\text{wbc}}\) defined by \[\begin{align} \mathcal{A}_{\text{wbc}}(u,v) &:= \int_\Omega \nabla u\cdot \nabla \overline{v} {\mathrm d}x - k^2 \int_\Omega u \overline{v} {\mathrm d}x \\ &\quad + 2\gamma_1 \int_{\Omega} ({\mathscr{L}}u) \overline{{\mathscr{L}}v} {\mathrm d}x + 2\gamma_2 \int_{\partial\Omega} \Bigl(\partial_{\mathbf{n}}u - i k u \Bigr) \overline{\Bigl(\partial_{\mathbf{n}}v - i k v \Bigr)} {\mathrm d}S . \end{align}\] Then \(w\) solves \[\label{fixed95p952} \mathcal{A}_{\text{wbc}} (w,v) = {ik} \int_{\partial \Omega}\widetilde{w} \, \overline{v}\, {\mathrm d}S + \int_{\Omega} f \overline{\bigl(v + 2\gamma_1 {\mathscr{L}}v\bigr)} {\mathrm d}x \qquad \text{for all } v \in {\mathscr V}.\tag{40}\] Notice that \(\mathcal{A}_{\text{wbc}}\) is Hermitian (conjugate-symmetric), i.e. \(\mathcal{A}_{\text{wbc}} (u,v)=\overline{\mathcal{A}_{\text{wbc}} (v,u)}\). The next theorem shows that \(w\) can be obtained as a minimiser of the real Lagrangian \[\label{Energy95H95gamma95wbc} \begin{align} \widetilde{\mathcal{F}}_{\gamma} (v) &= \int_{\Omega } \frac{1}{2} \Bigl( |\nabla v |^2 - k^2|v|^2 \Bigr) {\mathrm d}x - {\Re \;{ik} \int_{\partial \Omega}\widetilde{w} \, \overline{v}\, {\mathrm d}S} - {\Re \int_{\Omega } f \overline{v} {\mathrm d}x}\\ &\quad + \gamma_1 \int_{\Omega} \bigl| {\mathscr{L}}v - f \bigr|^2 {\mathrm d}x + \gamma_2 \int_{\partial \Omega} \Bigl| \partial_{\mathbf{n}} v - i k v \Bigr|^2 {\mathrm d}S \, . \end{align}\tag{41}\]

As before (i) follows by Theorem 1. Statements (ii) and (iii) then follow from the complex Lax-Milgram theorem. In fact, the sesquilinear form \(\mathcal{A}_{\text{wbc}}\) arises as the first variation of the real Lagrangian \(\widetilde{\mathcal{F}}_{\gamma}\). We conclude that \[\label{grzhliop} \begin{align} \langle D\widetilde{\mathcal{F}}_{\gamma}(u), v \rangle = \Re \Bigl\{ \mathcal{A}_{\text{wbc}}(u,v) - \;{ik} \int_{\partial \Omega}\widetilde{w} \, \overline{v}\, {\mathrm d}x - \int_{\Omega} f \overline{\bigl(v + 2\gamma_1 {\mathscr{L}}v\bigr)} {\mathrm d}x \Bigr\}. \end{align}\tag{42}\] Thus, the stationary points of \(\widetilde{\mathcal{F}}_{\gamma}\) are solutions of 40 . Theorem 1 implies \[\label{Energy95H95gamma95Re95495wbc} \Re \frac{1}{2} \mathcal{A}_{\text{wbc}}(v,v) \ge \tilde{c} \|v\|_{{\mathscr V}}^{2}, \qquad \text{for all } v\in {\mathscr V}.\tag{43}\] In addition the right-hand side of 40 defines a bounded antilinear functional on \({\mathscr V}\). Therefore, by Dautray_Lions_v2_1988?, the variational problem 40 admits a unique solution, and statement (ii) follows. The fact that \(T\) is a contraction in \(L^2(\partial \Omega)\) follows from the observation that the boundary term involving \(\widetilde{w}\) in 40 is independent of \(\gamma_1\) and \(\gamma_2\), together with the bound ?? , provided that \(\beta\) is chosen such that in \(\Big(2\beta - \alpha L - \frac{\alpha L^2}{\varepsilon_2}\Big)\) is sufficiently large, which in turn implies that \(\gamma_1\) must be taken sufficiently large. The remainder of statement (iii) is then immediate.

A convergent fixed point iteration↩︎

Let \(\rho >0\) and consider the following iterative scheme: \[\label{fixed95point95it95tilde} \widetilde{w} ^{n+1} = (1-\rho)\widetilde{w} ^{n} + \rho\, w^{n+1} = (1-\rho)\widetilde{w} ^{n} + \rho\, T (\widetilde{w} ^{n}) =: G (\widetilde{w} ^{n} )\, .\tag{44}\] Utilising (iii) of the above theorem we have \[\|G( \widetilde{w} _ 1) - G(\widetilde{w} _2) \| _{L^2 (\partial \Omega)} \leq \, \Big [(1-\rho) + \eta \rho \, \Big ]\, \| \widetilde{w} _ 1 - \widetilde{w} _2\| _{L^2 (\partial \Omega)} \, .\] Since \(0< \eta < 1\) it follows that \(G\) is a contraction for any \(\rho ,\) \(0< \rho < 1\, .\) Hence the sequence converges to the fixed point \(u\) in \(L^2 (\partial \Omega)\, .\) By the stability bound, \[\| u- T(\widetilde{w} ^{n}) \| _{{\mathscr V}} \leq C \, k \, \| u - \widetilde{w} ^{n} \| _{L^2 (\partial \Omega)} \, ,\] we conclude that \(T(\widetilde{w} ^{n}) \to u\) in \({\mathscr V}\, .\)

4 Proofs for strongly coercive variational principles↩︎

Design strategy↩︎

Our construction of the enhanced Lagrangians is intrinsically connected to the analysis establishing their coercivity. The key idea is to observe that the energy \(\mathcal{E}_P (u)\) contains the indefinite term \[\label{energy95wd} \frac{1}{2} \int_\Omega \!\left({\left|\nabla u\right|^2 - k^2 \left|u\right|^2}\right) \, {\mathrm d}x \, ,\tag{48}\] and at the same time a standard Rellich identity, see 51 , can be written as \[\label{eq:IntRellich95wd} \int_\Omega \!\left({(\nu -2) \left|\nabla u\right|^2 - \nu k^2 \left|u\right|^2}\right) \, {\mathrm d}x \\ + \int_\Omega 2 \Re\!\left({\!\left({\boldsymbol{x} \cdot \nabla u}\right) \overline{{\mathscr{L}}u}}\right) \, {\mathrm d}x + B_R(u)=0,\tag{49}\] where \(B_R(u)\) accounts for the boundary terms in 51 . Then, an appropriate linear combination of the principal part of the Lagrangian for \(f=0\) and the identity 49 has two key advantages: the energy changes only by a constant factor, since the contribution of 49 is zero, and the coefficients of \(\left|\nabla u\right|^2\) and \(k^2 \left|u\right|^2\) can both be made constants of the same sign. The least-squares regularisation term in the Lagrangian functional makes it possible to control the second term in 49 . Moreover, the boundary term in 49 can also be controlled, particularly through the use of another lower-order Morawetz identity—a fact that is not immediately apparent. It is important to emphasise that, within this framework, the structure of the regularised Largangian is preserved, and coercivity is ensured for all \(\gamma \geq \gamma_0\), where the constant \(\gamma_0\) can be determined explicitly.

Next, to simplify the presentation of our approach and the proof of Theorem 4.2, we first establish the coercivity bound under the assumption that the impedance condition is satisfied exactly; see Theorem 4.1. This is done for all elements of the space \({\mathscr V}_{BC}\subset {\mathscr V}\).

Rellich and Morawetz Identities↩︎

Let \({\mathscr{L}}u = -\Delta u - k^2 u\) be the Helmholtz operator applied to a sufficiently smooth function \(u \colon \mathbb{R}^\nu \to \mathbb{C}\). Then the following identity holds: \[\begin{align} \label{rellich95id} &\!\left({\boldsymbol{x} \cdot \overline{\nabla u}}\right) {\mathscr{L}}u + \!\left({\boldsymbol{x} \cdot \nabla u}\right) \overline{{\mathscr{L}}u} \nonumber\\ &= -\text{\rm div}\!\left({ (\boldsymbol{x} \cdot \overline{\nabla u}) \nabla u + (\boldsymbol{x} \cdot \nabla u) \overline{\nabla u} - \boldsymbol{x} \left|\nabla u\right|^2 + k^2 \boldsymbol{x} \left|u\right|^2 }\right) - (\nu -2) \left|\nabla u\right|^2 + \nu k^2 \left|u\right|^2 , \end{align}\tag{50}\] where \(\boldsymbol{x} \cdot \nabla u = \sum_{{ j=1}}^{\nu } x_j \frac{\partial u}{\partial x_j}\) is the directional derivative of \(u\) along \(\boldsymbol{x}\).

Although this identity is well known CummingsFeng2006?, Spence_Chandler-Wilde_integral_2011?, Moiola_Spence2014? as the simplest form of a Morawetz (or Rellich) identity, we briefly outline the main steps of its derivation to highlight the role of the Morawetz multiplier.

To begin, let \(\Phi\) denote the left-hand side of identity 50 . Using the definition of the Helmholtz operator \({\mathscr{L}}\), we have \[\Phi = \!\left({\boldsymbol{x} \cdot \overline{\nabla u}}\right) \!\left({-\Delta u - k^2 u}\right) + \!\left({\boldsymbol{x} \cdot \nabla u}\right) \!\left({-\Delta \overline{u} - k^2 \overline{u}}\right).\] Expanding each term gives \[\Phi = -\!\left({\boldsymbol{x} \cdot \overline{\nabla u}}\right) \Delta u - k^2 \!\left({\boldsymbol{x} \cdot \overline{\nabla u}}\right) u - \!\left({\boldsymbol{x} \cdot \nabla u}\right) \Delta \overline{u} - k^2 \!\left({\boldsymbol{x} \cdot \nabla u}\right) \overline{u}.\] The terms involving the Laplacian expand as follows: \[-\!\left({\boldsymbol{x} \cdot \overline{\nabla u}}\right) \Delta u - \!\left({\boldsymbol{x} \cdot \nabla u}\right) \Delta \overline{u} = -\text{\rm div}\!\left({\!\left({\boldsymbol{x} \cdot \overline{\nabla u}}\right) \nabla u + \!\left({\boldsymbol{x} \cdot \nabla u}\right) \overline{\nabla u} { - \boldsymbol{x} \left|\nabla u\right|^2}}\right) { - (\nu - 2)} \left|\nabla u\right|^2,\] where we have used that \[\text{\rm div}\!\left({\boldsymbol{x} \left|\nabla u\right|^2}\right) = { \nu }\, \left|\nabla u\right|^2 + \!\left({\boldsymbol{x} \cdot \nabla}\right)\!\left({\left|\nabla u\right|^2}\right).\] For the remaining two terms, we observe that they can be written in the form \[- k^2 \!\left({\!\left({\boldsymbol{x} \cdot \overline{\nabla u}}\right) u + \!\left({\boldsymbol{x} \cdot \nabla u}\right) \overline{u}}\right) = - k^2 \!\left({\boldsymbol{x} \cdot \nabla \!\left({\left|u\right|^2}\right)}\right).\] Rewriting this in divergence form, we obtain \[- k^2 \!\left({\boldsymbol{x} \cdot \nabla \!\left({\left|u\right|^2}\right)}\right) = -\text{\rm div}\!\left({k^2 \boldsymbol{x} \left|u\right|^2}\right) + \nu k^2 \left|u\right|^2.\] Combining all the contributions derived above, we recover identity 50 . Integration over \(\Omega\), together with the divergence theorem and a density argument (see, e.g., CummingsFeng2006?, Spence_Chandler-Wilde_integral_2011?, Moiola_Spence2014?) yields the integrated form of identity 50 , stated in the following proposition.

The following identity will be important for the analysis that follows.

Lagrangians and coercivity bounds↩︎

Let \(\Omega \subset \mathbb{R}^\nu\) be a bounded Lipschitz domain (for the Rellich identities we take the multiplier \(x\), which may be replaced by \(x-x_0\) to avoid assuming the origin lies in \(\Omega\)). Consider the Helmholtz equation with a source term: \[{\mathscr{L}}u = f, \qquad \text{where } {\mathscr{L}}u := -\Delta u - k^2 u,\] subject to the impedance boundary condition \[\label{bc95section4} \nabla u \cdot \boldsymbol{n} = i k u \quad \text{on } \partial \Omega.\tag{53}\] We shall first examine the coercivity of \[\label{Egamma} E_\gamma(u) := \frac{1}{2} \int_\Omega \!\left({\left|\nabla u\right|^2 - k^2 \left|u\right|^2}\right) {\mathrm d}x + \gamma \int_\Omega \left|{\mathscr{L}}u\right|^2 {\mathrm d}x,\tag{54}\] in the space \({\mathscr V}_{BC}\, .\)

To control the mixed volume term in ?? we use the following simple estimate.

Now consider the boundary energy term ?? , which we rewrite as \[\begin{align} \label{dimxEgamma95boundary} E_{\gamma,\partial\Omega}(\nu ,\alpha;u) =& \alpha \int_{\partial \Omega} \Big( \left|\nabla u\right|^2 (\boldsymbol{x}\cdot\boldsymbol{n}) - k^2 \left|u\right|^2 (\boldsymbol{x}\cdot\boldsymbol{n}) \Big) {\mathrm d}S - \alpha \int_{\partial \Omega} 2 \Re\!\left({(\boldsymbol{x}\cdot\overline{\nabla u}) (\nabla u\cdot\boldsymbol{n})}\right) {\mathrm d}S. \end{align}\tag{55}\] Next, utilising Assumption 1, we estimate \(E_{\gamma,\partial \Omega}\) as follows: \[\begin{align} \label{dimxEgamma95boundary2} E_{\gamma,\partial \Omega}(\nu ,\alpha;u) \ge& \alpha L_0 \int_{\partial \Omega} |\nabla u|^2 {\mathrm d}S - {\alpha} L k^2 \int_{\partial \Omega} |u|^2 {\mathrm d}S - \alpha \int_{\partial \Omega} 2 \Re \big((\boldsymbol{x}\cdot \overline{\nabla u}) (\nabla u\cdot \boldsymbol{n})\big) {\mathrm d}S . \end{align}\tag{56}\]

First consider the case when \(u\) exactly satisfies the boundary condition 53 . With the notation of Proposition [the:boundarymorawetz] we have \[\label{eq:M095identity95strong} - \int_\Omega 2 \Re \big(\overline{M_0 u} {\mathscr{L}}u\big) {\mathrm d}x + 2k^2\beta \int_{\partial\Omega} |u|^2 {\mathrm d}S = 0,\tag{57}\] where \(M_0 u := -ik\beta u\). Adding the null quantity 57 to the right-hand side of 56 and using 53 gives \[\begin{align} \label{dimxEgamma95boundary3} E_{\gamma,\partial \Omega}(\nu ,\alpha;u) \ge & \alpha L_0 \int_{\partial \Omega} |\nabla u|^2 {\mathrm d}S + \bigl(2\beta - L\bigr) k^2 \int_{\partial \Omega} |u|^2 {\mathrm d}S \nonumber\\ & - \alpha \int_{\partial \Omega} 2 \Re \big((\boldsymbol{x}\cdot \overline{\nabla u}) (ik u)\big) {\mathrm d}S - \int_\Omega 2 \Re \big(\overline{M_0 u} {\mathscr{L}}u\big) {\mathrm d}x . \end{align}\tag{58}\] Using Assumption 1 and Young’s inequality, \[\label{dimxEgamma95boundary95inder951} \begin{align} \alpha \int_{\partial \Omega} 2 \Re \big((\boldsymbol{x}\cdot \overline{\nabla u}) (ik u)\big) {\mathrm d}S &\le \alpha \int_{\partial \Omega} 2L |\nabla u| k|u| {\mathrm d}S \leq \frac{\alpha L_0}{2}\int_{\partial \Omega} |\nabla u|^2 {\mathrm d}S + \frac{2\alpha L^2}{L_0} k^2 \int_{\partial \Omega} |u|^2 {\mathrm d}S . \end{align}\tag{59}\] Furthermore, by Young’s inequality (for any \(\varepsilon_3>0\)), \[\label{dimxEgamma95boundary95inder952} \begin{align} \int_\Omega 2 \Re \big(\overline{M_0 u} {\mathscr{L}}u\big) {\mathrm d}x &\le \int_\Omega 2 \beta k |u| |{\mathscr{L}}u| {\mathrm d}x \le \varepsilon_3 k^2 \int_{\Omega} |u|^2 {\mathrm d}x + \frac{\beta^2}{\varepsilon_3}\int_{\Omega} |{\mathscr{L}}u|^2 {\mathrm d}x . \end{align}\tag{60}\]

Combining 58 60 with the bulk estimate and an appropriate choice of parameters \((\alpha,\beta,\varepsilon_1,\varepsilon_3)\) yields the boundary control needed for the global coercivity. This leads to the theorem stated next.

Proof of Theorem 3.1↩︎

We now consider the case where \(u\) satisfies the impedance condition 53 weakly and examine the coercivity of \(F_\gamma(u):=\mathcal{F}_\gamma(u)\big|_{f=0}= \, \Re \, \mathcal{F}_{\gamma,\, {\mathbb{C}}\, } (v) \big|_{f=0}\).

We have \[\label{dimxFgamma} \nu F_\gamma(u) = F_{\gamma_1}(u) + F_{\gamma_2}(u),\tag{61}\] where \(F_{\gamma_1}(u) := E_{\gamma_1,\Omega}(u)+E_{\gamma_1,\partial\Omega}(u)\) and \(F_{\gamma_2}(u) := \nu \gamma_2\int_{\partial\Omega}\Bigl|\partial_{\mathbf{n}}u-ik u\Bigr|^2 {\mathrm d}S.\) The bulk estimate for \(E_{\gamma_1,\Omega}\) is identical to that for \(E_{\gamma,\Omega}\) from the strongly imposed case with \(\gamma\) replaced by \(\gamma_1\).

For the boundary term \(E_{\gamma_1,\partial\Omega}\) we start from 56 (with \(\gamma=\gamma_1\)) and split the normal derivative: \[\begin{align} -\alpha\int_{\partial\Omega}2 \Re \big((\boldsymbol{x} \cdot \overline{\nabla u})(\nabla u \cdot \boldsymbol{n})\big) {\mathrm d}S &= -\alpha\int_{\partial\Omega}2 \Re \big((\boldsymbol{x} \cdot \overline{\nabla u}) (\partial_{\mathbf{n}}u-ik u)\big) {\mathrm d}S -\alpha\int_{\partial\Omega}2 \Re \big((\boldsymbol{x} \cdot \overline{\nabla u}) ik u\big) {\mathrm d}S. \end{align}\] By Cauchy-Schwarz/Young’s inequality, \(| \boldsymbol{x} \cdot \nabla u |\le L|\nabla u|\), so for any \(\varepsilon_3>0\), \[-\alpha\int_{\partial\Omega}2 \Re \big((\boldsymbol{x} \cdot \overline{\nabla u}) (\partial_{\mathbf{n}}u-ik u)\big) {\mathrm d}S \ge -\alpha L^2\varepsilon_3 \|\nabla u\|^2_{\partial\Omega} - \frac{\alpha}{\varepsilon_3} \|\partial_{\mathbf{n}}u-ik u\|^2_{\partial\Omega}.\] Likewise, for any \(\varepsilon_2>0\), \[-\alpha\int_{\partial\Omega}2 \Re \big((\boldsymbol{x} \cdot \overline{\nabla u}) ik u\big) {\mathrm d}S \ge -\alpha\varepsilon_2 \|\nabla u\|^2_{\partial\Omega} - \frac{\alpha L^2}{\varepsilon_2} k^2\|u\|^2_{\partial\Omega}.\] To handle the bulk coupling that appears in the strong case via \(M_0\), we estimate directly by Young’s inequality: for any \(\beta>0\) and any \(\varepsilon_4>0\), \[\int_\Omega 2 \Re \big(\overline{(-ik\beta u)} {\mathscr{L}}u\big) {\mathrm d}x \le \varepsilon_4 k^2\|u\|^2 + \frac{\beta^2}{\varepsilon_4} \|{\mathscr{L}}u\|^2,\] which yields the same volume contributions as in the strong case after redistribution into the bulk bound.

Collecting the above bounds with 56 gives \[\label{coerc95195wbc} \begin{align} E_{\gamma_1,\partial\Omega}(u) \ge & \alpha L_0 \|\nabla u\|^2_{\partial\Omega} + (2\beta-\alpha L) k^2\|u\|^2_{\partial\Omega} - \alpha L^2\varepsilon_3 \|\nabla u\|^2_{\partial\Omega} - \frac{\alpha}{\varepsilon_3} \|\partial_{\mathbf{n}}u-ik u\|^2_{\partial\Omega}\\ &- \alpha\varepsilon_2 \|\nabla u\|^2_{\partial\Omega} - \frac{\alpha L^2}{\varepsilon_2} k^2\|u\|^2_{\partial\Omega} - \frac{(\alpha-\tfrac12) \nu }{2} k^2\|u\|^2 - \frac{2\beta^2}{(\alpha-\tfrac12) \nu } \|{\mathscr{L}}u\|^2. \end{align}\tag{62}\] Adding \(F_{\gamma_2}(u)\) from 61 yields the net boundary-residual coefficient \({\big(\nu \gamma_2 - \tfrac{\alpha}{\varepsilon_3}\big) \|\partial_{\mathbf{n}}u-ik u\|^2_{\partial\Omega}.}\) Combining the bulk estimate (with \(\gamma = \gamma_1\)) and 62 we obtain the result of Theorem 1.

One admissible parameter choice is obtained by taking \(\alpha\in(\tfrac12,\frac{\nu }{2(\nu -2)})\) (or any \(\alpha>\tfrac12\) if \(\nu \in\{1,2\}\)), then fixing small \(\varepsilon_1,\varepsilon_2,\varepsilon_3>0\) so that the boundary coefficients are positive, choosing \(\beta\) to make \(2\beta-\alpha L - \frac{\alpha L^2}{\varepsilon_2}>0\), and finally taking \(\gamma_1,\gamma_2\) above the indicated thresholds.

5 Variational discretisation methods↩︎

This section describes the discrete formulations used in our numerical examples. We first consider a conforming finite element discretisation of the augmented variational principle introduced in Section 3. The neural-network formulation is discussed in the following subsection. The purpose of the experiments is verification and illustration, a systematic computational comparison is left for future work.

5.1 Finite element approximation↩︎

Discrete formulation↩︎

Let \(V_h\subset {\mathscr V}\) be a finite-dimensional conforming space. Since the augmented formulation contains the bulk and the impedance residual, all terms are well defined if, for example, \(V_h\subset H^2(\Omega)\). In the computations below we use \(H^2\)-conforming Argyris elements, although other \(C^1\) finite element spaces could also be used.

The finite element approximation is obtained by restricting the stationary formulation associated with \(\mathcal{F}_{\gamma,{\mathbb{C}}}\) to \(V_h\): find \(u_h\in V_h\) such that \[\label{eq:fem95augmented95form} \mathcal{A}_{\gamma,{\mathbb{C}}}(u_h,v_h) = \int_{\Omega} f \overline{\!\left({v_h+2\gamma_1{\mathscr{L}}v_h}\right)} {\mathrm d}x, \qquad \text{for all } v_h\in V_h,\tag{63}\] where \(\mathcal{A}_{\gamma,{\mathbb{C}}}\) is the sesquilinear form defined in 36 .

Indeed, coercivity of the real part of \(\mathcal{A}_{\gamma,{\mathbb{C}}}\) on \({\mathscr V}\) implies coercivity on the subspace \(V_h\), while boundedness follows directly from the definition of \({\mathscr V}\). The result is therefore an immediate finite-dimensional consequence of the complex Lax–Milgram theorem.

The same argument gives the usual quasi-optimality estimate in the augmented norm.

This estimate is a stability and best-approximation result in the norm of \({\mathscr V}\). It should not be interpreted as a pollution-free estimate in \(L^2(\Omega)\) or \(H^1(\Omega)\) norms. For this reason, the verification tests below also report relative \(L^2\) errors, relative wavenumber-weighted \(H^1\) errors and residual diagnostics.

In the finite element implementation, the complex-valued problem is written as a mixed real-valued problem for \((u_R,u_I)\). All integrals are assembled in Firedrake using quadrature degree \(12\). Unless otherwise stated, the resulting sparse systems are solved by a direct LU factorisation using MUMPS through PETSc.

For the numerical tests we use a mesh-scaled boundary-residual weight \[\gamma_{2,h}=\frac{\tau}{h},\] where \(h\) is the mesh parameter and \(\tau>0\) is fixed. Thus the boundary contribution is evaluated with \(\gamma_2=\gamma_{2,h}\). This is a discrete stabilisation choice used to balance the boundary residual with the interior residual on the finite element space. Unless otherwise stated, we take \(\gamma_1=2\) and \(\tau=50\).

Manufactured-solution verification↩︎

We first consider a manufactured-solution test on the unit square \(\Omega=(0,1)^2\). The exact solution is \[u_{\mathrm{ex}}(x,y) = \exp\!\left({ ik\!\left({(x-\tfrac 12)^2+(y-\tfrac 12)^2}\right) }\right).\] On each side of the square, the outward normal derivative satisfies \[\partial_{\mathbf{n}}u_{\mathrm{ex}}-iku_{\mathrm{ex}}=0, \qquad \text{on } \partial\Omega.\] The corresponding right-hand side is defined by \[f_{\mathrm{ex}}=-\Delta u_{\mathrm{ex}}-k^2u_{\mathrm{ex}},\] namely \[f_{\mathrm{ex}}(x,y) = \!\left({ -4ik +4k^2\!\left({(x-\tfrac 12)^2+(y-\tfrac 12)^2}\right) -k^2 }\right) u_{\mathrm{ex}}(x,y).\]

We report the relative \(L^2(\Omega)\) error, the relative \(H^1_k(\Omega)\) error, the relative bulk residual \[\frac{ \|-\Delta u_h-k^2u_h-f_{\mathrm{ex}}\|_{L^2(\Omega)} }{ \|f_{\mathrm{ex}}\|_{L^2(\Omega)} },\] and the relative impedance residual \[\frac{ \|\partial_{\mathbf{n}}u_h-iku_h\|_{L^2(\partial\Omega)} }{ \|\partial_{\mathbf{n}}u_h\|_{L^2(\partial\Omega)} +k\|u_h\|_{L^2(\partial\Omega)} }.\]

Table 1 shows the refinement behaviour at fixed wavenumber \(k=25\). The relative \(L^2\) and \(H^1_k\) errors decrease under mesh refinement, and the bulk and boundary residuals decrease consistently. Table 2 reports a fixed-mesh wavenumber sweep with \(h=2^{-6}\). As expected, the errors increase as the wavenumber grows on a fixed mesh.

Localised-source illustration↩︎

We next consider a qualitative test with a smooth source concentrated near the centre of the domain: \[f(x,y) = \exp\!\left({ -\frac{(x-0.5)^2+(y-0.5)^2}{\varepsilon} }\right), \qquad \varepsilon=10^{-4}.\] The same Argyris discretisation and parameters \(\gamma_1=2\) and \(\tau=50\) are used. The computations are performed on \(\Omega=(0,1)^2\) using a uniform triangular mesh with \(h\approx 2^{-7}\). Figure 1 shows the real part of the computed solution for increasing wavenumber.

Nonconvex scattering benchmark↩︎

Finally, we include a nonconvex scattering benchmark. This example lies outside the star-shaped setting used in the coercivity proof, and is included as an illustrative stress test of the discrete formulation.

The computational domain is the unit square with a U-shaped sound-soft obstacle removed: \[\Omega=(0,1)^2\setminus\overline{D_U},\] where \[D_U = \!\left({[0.38,0.45]\times[0.30,0.75]}\right) \cup \!\left({[0.55,0.62]\times[0.30,0.75]}\right) \cup \!\left({[0.38,0.62]\times[0.30,0.38]}\right).\] The geometry contains re-entrant corners and a cavity opening.

We prescribe the incident plane wave \[u_{\mathrm{inc}}(x,y) = \exp\!\left({ik d\cdot(x,y)}\right), \qquad d=(0,-1),\] so that the wave is directed into the opening of the obstacle. The total field \(u\) satisfies \[-\Delta u-k^2u=0 \qquad \text{in } \Omega.\] On the obstacle boundary \(\Gamma_{\mathrm{obs}}:=\partial D_U\) we impose the sound-soft condition \[u=0 \qquad \text{on } \Gamma_{\mathrm{obs}}.\] On the exterior boundary \(\Gamma_{\mathrm{out}}:=\partial(0,1)^2\) we impose the first-order absorbing condition on the scattered field. Equivalently, the total field satisfies \[\partial_{\mathbf{n}}u-iku = \partial_{\mathbf{n}}u_{\mathrm{inc}}-iku_{\mathrm{inc}} \qquad \text{on } \Gamma_{\mathrm{out}}.\] In the implementation, the outer boundary residual is shifted by the known datum \[g_{\mathrm{out}} := \partial_{\mathbf{n}}u_{\mathrm{inc}}-iku_{\mathrm{inc}},\] while the sound-soft condition on \(\Gamma_{\mathrm{obs}}\) is imposed weakly through a mesh-scaled penalty term.

Figure 2 shows the augmented finite element solution for \(k=40\) on a mesh with \(h\approx0.045\). The incident wave enters the cavity from above, producing interference around the re-entrant corners and enhanced field intensity inside the U-shaped region.

5.2 Neural network approximation↩︎

Deep fixed point flux formulation↩︎

We also consider a neural network discretisation of the augmented variational formulation. The main difference from the finite element method is that the impedance contribution is represented through a boundary flux variable rather than assembled directly into a conforming trial space.

We introduce a complex boundary variable \(\lambda\) and write the weak Helmholtz equation as \[\label{eq:nn95flux95weak95form} \int_\Omega \nabla u\cdot\nabla\overline{v} {\mathrm d}x -k^2\int_\Omega u\overline{v} {\mathrm d}x -\int_{\partial\Omega}\lambda\overline{v} {\mathrm d}S = \int_\Omega f\overline{v} {\mathrm d}x .\tag{64}\] The impedance condition is encoded by the boundary relation \(\lambda=iku \qquad \text{on } \partial\Omega.\) At a fixed point, we recover the impedance term appearing in the variational formulation. For a fixed boundary flux \(\lambda\), the primal network is obtained by approximately minimising \[\label{eq:nn95deep95uzawa95functional} \begin{align} \mathcal{J}_{\gamma,\lambda}(v) &= \frac{1}{2}\int_\Omega\!\left({|\nabla v|^2-k^2|v|^2}\right){\mathrm d}x -\Re\int_\Omega f\overline{v}{\mathrm d}x +\gamma_1\int_\Omega|{\mathscr{L}}v-f|^2{\mathrm d}x \\ &\quad +\gamma_2\int_{\partial\Omega} |\partial_{\mathbf{n}}v-ikv|^2{\mathrm d}S -\Re\int_{\partial\Omega}\lambda\overline{v}{\mathrm d}S . \end{align}\tag{65}\] The last term supplies the weak flux contribution. Since this term is linear in \(v\), it does not affect the coercivity of the quadratic part, see Theorem 2.

The complex-valued approximation is represented as \[u_\theta=u_{\theta,R}+iu_{\theta,I},\] where \(u_{\theta,R}\) and \(u_{\theta,I}\) are the two real-valued outputs of a SIREN network. Derivatives in \({\mathscr{L}}u_\theta\) and on \(\partial\Omega\) are computed by automatic differentiation. On a fixed boundary collocation set \(\{y_m\}_{m=1}^{M_\partial}\subset\partial\Omega\), the flux is updated by the damped fixed point step, compare to Theorem 2 and 44 , \[\label{eq:nn95uzawa95update} \lambda^{m+1}(y) = (1-\rho)\lambda^m(y) +\rho\,ik\,u_{\theta^{m+1}}(y), \qquad y\in\partial\Omega,\tag{66}\] with \(0<\rho\le 1\). Equivalently, \[\lambda_R^{m+1} = (1-\rho)\lambda_R^m-\rho k u_I^{m+1}, \qquad \lambda_I^{m+1} = (1-\rho)\lambda_I^m+\rho k u_R^{m+1}.\]

Implementation details↩︎

All integrals in 65 are approximated by Monte Carlo quadrature using independently sampled interior and boundary points, with the appropriate domain and boundary measure factors. Unless otherwise stated, we use a SIREN network with two outputs, width \(96\), depth \(4\), Adam optimisation and random seed \(67\).

For the two-dimensional manufactured-solution tests we use \(4096\) interior points and \(2048\) boundary points per primal update. The Deep Fixed Point runs use \(200\) outer iterations, \(20\) primal Adam steps per outer iteration and a preliminary residual least-squares warm start of \(500\) Adam steps all in single precision. The relaxation parameter is fixed at \(\rho=0.2\).

As a baseline, we train a residual PINN with the same architecture, optimiser, sampling budget and approximately the same number of Adam updates. This baseline minimises the bulk residual and the impedance boundary residual directly, without the physical Helmholtz energy and without the flux update 66 .

Manufactured-solution comparison↩︎

We use the same two-dimensional manufactured solution as in the finite element verification test, \[u_{\mathrm{ex}}(x,y) = \exp\!\left({ ik\!\left({(x-\tfrac 12)^2+(y-\tfrac 12)^2}\right) }\right),\] with \(f_{\mathrm{ex}}=-\Delta u_{\mathrm{ex}}-k^2u_{\mathrm{ex}}\). The errors are evaluated on validation points not used in the corresponding training step. We report the relative \(L^2(\Omega)\) error, the relative \(H^1_k(\Omega)\) error, the relative bulk residual and the relative boundary residual when available.

Table 3 shows that the augmented formulation gives substantially smaller validation errors across this sweep. The comparison is merely a verification test for the proposed formulation not a complete performance study.

Three-dimensional obstacle scattering↩︎

We finally include a three-dimensional scattering example to illustrate the neural formulation on a non-star-shaped geometry. The computational domain is the unit cube with one or more vertical cylindrical sound-soft obstacles removed, \(\Omega=(0,1)^3\setminus\overline{D} .\) In the reported configuration, \(D\) is a staggered pair of cylinders, \(D = D_1 \cup D_2 ,\) where \[D_j = \{(x,y,z)\in(0,1)^3: (x-c_{j,1})^2+(y-c_{j,2})^2<R_j^2\}, \qquad j=1,2 .\] We prescribe the incident plane wave \[u_{\mathrm{inc}}(\boldsymbol{x}) = \exp\!\left({ik d\cdot \boldsymbol{x}}\right), \qquad d=\frac{(0.55,-1,0)}{\sqrt{0.55^2+1}} ,\] so that the incoming wave is oblique in the \(xy\)-plane. The total field \(u\) satisfies \(-\Delta u-k^2u=0\) in \(\Omega .\) On the cylindrical obstacle boundary \(\Gamma_{\mathrm{obs}}:=\partial D\) we impose the sound-soft condition \(u=0\) on \(\Gamma_{\mathrm{obs}}.\) This condition is imposed strongly in the neural ansatz.

The scattered field is defined by \(u_s=u-u_{\mathrm{inc}}.\) On the four vertical faces of the cube we impose the first-order absorbing condition on the scattered field, \(\partial_{\boldsymbol{n}}u_s-iku_s=0 .\) On the top and bottom faces, \(z=0\) and \(z=1\), we impose the corresponding axial consistency condition \(\partial_z u = \partial_z u_{\mathrm{inc}},\) which reduces to a homogeneous Neumann condition for the incident direction as \(d_3=0\).

The reported run uses the same neural residual least-squares structure as in the two-dimensional circular benchmark. Figure 3 shows two cutaway visualisations of the computed field. The plots display the three-dimensional extruded scattering pattern generated by the oblique incident wave interacting with the cylindrical sound-soft obstacles.

---

1. \(^{1}\) DMAM, University of Crete, Greece
   \(^{2 }\) IACM-FORTH, Greece
   \(^{3 }\) MPS, University of Sussex, United Kingdom
   \(^{4 }\) DNA, University of West Attica, Greece
   \(^{5 }\) DMS and Institute of Mathematical Innovation, University of Bath, United Kingdom↩︎