Defining a goal

The goal of this incipient research is to be able to compute the eigenvalues (specifically, the imaginary part) of the Helmholtz equation \[ \Delta u +k^2 n^2 u = \beta^2 u. \] for \( k \) a given wavenumber. Notice that, for the exact solution and any domain \(\omega \subset \Omega \), the following holds \[ \Im(\beta^2) = \frac{1}{2\imath \lVert u\rVert_\omega^2} \int_{\partial \omega} (\partial_n \overline{u})u - (\partial_n u)\overline{u} \, \mathrm{d}s, \] where this expression is obtained by means of testing the previous equation against the conjugate of \( u \), conjugating the equation and testing it against the solution, and taking diferences.

Thus, we define the functional \[ j_\omega(u) := \frac{1}{2\imath \lVert u\rVert_\omega^2} \int_{\partial \omega} (\partial_n \overline{u})u - (\partial_n u)\overline{u} \, \mathrm{d}s, \] for all subdomains \( \omega \subset \Omega \). In particular, we set \( J(u) := j_\Omega(u) \).

Local estimators?

Consider \( \Omega_h \) a triangulation of \( \Omega \). Notice we could rewrite \( J(u) \) as \[ J(u) = \sum_{ T\in \Omega_h } \kappa_T j_T(u), \] for appropiated weights \(\kappa_T = \frac{\lVert u \rVert_T^2}{\lVert u \rVert_\Omega^2} \).

Potential problems

Despite the theory, seems hard to believe these functionals could work as their derivation holds for the exact solution. Might be worth to analize the weak formulation and the functional setting. A restriction operator \(r_\omega: H^1(\Omega) \to H^1(\omega) \) could play a role in this context. On the other hand, a global phenomenon seems unlikely to be described by a local behavior. We might want to experiment and analize the behaviour of the computed values.

Error control for elliptic problems

There is a general framework for optimal control to a posteriori error estimations in FEM by Rannacher and Becker, and further work in elliptic eigenvalue problems by Rannacher and Heuveline that suggesting how to proceed.

Setting of the problem

Let \(\mathcal{a}(\cdot,\cdot)\) the bilinear form associated to the weak formulation of the Helmholtz operator. We want to find a pair \((v,\lambda) \in H \times \mathbb{C} \) such that \[ a(v,\phi) = \lambda (v,\phi), \text{ for all } \phi \in H, \] where \( (\cdot, \cdot) \) is the standard \(L^2\) inner product, and \( H \) is an appropiate Hilbert space. There is a natural counterpart of previous equation, that is given by the Galerkin finite element approximation of the previous equation (i.e., replace \(H\) with \(H_h \subset H\) ). We can add a normalization constrain \( \lVert v \rVert = 1 \) to the above problem.

Define \( V:= H\times \mathbb{C} \) and \( A: V \times V \to \mathbb{C} \) given by \[ A( (v,\lambda); (\psi,\chi) ) := -a(v,\psi) + \lambda(v,\psi) + \overline{\chi} \left( \lVert v \rVert^2 -1 \right). \] Provided a functional \( J: V \to \mathbb{R} \), we want to compute a target value \(J(v,\lambda)\) subject to \(A((v,\lambda);(\psi, \chi))=0\).

Lagrangian!

We seek to characterize the stationary points of \( L: V\times V \to \mathcal{C} \) given by \[ L((v,\lambda);(\psi, \chi)) := J((v,\lambda)) - A((v,\lambda);(\psi, \chi)). \]

This induces a second problem \( DA((v,\lambda);(\Psi, \Chi), (\psi, \chi))= DJ( (v,\lambda);(\Psi, \Chi)) \) for all \( (\Psi, \Chi)\in V \), that requires to be analized for this context. This problem is the dual problem. We need to ensure the solvability of this dual problem, thus we might need to analyze the structure of \(J\).

Main result estimators

There is a big family of residuals that need to be considered.

Global residuals

Consider, for simplicity \( u := (v,\lambda) \) and similarly for (z \in V). Then define \[ \rho(u_h; \cdot) := F(\cdot) - A(u_h;\cdot), \] \[ \rho^\star(z_h; \cdot) := DJ(u_h;\cdot) - DA(u_h;\cdot, z_h). \]

This give us an error representation \( J(u) - J(u_h) \approx \frac{1}{2} \rho(u_h;z-\psi_h) + \frac{1}{2} \rho^\star(z_h;u-\phi_h) \).

Local residuals

For practical purposes, we define local residuals. For this goal, consider \(\mathcal{A}\) to be the differential operator defined by the governing equation. Then, \[ \rho_T := \left( \lVert\mathcal{A}v_h - \lambda_h v_h \rVert_T^2 + \frac{1}{2}h_T^{-1}\lVert a[\partial_n v_h] \rVert_{\partial \Omega \setminus \partial T}^2 \right)^2, \] \[ \rho^\star_T := \left( \lVert\mathcal{A}^\star v^\star_h - \lambda^\star_h v^\star_h \rVert_T^2 + \frac{1}{2}h_T^{-1}\lVert a[\partial_n v^\star_h] \rVert_{\partial \Omega \setminus \partial T}^2 \right)^2. \] These residuals are obtained after imposing the normalization constrains and breaking the global residuals into element-wise integrands.

For employing these estimates, we also require an interpolator \( i_h: H \to H_h \) verifying some interpolation inequalities.

Estimators

Consider a saturation assumption (i.e., we are able to approximate the eigenvalue problem) given by \[ \Delta_h := \max\left( \lvert \lambda - \lambda_h \rvert, \lVert v - v_h \rVert, \lVert v^\star - v^\star_h \rVert\right) < 1. \] Then under a specific choice of the goal functional \(J\) we obtain \[ \lvert\lambda - \lambda_h\rvert \lesssim \sum_T h_T^2\left( \rho_T^2 + {\rho_T^{\star}}^2 \right). \]

Nevertheless, there is a more important control for the convergence of the eigenvalues. \[ \lvert\lambda - \lambda_h\rvert \leq \sum_T h_T^2\left( \rho_T w_T^{\star} + \rho_T^{\star} w_T \right). \]

The quantities \( w_T, w_T^\star\) are bounded by the Hessians of the eigenvalue problems. This means \( w_T \lesssim \lVert \nabla^2 v \rVert \) and \( w_T^\star \lesssim \lVert \nabla^2 v^\star \rVert \).

Heuristically, it makes sense to replace these weights with \( \tilde w_T := \lVert \nabla^2_h v_h \rVert_T \) and \( \tilde {w_T}^\star := \lVert \nabla_h^2 {v_h}^\star \rVert_T \).

Error-balancing stratergy

The main idea is to balance the error among the elements. Define \[ \eta_T := h_T^2 \left( \rho_T \tilde w_T^\star + \rho_T^\star \tilde w_T \right). \] Then, as \( \lvert J(v,\lambda) - J(v_h, \lambda_h)\rvert \approx \sum_T \eta_T\), we mark/flag the elements \(T\) such that \(\eta_T\) does not meet certaing tolerance criteria. For example, we will refine the triangle \(T\) if \[ \eta_T \geq \frac{\mathrm{TOL}}{ \lvert \mathbb{T}_h \rvert }. \]