A rational approximation method for solving acoustic nonlinear eigenvalue problems

TL;DR

This paper introduces two efficient rational approximation methods, using Cauchy integral and Chebyshev interpolation, for solving large-scale nonlinear acoustic eigenvalue problems derived from boundary element methods.

Contribution

It develops two novel approximation techniques with a parallelizable Rayleigh-Ritz procedure for large-scale acoustic eigenvalue problems.

Findings

Effective for problems with up to two million degrees of freedom.

Demonstrated high accuracy on benchmark and industrial applications.

Methods outperform traditional approaches in computational efficiency.

Abstract

We present two approximation methods for computing eigenfrequencies and eigenmodes of large-scale nonlinear eigenvalue problems resulting from boundary element method (BEM) solutions of some types of acoustic eigenvalue problems in three-dimensional space. The main idea of the first method is to approximate the resulting boundary element matrix within a contour in the complex plane by a high accuracy rational approximation using the Cauchy integral formula. The second method is based on the Chebyshev interpolation within real intervals. A Rayleigh-Ritz procedure, which is suitable for parallelization is developed for both the Cauchy and the Chebyshev approximation methods when dealing with large-scale practical applications. The performance of the proposed methods is illustrated with a variety of benchmark examples and large-scale industrial applications with degrees of freedom varying from several hundred up to around two million.

Tables

Table 1. Table 1: Approximations of the 6 6 6 smallest eigenvalues (including multiplicities) of the 3D Laplace eigenvalue problem on Ω = [ 0 , 1 ] 3 Ω superscript 0 1 3 \Omega=[0,1]^{3} with homogeneous Dirichlet boundary conditions  [ 7 , Table 1] .

no.	eigenvalue	multiplicity
1	5.441398	1
2	7.695299	3
3	9.424778	3
4	10.419484	3
5	10.882796	1
6	11.754763	6

Table 2. Table 2: Exact eigenvalues of the 3D Laplace eigenvalue problem on a unit sphere with homogeneous Dirichlet boundary conditions.

no.	eigenvalue	multiplicity
1	3.1416	1
2	4.4934	3
3	5.7634	5
4	6.2831	1
5	6.9879	7
6	7.7252	3

Equations

$$\Delta u(x)+k^{2}u(x)=0,x\in \Omega \subset {\mathbb{R}}^3,$$

$$a(x)u(x)+b(x)\frac{\partial u(x)}{\partial n}=0,x\in \partial \Omega ,$$

$$C(x)u(x)={\int }_{\partial \Omega }\left(\frac{\partial g}{\partial n_y}u(y)-g(\Vert x-y\Vert )\frac{\partial u}{\partial n_y}\right)dy,$$

$$g(\Vert x-y\Vert )=\frac{e^{iz\Vert x-y\Vert }}{4\pi \Vert x-y\Vert }.$$

$$T(z)u=0,T(z):=AH(z)-BG(z),$$

$$T(\lambda )u=0,$$

$$T(z)=\frac{1}{2ı\pi }{\int }_{\Gamma }\frac{T(t)}{t-z}dt.$$

$$\tilde{T}(z)\approx \sum _{i=0}^m\frac{{\omega }_{i}T({\sigma }_i)}{z-{\sigma }_i},$$

$$\tilde{T}(z)$$

$$=B_0{f}_0(z)+B_1{f}_1(z)+\dots +B_m{f}_m(z),$$

$$v_i:=f_i(z)u,\text{mbox}fori=0,\dots ,m.$$

$$\tilde{T}(z)u=B_0{v}_0+B_1{v}_1+\dots +B_m{v}_m=0.$$

$$\tilde{T}(z)=B_0{\tau }_0(z)+B_1{\tau }_1(z)+\dots +B_m{\tau }_m(z),$$

$$\mathcal{A}w=\lambda \mathcal{M}w.$$

$$\mathcal{A}=\left[\begin{array}{ccccc}{\sigma }_{0}I& & & & I\\ & {\sigma }_{1}I& & & I\\ & & \ddots & & ⋮\\ & & & {\sigma }_{m}I& I\\ -B_0& -B_1& \cdots & -B_m& 0\text{par}\end{array}\right],\mathcal{M}=\left[\begin{array}{ccccc}I& & & & \\ & I& & & \\ & & \ddots & & \\ & & & \ddots & \\ & & & & 0\end{array}\right],$$

$$\mathcal{A}=\left[\begin{array}{ccccc}0& I& & & \\ I& 0& I& & \\ & \ddots & \ddots & \ddots & \\ & & I& 0& I\\ -B_0& \cdots & -B_{m-3}& C_{m-2}& -B_{m-1}\text{par}\end{array}\right],\mathcal{M}=2\left[\begin{array}{ccccc}\frac{1}{2}I& & & & \\ & I& & & \\ & & \ddots & & \\ & & & I& \\ & & & & B_m\end{array}\right],$$

$${\hat{B}}_0{f}_0(z)y+{\hat{B}}_1{f}_1(z)y+\dots +{\hat{B}}_m{f}_m(z)y=0,$$

$$T_U(z)={\hat{B}}_0{f}_0(z)+{\hat{B}}_1{f}_1(z)+\dots +{\hat{B}}_m{f}_m(z).$$

$${\hat{B}}_0{\hat{v}}_0+{\hat{B}}_1{\hat{v}}_1+\dots +{\hat{B}}_m{\hat{v}}_m=0,$$

$$\mathcal{A}w=\lambda \mathcal{M}w,$$

$$\mathcal{A}=\left[\begin{array}{ccccc}{\sigma }_{1}I& & & & I\\ & {\sigma }_{2}I& & & I\\ & & \ddots & & ⋮\\ & & & {\sigma }_{m}I& I\\ -{\hat{B}}_1& -{\hat{B}}_2& \cdots & -{\hat{B}}_m& 0\text{par}\end{array}\right],\mathcal{M}=\left[\begin{array}{ccccc}I& & & & \\ & I& & & \\ & & \ddots & & \\ & & & \ddots & \\ & & & & 0\end{array}\right]$$

$$\mathcal{A}=\left[\begin{array}{ccccc}0& I& & & \\ I& 0& I& & \\ & \ddots & \ddots & \ddots & \\ & & I& 0& I\\ -{\hat{B}}_0& \cdots & -{\hat{B}}_{m-3}& {\hat{C}}_{m-2}& -{\hat{B}}_{m-1}\end{array}\right],\mathcal{M}=\left[\begin{array}{ccccc}I& & & & \\ & 2I& & & \\ & & \ddots & & \\ & & & 2I& \\ & & & & 2{\hat{B}}_m\end{array}\right]$$

$$\mathcal{H}w=\frac{1}{\lambda -\sigma }w,\mathcal{H}:={\left(\mathcal{A}-\sigma \mathcal{M}\right)}^{-1}\mathcal{M}.$$

$$\left[\begin{array}{cc}D& F\\ B& C\end{array}\right].$$

$$L=\left[\begin{array}{cc}I& 0\\ B{D}^{-1}& I\end{array}\right],U=\left[\begin{array}{cc}D& F\\ 0& S\end{array}\right],$$

$$\left[\begin{array}{cc}D& F\\ B& C\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}a\\ b\end{array}\right].$$

$$\mathcal{A}{w}^{(j+1)}=y^{(j)}\text{mbox}with{y}^{(j)}=\mathcal{M}{w}^{(j)}\text{mbox}and{w}^{(j)}=[v^{(j)};u^{(j)}],$$

$$S{u}^{(j+1)}$$

$$D{v}^{(j+1)}$$

$$v_i^{(j+1)}=\frac{v_i^{(j)}-u^{(j+1)}}{{\sigma }_i},i=0,\dots ,m.$$

$$S=-\sum _{i=0}^m\frac{B_i}{{\sigma }_i},\text{and}b=u^{(j)}-\sum _{i=0}^m\frac{B_i}{{\sigma }_i}{v}_i^{(j)}.$$

Full text

1 Background and Introduction

The Boundary Element Method (BEM) is a powerful approach developed to solve integral equations [16]. The idea of applying the BEM in many branches of science and engineering has gained popularity over the past few years, e.g., in elasticity, ground and water flow, wave propagation and in electromagnetic problems [19]. The most commonly used approaches for numerically solving PDEs are the Finite Difference Method (FDM) and the Finite Element Method (FEM). A standard finite difference method is suitable when dealing with simple domains (e.g. rectangular grids), while the finite element method can handle more complex domains. However, much work has to be done to numerically dicretize a whole computational domain (generate meshes) and this task becomes even more difficult when dealing with complicated domains in higher dimensions, i.e., $d\geq 3$. This is where BEM becomes appealing because it allows to significatly reduce the overall computational complexity of the solution process. Instead of solving a problem for the partial differential operator defined on the whole domain $\Omega$, the boundary element method uses an associated boundary integral equation reducing the domain of the problem to the boundary $\partial\Omega$. This comes at a cost since the matrix problem to solve in the approximation becomes dense.

In the following, we are interested in the efficient solution of nonlinear eigenvalue problems (NLEVPs) resulting from the boundary element (BE) discertization of the acoustic problems. Although a finite element discretization of the problem yields a generalized (linear) eigenvalue problem, it requires a discretization of the whole domain $\Omega$ which is not always feasible, e.g., if the domain is unbounded. Though the topic of NLEVPs built upon the boundary element method (BEM) has been around for a number of years, the lack of efficient eigensolvers has delayed a full exploration of BE–based approaches. Recently, eigenvalue solvers based on contour integrals were developed and this made BEM an attractive alternative to the usual contenders when solving challenging nonlinear eigenvalue problems [15, 22, 23]. Contour based methods have the ability to solve NLEVPs when the eigenvalues of interest lie inside a given closed contour in the complex plane using rational or polynomial approximation [7, 9, 4, 21, 13]. Despite these efforts, solving NLEVPs is still a computationally intensive task. Assembling interpolation matrices and solving linear systems in the BE framework are already very expensive due to the unstructured, dense and complex nature of the resulting matrices. For example, the Chebyshev interpolation of the BE formulation of the large-scale accoustic problem discussed in [7] results in a generalized eigenvalue problem which cannot be easily handled with the state-of-the-art linear solvers. Another drawback of this method is that the quality of the approximations quickly deteriorates when dealing with complex eigenvalues.

It is the purpose of this paper to overcome the aforementioned difficulties and develop eigensolvers suitable for calculations of eigenvalues of NLEVPs arbitrarly located in the complex plane. The paper illustrates the performance of the proposed method with a problem that arises in the modal analysis of large-scale acoustic problems.

Consider the three-dimensional (3D) acoustic Helmholtz equation

[TABLE]

where $\Delta$ is the Laplace operator, $u(x)$ is the sound pressure at point $x$, $k=\omega/c$ is the wave number with the circular frequency $\omega$ and the speed of sound $c$ through the fluid medium. Equation (1) is subject to a homogeneous condition on its boundary $\partial\Omega$ of the form

[TABLE]

where $\frac{\partial}{\partial n}$ denotes the outward normal to the boundary at point $x$.

Using BEM yields the following Helmholtz integral equation [11]

[TABLE]

where $C$ denotes the solid angle at point $x$, $n_{y}$ the surface unit normal vector at point $y$ and $g(\cdot)$ the free-space Green’s function [6, 3]

[TABLE]

The continuous Helmholtz integral equation (3) can be discretized to form the following discrete problem from which the unknown boundary node values $z$ can be determined,

[TABLE]

where $A$ and $B$ are diagonal matrices related to the functions $a(\cdot)$ and $b(\cdot)$ in (2), $H$ and $G$ are the matrices containing the coefficients related to the integrals on the surface of $\frac{\partial g}{\partial n_{y}}$ and $g$, respectively [11]. The boundary integrals are discretized by Gauss-Legendre quadrature where the singularities of Green’s function and its derivative are isolated in the integral of revolution, and the integrations are performed analytically using sums of elliptic integrals [12]. Here, $T(z)\ \in\mathbb{C}^{n\times n}$ is a matrix function that is nonlinear is $z$ and holomorphic since the free-space Green’s functions are holomorphic functions of $z$. Obviously, equation (5) is a NLEVP of the general form

[TABLE]

and the objective of this paper is to develop methods for finding all eigenvalues $z$, satisfying (5), that are located inside a certain region of the complex plane enclosed by the contour $\Gamma$.

2 Rational and Chebyshev approximation methods for NLEVPs

The first method we consider is adapted from [14] and it is based on the Cauchy’s integral formula. Given a Jordan curve $\Gamma$ that surrounds the eigenvalues of interest, we express the matrix function $T(z)$ as follows:

[TABLE]

By replacing both occurrences of $T(\cdot)$ in (7) by $T_{ij}(\cdot)$, one can see that the above expression is equivalent to expressing each individual entry $T_{ij}(z)$ of $T(z)$ by the Cauchy integral formula. Equality (7) is valid for $z$ inside the contour $\Gamma$ and the only requirement is that $T(z)$ be analytic inside the contour. As is classically done [10] we use a numerical quadrature formula to obtain the following Cauchy integral approximation $\widetilde{T}(z)$ of $T(z)$

[TABLE]

where the $\sigma_{i}$’s are quadrature points located on the contour $\Gamma$ and the $\omega_{i}$’s the corresponding quadrature weights.

Setting $B_{i}=\omega_{i}T(\sigma_{i}),$ equation (8) can be rewritten as

[TABLE]

with $f_{i}(z)=\frac{1}{z-\sigma_{i}},i=0,\ldots,m$. For a given vector $u$ we now define

[TABLE]

Then the approximate nonlinear eigenvalue problem $\widetilde{T}(z)u=0$ yields

[TABLE]

Chebyshev interpolation of order $m$ can also be used to obtain the same form as (10) of the approximation of the matrix-valued function $T(z)$. In this method, proposed in [7], the function $T(z)$ is expanded using a degree $m$ Chebyshev polynomial expansion of the form [1] :

[TABLE]

where $B_{i}$ and $\tau_{i}(z)$ are coefficient matrices and Chebyshev basis functions, respectively. The corresponding nonlinear eigenvalue problem is of the same form as (11) with the vectors $v_{i}$ now defined by $v_{i}=\tau_{i}(z)u$.

The problem (11) for the Cauchy interpolation, and its Chebyshev interpolation counterpart, can be reformulated as a generalized linear eigenvalue problem:

[TABLE]

For the Cauchy rational approximation we have:

[TABLE]

and $w=\begin{bmatrix}v_{0}^{T},v_{1}^{T},\ldots,v_{m}^{T},u^{T}\end{bmatrix}^{T}$ whereas for the Chebyshev interpolation

[TABLE]

where $C_{m-2}\equiv B_{m}-B_{m-2}$ and $w=\begin{bmatrix}u^{T},v_{1}^{T},\ldots,v_{m-1}^{T},v_{m}^{T}\end{bmatrix}^{T}$.

With regards to the rational approximation described above, we note that an alternative that has been used with some success in the literature is the Barycentric approximation formula [4]. However, our tests with this technique showed no significant improvement in our context over the simple Cauchy formula used above. Note that it is also possible to exploit other polynomials, using different classes of orthogonal polynomials but we will restrict our attention to Chebyshev polynomials of the first kind. Finally note that Chebyshev approximation works best for eigenvalues located in an interval while the Cauchy rational approximation is suitable for general complex spectra.

3 Rayleigh–Ritz procedure for BEM eigenvalue problem

Let $U$ be a basis of dimension $\nu$ of a subspace that contains good approximations of the eigenvectors of the NLEVP problem (5). Then, it is possible to apply a Rayleigh-Ritz procedure to (5) to obtain approximate eigenpairs. The approximate eigenvector will be of the form $u=Uy$  with  $y\in\mathbb{C}^{\nu}$. Then expressing that $T(z)u$ is orthogonal to the range of $U$ yields the projected problem $U^{H}T_{U}(z)u=0$ or,

[TABLE]

where $\widehat{B}_{i}=U^{H}B_{i}U$. We will denote $T_{U}(z)$ the projected operator, namely,

[TABLE]

Then, applying the same procedure as before to the projected problem we see that (16) becomes:

[TABLE]

with $\widehat{v}_{i}=\frac{y}{z-\sigma_{i}}$ in the case of rational approximation and $\widehat{v}_{i}=\tau_{i}(z)y$ when a Rayleigh-Ritz procedure is applied to the Chebyshev interpolation.

3.1 Solution of the reduced NLEVP

Analogously to what was discussed in Section 2, the problem (18) for the Cauchy rational approximation, as well as its Chebyshev interpolation counterpart, can be written down in a block form (13), but now of much smaller dimension. The projected nonlinear problem (18) yields the following linearized problem

[TABLE]

with $w=\begin{bmatrix}\widehat{v}_{0}^{T},\widehat{v}_{2}^{T},\ldots,\widehat{v}_{m}^{T},y^{T}\end{bmatrix}^{T}$ and

[TABLE]

for the Cauchy rational approximation, and $w=\begin{bmatrix}y^{T},\widehat{v}_{1}^{T},\ldots,\widehat{v}_{m-1}^{T},\widehat{v}_{m}^{T}\end{bmatrix}^{T}$,

[TABLE]

for the Chebyshev interpolation where $\widehat{C}_{m-2}=\widehat{B}_{m}-\widehat{B}_{m-2}$. If $\nu$ is fairly small, the problem (19) can be solved directly, i.e., using standard dense packages. When $\nu$ is larger, the problem must be handled differently by some iterative procedure. Since for BEM problems the matrices $B_{i}$ are generally complex, dense and unstructured, solving these linear eigenvalue problems can be computationally expensive. Therefore, it may be advantageous to rely on subspace iteration or an Arnoldi-type method to solve (19).

Note that for both the Cauchy rational approximation and the Chebyshev interpolation method, the matrices ${\cal A}$ and ${\cal M}$ don’t have to be formed explicitly. If the partial solution of the problem (19) are of interest, effective methods such as the Implicitly Restarted Arnoldi method can be used to find a few of the extremal eigenvalues. Unfortunately, these methods become expensive when the eigenvalues of interest are deep inside the spectrum.

Alternatively, we can solve the interior eigenvalue problem with the help of the shift-and-invert technique, which replaces the solution of the generalized eigenvalue problem (19) by the following problem

[TABLE]

Using the Arnoldi or the subspace iteration method to extract extremal eigenvalues of (22) will result in approximations of the eigenvalues of (19) closest to $\sigma$. Again, the matrix $\mathcal{H}$ need not be formed explicitly to compute the matrix-vector product $y=\mathcal{H}x$. Instead, we can use a simple $LU$ factorization that takes advantage of the sparsity of ${\cal A}$ and ${\cal M}$. First, note that the matrix $\left({\cal A}-\sigma{\cal M}\right)$ is of the form

[TABLE]

By exploiting the sparsity of the matrices $D$ and $F$, we can easily form the following $LU$ factorization

[TABLE]

where $S=C-BD^{-1}F$ is known as the Schur complement of the block $C$. With the use of matrix $S$, we can use the Arnoldi algorithm on vectors of shorter length. Solving the shifted and inverted problem (22) with Arnoldi algorithm requires solving linear systems of the form

[TABLE]

Using the Schur complement $S$, $y$ can be easily obtained by solving $Sy=b-BD^{-1}a$ and since $D$ is a diagonal matrix and $F$ is a block of identity matrices, one can determine $x$ by using the relation $Dx+Fy=a$.

3.2 Construction of the subspace of approximants

We begin this section by noting that the Arnoldi-type or subspace iteration methods discussed in the previous section can be applied to a linear eigenvalue problem ${\cal A}w=\lambda{\cal M}w$ obtained directly from (11). However, proceeding in this way would require either solving linear eigenvalue problems of size $mn+n$ when using Arnoldi-type methods, or storing vectors of length $mn+n$ in the subspace iteration method, and this can be computationally expensive when $m$ is large. Therefore, it is important to develop a technique that allows to work with subspaces of smaller dimensions that requires storing shorter vectors. A procedure of this type, which works with subspaces of dimention $m$ is presented next.

Let us first consider a large linear eigenproblem of the form (13) obtained from (11) without a projection. To introduce the approach that works with vectors of dimension $n$, we first point out that for an approximate eigenpair $(\lambda,u)$, $u$ is the bottom (resp. top part) of an approximate eigenvector $w$ of the large linear eigenvalue problem (13) associated with (10) for the Cauchy rational approximation, (resp. (12) for the Chebyshev interpolation). Let $W^{(0)}$ be a random initial set of $\nu$ basis vectors of a certain subspace, where each of the $\nu$ columns is of the form 111Here we use Matlab notation: $[v;u]$ is a vector that stacks $v$ on top of $u$. $w=[v;u]$ (resp. $w=[u;v]$) for the splitting associated with the Cauchy rational approximation (resp. Chebyshev interpolation). Next, in order to make these initial random vectors close to the eigenvectors of interest, we apply $q$ steps of the inverse power method with matrix $\mathcal{H}$ in (22) to each column of $W^{(0)}$ separately. A subspace of dimension $n$ that approximates the eigenvectors of (11) is then obtained from the bottom parts (resp. top parts) of the processed columns for the Cauchy rational approximation (rep. Chebyshev interpolation). Although this process involves the column vectors of $W^{(0)}$, only vectors of length $n$ need to be saved and the iterates $v$ can be discarded. The accuracy of the extracted eigenpairs obtained from applying a Rayleigh-Ritz projection can be further refined by updating $U$ in a process that takes advantage of the structure of the approximate eigenvectors. Let $(\lambda,u)$ be an approximate eigenpair of (11) obtained from applying a Rayleigh-Ritz projection using $U$. The new redefined vector $w$ for each interpolation method is discussed next. For the Cauchy rational approximation, the vector $v=[v_{1};v_{2};...;v_{m}]^{T}$ (the top part of vector $w$), is obtained by setting $v_{i}=\frac{u}{\lambda-\sigma_{i}}$, whereas for the Chebyshev interpolation (the bottom part of vector $w$) it is defined by setting $v_{i}=\tau_{i}(\lambda)u$.

3.3 The inverse power method

The straightforward linearizations (14) of the Cauchy rational approximation and (15) of the Chebyshev interpolation, discussed in Section 2, are high dimensional problems and they become computationally demanding as the order $m$ of the approximations grows. The Rayleigh-Ritz approach discussed above is inexpensive even if $m$ is large. The biggest computational task of the presented Rayleigh-Ritz projection lies in performing $q$ steps of the inverse power method with the matrix $\left({\cal A}-\sigma{\cal M}\right)^{-1}{\cal M}$. It is the purpose of the following discussion to show how each step of the inverse power method can be carried out inexpensively. For simplicity, we will assume that the shift $\sigma$ is the center of the unit circle (resp. interval $[-1,1]$) for the Cauchy rational approximation (resp. Chebyshev interpolation). This is a natural choice, since any circle in the complex plane can be scaled to the unit circle and any real interval $[a,b]$ can be scaled to the interval $[-1,1]$. Throughout this discussion, the superscript $j$ will correspond to the iteration number, while the subscript $i$ will correspond to the blocks of the vectors $v^{(j)}$. We begin by discussing the inverse power method for the Cauchy rational approximation.

Inverse power iteration for the Cauchy rational approximation

For the Cauchy interpolation, each step of the inverse power iteration method requires solving a linear system

[TABLE]

which is of the form (25). Therefore, the iterates of the inverse power method can be determined by solving

[TABLE]

Since $D$ is a diagonal matrix and $F$ is a block vector of identity matrices, $v_{i}^{({j+1})}$ are determined by

[TABLE]

Again, exploiting the structure of $D$ and $F$, the iterate $u^{({j+1})}$ can be obtained by solving (27) with

[TABLE]

Algorithm 1 performs one step of the inverse power iteration for the Cauchy rational approximation.

Inverse power iteration for the Chebyshev interpolation

Recall that the iterates obtained from the inverse power method for the Chebyshev interpolation can be written as $w^{({j})}=[u^{({j})};\ v^{({j})}]$. Similarly to the Cauchy rational approximation, each step of the inverse power method requires solving the linear system

[TABLE]

Since $\mathcal{M}$ is a block diagonal matrix, $y^{({j})}=\mathcal{M}w^{({j})}$ can be easily evaluated. The question that remains to be answered is how to solve efficiently the linear system $\mathcal{A}w^{({j+1})}=y^{({j})}$. By taking advantage of the block structure of $\mathcal{A}$ for the Chebyshev interpolation, it follows naturally that this problem can be treated by performing the following steps, see [7, Section 2.3]. To compute the bottom part $v^{({j+1})}$ of $w^{({j+1})}$ we will use the recursion

[TABLE]

for odd-numbered blocks and

[TABLE]

for even-numbered blocks. Since the blocks $v_{2i-2}^{(j+1)}$ in (33) are even-numbered, we can further expand the recurrence relation, i.e.,

[TABLE]

where

[TABLE]

Since $v_{0}=\tau_{0}(z)u$ and $\tau_{0}=1$ (zeroth Chebyshev polynomial), $u^{({j+1})}$ is the top part of vector $w^{(j+1)}$, i.e., $u^{(j+1)}=v_{0}^{(j+1)}$ and it can be obtained by solving

[TABLE]

Given the number of quadrature nodes $m$, let us consider the Euclidean division of $m$ by $2$, i.e., $m=2\cdot q+r$, with quotient $q$ and remainder $r$. Then the matrix $G$ has the following form

[TABLE]

The vector $b$ depends on the parity of $m$. If $m$ is odd

[TABLE]

and when it is even, then

[TABLE]

Algorithm 2 implements one step of inverse power method for the Chebyshev interpolation.

To this end, only one LU factorization is required – of the Schur complement matrix $S$ in the case of the Cauchy rational approximation or matrix $G$ for the Chebyshev interpolation – in the preprocessing step for all $q$ steps of the inverse power method.

4 Numerical Experiments

This section will illustrate the behavior of the approaches presented in this paper for solving nonlinear eigenvalue problems in the form (5) resulting from boundary element disretization of (1) – (2). All experiments were performed with Matlab R2018a. Furthermore, computations in Example 3 were performed in parallel on a Linux cluster at the Minnesota Supercomputer Institute that has $32$ cores and $31.180$ GB per-core memory.

For the presented examples, the contour $\Gamma$ is either circular or elliptic and the eigenvalues of interest are those closest to the center of $\Gamma$, i.e., in Algorithms 1 and 2 the shift $\sigma$ is selected to be the center of the region enclosed by the contour $\Gamma$.

In the case of a circular contour, the $m$ quadrature nodes and weights used to perform the numerical integration to approximate the functions $f_{j}$ inside the contour $\Gamma$ were generated using the Gauss-Legendre quadrature rule. To illustrate the effectiveness of the proposed approaches, we compare the eigenvalues obtained by each algorithm either with exact eigenvalues or the approximations obtained by the Beyn’s method [5] or/and via a corresponding linearization.

Example 1

As our first example, we consider the 3D Laplace eigenvalue problem (1) on the unit cube $\Omega=[0,1]^{3}$ with homogeneous Dirichlet boundary conditions, i.e., (2) with $a(x)=1$ and $b(x)=0$. The exact eigenvalues for this problem are known and given by

[TABLE]

We are interested in the six smallest eigenvalues (including multiplicities) of (1) presented in Table 1.

To determine these eigenvalues using the Cauchy approximation technique, we will build the rational approximation of the matrix-valued function $T(\cdot)$ using circular and elliptic contours. First, we compare the accuracy between the Cauchy rational approximation and the Chebyshev interpolation of $T(\cdot)$. Note that since the eigenvalues of (1) with homogenous Dirichlet boundary condition are real we can use the Chebyshev interpolation technique which target situations when the eigenvalues of interest lie in an interval. Figure 1 shows the errors of each approximation versus the order of approximation $m$ for both circular and elliptic contour. For simplicity, all the errors are evaluated on a fine mesh in $[-1,1]$, since arbitrary curves in the complex plane can be parametrized using this interval. From this figure, we can easily see that the errors in the Cauchy rational approximation and Chebyshev interpolation decay exponentially with the order of the approximation $m$, which implies that a moderate $m$ is usually sufficient to reach a good accuracy. To capture the eigenvalues of interest, we first consider a circle of radius $r=3.5$ centered at $c=8.5$. We can then solve the linear eigenvalue problem (13) associated with (11) with $m=25$ trapezoidal quadrature nodes by performing as many steps of shift-and-invert Arnoldi algorithm as needed to extract the $17$ eigenvalues closest to the center $c$. The left hand side of Figure 2 presents the eigenvalues computed by Cauchy approximation and those computed by Chebyshev interpolation on a uniform mesh with $864$ triangles. Note that, in order to make a fair comparison between the two methods, the number of interpolation points for Chebyshev interpolation method is chosen to be the same as the number of quadrature nodes $m$. For Chebyshev interpolation method the real interval enclosing the eigenvalues of interest is chosen as $[5,12]$. The right hand side of Figure 2 illustrates the comparison between the accuracy of the rational approximation and the Chebyshev interpolation. The accuracy of an eigenpair $(\lambda,u)$ is measured by the relative residual $\|T(\lambda)u\|_{2}/\|u\|_{2}$. Note that the accuracy of the rational approximation can be considerably improved by using an elliptic contour instead of a circle. A rational approximation (8) is then built using an elliptic contour centered at $c=8.5$ with semi-major axis $r_{x}=3.5$ and semi-minor axis $r_{y}=0.1$. The left hand side of Figure 3 presents the eigenvalues computed by the Cauchy rational approximation and those computed by the Chebyshev interpolation, and the right hand side of Figure 3 compares the relative residuals of the two methods.

We now repeat the same experiment using the Rayleigh-Ritz procedure for the Cauchy and Chebyshev approximations. To extract the $17$ eigenvalues listed in Table 1, we start with a random subspace $W$ of dimension $\nu=20$. We then apply $q=10$ steps of inverse power method to $W$ to build a subspace of dimension $\nu$, where each column vector is of size $n$. We recall that these vectors are the resulting top parts and bottom parts of the iterates of Algorithm 1 and 2 for the Cauchy and Chebyshev approximations, respectively. The resulting subspace $U$ is then orthogonalized to obtain an orthonormal basis $U$ that can be used to perform Rayleigh-Ritz projection that leads to a small nonlinear eigenvalue problem of size $\nu$. This small problem is then solved by computing the eigenvalues and eigenvectors of the expanded linear eigenvalue problem (19) of size $(m+1)\nu$. The outer iterations of the reduced procedure for the Cauchy approximation are stopped when

[TABLE]

and for the Chebyshev approximation when

[TABLE]

where $X$, $\Lambda$ are the extracted eigenpairs at each iteration, $\|\cdot\|_{F}$ denotes the Frobenius norm and $tol$ the desired tolerance for the convergence. In our experiments, we run as many outer iterations as needed to achieve convergence with a tolerance $tol=10^{-12}$ for both Cauchy and Chebyshev approximations. This tolerance is achieved after $10$ outer iterations for the Cauchy approximation and after $7$ outer iterations for Chebyshev approximation. Furthermore, Figure 4 presents the relative residuals $\|T(\lambda)u\|_{2}/\|u\|_{2}$ for the $17$ computed eigenvalues obtained using each approximation method. We ephasise that the Rayleigh-Ritz approach combined with Cauchy and Chebyshev approximations delivers more accurate eigenpair approximations than those computed by solving the linearized problem obtained directly from the Cauchy and Chebyshev approximation with projection.

Example 2

As a second example, we consider the 3D Laplace eigenvalue problem (1) on a unit sphere with homogeneous Dirichlet boundary conditions. The analytic expressions for the eigenvalues for this geometry are well-known and given as the zeros of the spherical Bessel function of order $\ell$. We are interested in the $6$ eigenvalues of (1) listed in Table 2.

In order to compute the eigenvalues of interest using the rational approximation technique, we consider an elliptic contour centered at $c=5.5$ with semi-major axis $r_{x}=2.5$ and semi-minor axis $r_{y}=0.1$. The left-hand side of Figure 6 presents the eigenvalues computed by the Cauchy rational approximation and those computed by the Chebyshev interpolation with $m=25$ on a uniform mesh with $384$ triangles, whereas the right-hand side of Figure 6 illustrates the accuracy of the two methods. Also in this example, we have tested the reduced subspace iteration given by Algorithm 3. To extract the $20$ eigenvalues of interest, we consider the Cauchy approximation on a circle centered at $c=5.5$ with radius $r=2.5$. The first outer iteration was carried out with a random subspace $W$ of dimension $\nu=25$ to which $q=10$ steps of inverse power method, given in Algorithm 3, were applied. As for Example 1, we run as many outer iterations as needed to achieve convergence with a tolerance $tol=10^{-12}$ for both approximation methods. Note that $q=10$ steps of the inverse power method were applied at each outer iteration. The Cauchy approximation and Chebyshev interpolation methods achieved desired tolerance after $17$ and $11$ outer iterations, respectively. For completeness, we have also computed the corresponding relative residuals $\|T(\lambda)u\|_{2}/\|u\|_{2}$ for the resulting eigenpairs. These relative residuals are shown in Figure 7 for each eigenpair index.

Example 3

In this example, we illustrate the efficiency of the Cauchy approximation technique applied to the nonlinear eigenvalue problem resulting from BE discretization of a real-world problem of industrial relevance. We consider the geometry corresponding to a a pump casing model created by using the Gmsh tool [8]. Several methods have been proposed in the literature to comprehensively study the acoustic behaviors of the pump casing  [20, 17]. The boundary of the pump model displayed in Figure 8 is partitioned into $3\ 479\ 652$ triangles, leading to a nonlinear eigenvalue problem with $1\ 728\ 508$ DoFs. Problems of such large size add another level of difficulty to our methods, for example, we are unable to store the underlying matrices $B_{i}$ in memory. To overcome this we resort to the $\mathcal{H}$-matrix based compression techniques. Specifically, we will use the Gypsilab toolbox library openHMX [2] in order to directly assemble $\mathcal{H}$-matrix compressed versions of matrices $B_{i}$. Here, we consider the boundary element discretization of problem (1) with a rigid boundary, i.e., $b(x)=0$. For this problem, as mentioned before, complex eigenvalues may occur. Hence, only the Cauchy approximation can be used to determine the eigenvalues associated with this problem. Since for this example the analytic expressions of the eigenvalues are not available, the relative residuals of the computed eigenpairs will be used to verify the accuracy of the obtained approximations.

Let the domain for the Cauchy approximation be given as a circular contour $\Omega$ centered at $c=-15i$ with radius $r=12$. To choose a suitable order of approximation, we consider another circle $\Omega_{1}$ inside $\Omega$ with the same center and radius $r_{1}=r/2$ and then increase $m$ until the resulting rational approximation inside $\Omega_{1}$ is accurate enough. The eigenvalue approximations inside $\Omega$ can be obtained using a different, much coarser triangular mesh and running as many steps of Arnoldi algorithm as needed to accurately solve the expanded linear eigenvalue problem (13). We recall that only one $LU$ factorization of the Schur complement $\mathcal{H}$-matrix $S$ is required in a preprocessing step before the actual Arnoldi algorithm is invoked. In Figure 9, we present the approximation errors versus the order of the Cauchy approximation on a fine mesh on $\Omega_{1}$. The right-hand side of Figure 9 shows that a high accuracy of the rational approximation can be reached for $m=24$. We can therefore solve the eigenvalue problem (5) using $m=24$ trapezoidal quadrature nodes. Forming the $24$ matrices $B_{i}$ and performing the matrix-vector multiplications with $B_{i}$ are efficiently parallelized on $32$ cores, where the per-core memory limit is $\approx 31$ GB. The overall computational time is $7.3$ hours. The left-hand side of Figure 9 shows the computed eigenvalues. It turns out that there are $33$ eigenvalues inside the contour $\Omega$. The relative residuals $\|T(\lambda)u\|_{2}/\|u\|_{2}$ associated with the computed eigenvalues are presented in Figure 10 and Figure 11 shows $4$ different modes of the pump model.

We now consider the reduced subspace iteration approach to solve the same problem with $10704$ triangles. To extract the same $12$ eigenvalues displayed in the left-hand side of Figure 9, we start with a random subspace $W$ of size $\nu=20$ and carry out $20$ outer iterations of Algorithm 3 with $q=10$ steps of inverse power method (Algorithm 1) performed at each single outer iteration. Note that Algorithm 3 is suitable for parallelization. In our tests, we have exploited parallelism for computing the matrices $B_{i}$ and $\widehat{B}_{i}=U^{T}B_{i}U$ using $32$ cores. Since Algorithm 1 is applied to each vector separately to obtain a block of vectors $U$, the construction of the approximate subspace at each outer iteration can also be performed in parallel. The relative residuals reached at the end of the $20$ outer iterations of Algorithm 3 are displayed in Figure 10. Figure 11 shows $4$ modes of the model.

Acknowledgements

The authors would like to thank Mohammed Seaid for providing the mesh data for the test in Example 3 of the experiments and to Gmsh team for making their software available. Similarly, calculations for Example 3 in the experiments could not have been carried out without the availability of the Gypsilab toolbox library. The authors benefitted from the hardware resources and support from the Minnesota Supercomputing Institute