Mathematical background

Full description of the method including proofs and examples you can find in PDF.

1. Finite element method

We seek eigenpairs \((\omega^2, u)\) (note that we use the notation \(\omega^2\) for the eigenvalue) of a linear partial derivative operator \(L\) on a domain \(\Omega\) with given boundary conditions. This means we aim to find \(u\) such that \(Lu = \omega^2 u\). We begin with the standard Finite Element Method (FEM) approach: defining a mesh, fixing a finite-dimensional solution space of hat functions (or optionally, splines) \(V_h\), and reformulating the problem in its weak form.

As an example, we consider the negative Laplacian eigenvalue problem with Neumann boundary conditions:

\[-\Delta u = \omega^2 u \quad \text{ in } \Omega,\]
\[\frac{\partial u}{\partial \nu} = 0 \quad\text{ on } \partial\Omega.\]

Its weak form is:

\[\int_\Omega \nabla u \cdot \nabla \varphi \, dx = \omega_h^2 \int_\Omega u \varphi \, dx.\]

Using a fixed basis \(\varphi_1, \dots, \varphi_N\) of the solutions space \(V_h\), we define the discretization matrices \(S\) and \(M\) as follows:

\[s_{ij} := \int_\Omega \nabla \varphi_i \cdot \nabla \varphi_j \, dx \quad \text{ and } \quad m_{ij} := \int_\Omega \varphi_i \varphi_j \, dx.\]

This leads to the discrete matrix eigenvalue problem:

\[Sv = \omega^2 Mv,\]

or, equivalently,

\[M^{-1}Sv = \omega^2 v,\]

where \(v\) denotes the coordinate vector of the eigenfunction.

2. Filtered time-domain solutions

Krylov eigenvalue solver finds eigenvalues within a specified region of interest, denoted as \((\omega_\min^2, \omega_\max^2)\). We aim to construct a linear operator \(C\) that shares the same eigenspaces as \(M^{-1}S\), but with different eigenvalues. Crucial to the contruction of the operator \(C\) is the weight vector \(\alpha\) and induced discrete filter function (dff) \(\beta(\omega)\).

  • If \(v\) is an eigenvector of \(M^{-1}S\) corresponding to an eigenvalue \(\omega^2 \in (\omega_\min^2, \omega_\max^2)\), then \(v\) is an eigenvector of \(C\) corresponding to a large eigenvalue \(\beta(\omega)\).

  • If \(v\) is an eigenvector of \(M^{-1}S\) corresponding to an eigenvalue \(\omega^2 \notin (\omega_\min^2, \omega_\max^2)\), then \(v\) is an eigenvector of \(C\) corresponding to a close to zero eigenvalue \(\beta(\omega)\).

_images/dff.png

For a detailed description of this operator, see PDF. It is impossible to push values of \(\beta\) outside the region of interest close to 0, so we control the values of dff withinin a control interval denoted as \((0, \omega_{\mathrm{end}}^2)\) only. All eigenvalues of \(M^{-1}S\) should lie in this interval: \((0, \omega_{\mathrm{end}}^2)\).

In constructing \(C\), we perform simple time-stepping with \(L\) time-steps of size \(\tau\) to the end-time \(T = L\tau\). The CFL condition requires, that \(\tau \leq 2/\omega_{\mathrm{end}}\). For computational efficiency, we recommend using \(\tau \lessapprox 2/\omega_{\mathrm{end}}\). A higher number of time-steps \(L\) (or larger end-time \(L\)) increases the reliability of the algorithm and improves the behavior of the dff, but linearly increases computation costs.

3. Krylov iteration

To find the eigenpairs of the operator \(C\), we use the Krylov iteration method, starting from a random vector. With each iteration step, the algorithm “catches” one eigenpair of \(C\), beginning with the largest eigenvalues. increases accuracy of already catched eigenpairs. Simultaneously, it refines the accuracy of the previously computed eigenpairs. Once the eigenvectors of \(C\) are determined, the eigenvectors of the original problem, \(M^{-1}S\), are also known, allowing for the computation of the corresponding eigenvalues, \(\omega^2\).

The properties of the Krylov iteration, along with an appropriate filtering method, ensure that the desired eigenvalues are captured early in the process. This significantly reduces the computational cost compared to other methods.

_images/convergence.png

The results of the first 50 steps of the Krylov iteration for the negative Laplacian problem with the region of interest \((\omega_\min, \omega_\max) = (11, 13)\). In the plot, each horizontal level corresponds to one iteration step (\(k\) denotes the step number). Each point represents a computed value of \(\omega\), with its color and marker indicating the accuracy of the approximation. Eigenvalues computed with an accuracy better than \(10^{-5}\) are marked with an “o”. The vertical dotted lines represent true eigenvalues of \(M^{-1}S\) for reference.

The red line (scaled on the right axis) represents the dff \(\beta(\omega)\). Faster convergence is observed in regions where the dff is larger, while little to no convergence occurs in other regions. After 50 iterations, all eigenvalues within the target interval are approximated with an accuracy exceeding \(10^{-5}\). A direct approach would require solving a 629-dimensional matrix eigenvalue problem.