# Well-Conditioning Map for Tensor Spectral Solutions for Eigenvalues in 2D Irregular Planar Domains

Authors: O. Guimarães; W. Labecca; José R. C. Piqueira
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This study proposes a diffeomorphic map for planar irregular domains, which transforms them into a unitary circle suitable for implementing spectral solutions on a complete and orthogonal basis. The domain must be enclosed by a Jordan curve that is univocally described in polar coordinates by aC2 class function. Because of the achieved precision, general eigenvalue problems can be solved, yielding hundreds of eigenvalues and eigenfunctions in a single shot with at least ten-digit precision. In a constructive approach, there is no need for special functions or roots. The examples focus on the Helmholtz equation, but the method applies to other eigenvalue problems such as Hamiltonians. For validation, we compare all the results with the values obtained by finite elements (FE), some with other studies, and discuss sharp error estimate criteria and the asymptotic behavior of the eigenvalues. To demonstrate the generality of the approach, we present several examples with different geometries in which solutions emerge naturally. Considering the spectral accuracy, half of the examples are unpublished solutions, and the other half agree with or expand the precision of previous studies using different approaches.


---

## Body

## 1. Introduction

Eigenvalue problems have applications in several areas of mathematical physics, engineering, and, more recently, transport phenomena in the biological sciences. [1]. Since Orszag’s pioneering work k2, k3, the application of spectral methods to obtain numerical solutions of equations to partial derivatives (PDE) [2, 3] has established itself as the standard for providing, in rectangular domains, solutions with relative errors practically at the limit of the precision of the standard floating-point [4].For irregular domains, the flexible finite element (FE) method, although not spectrally accurate, provides solutions with sufficient precision for many practical engineering problems, solutions in nontensor domains, and is available in proprietary packages such as the MATLAB® 2020 toolbox [5] for PDEs. This package was used to validate the results. For the Helmholtz equation in particular, even in nontensor domains, some researchers [6–11] have determined eigenvalues for the operator with the homogeneous Dirichlet condition with dozens of digits using the fundamental solutions method (MFS), taking advantage of the inequality demonstrated by Fox et al. [12]. In the early 1980s, Orszag [13] proposed a constructive approach to PDEs in irregular domains using diffeomorphism and presented some spectral solutions in boundary value problems with polar coordinates via Chebyshev polynomials for the radial basis and the Fourier series for the azimuthal coordinates.Once the diffeomorphic transformation is established, the equation is transformed by evaluating the effect of the diffeomorphism on the function and its derivatives. Often, this work is lengthy and tedious. Depending on the transformation, mixed derivatives can appear in a simple Laplacian, variable coefficients can emerge by multiplying derivatives, and a sign error can ruin the final result. For example, equation 5.5 [in [13], p.82] presents this problem, which, to the best of our knowledge, has not been discussed in literature. Nowadays, we can circumvent these snags using symbolic computation to automate the process, as we will discuss here.The approach developed here is not intended to compete with the MFS, the flexibility of the FE, or even the extensive method LNKM (local nonsingular knot method) [14]. Instead, we aim to fill a gap in the numerical solutions of PDEs, reaching spectral accuracy through a constructive process for irregular domains delimited by a Jordan curve, uniquely describable in polar coordinates, and sufficiently derivable (class C2) because we restricted the study to 2nd order PDEs.Since the derivative coefficients are considered to be variable, the method applies to several classes of PDEs (elliptical, hyperbolic, or parabolic).Section2.1 describes the diffeomorphic map of the transformation, which bi-univocally associates all points of the irregular domain to the unit circle, except for the origin. We call this map the rational radial map (RRM) and discuss the impact of this transformation on the function and its derivatives.In Section3, we briefly review the well-known orthogonal bases of Chebyshev and Fourier [15], with their operational matrices, since they will be applied in the following section.The resolution of2nd order PDEs in the unit circle by the pseudospectral method is shown in Section 4, explaining the expressions to be applied in the boundary conditions and the construction of an optimal grid so that we have the same resolution, that is, points per wavelength, for both coordinates in the numerical approach. In Section 5, we discuss the analysis of the error, both in eigenvalues and eigenfunctions, as derived by Still [16] with some adaptation. The asymptotical behavior is commented upon for the Dirichlet condition, in addition to some inequalities for the eigenvalues.Finally, in section6, we present examples compared with some results from other authors and with the values obtained by FE using the MATLAB® 2020 toolbox. In addition, we analyze the asymptotic behavior of the numerical solutions and their coherence with some theoretical inequalities regarding eigenvalues in the Dirichlet case.


## 2. Analytical Foundation

The following two sections discuss the application of differential geometry to arrive at the numerical structure for eigenvalue problems in irregular domains.


### 2.1. Radial Reciprocal Map (RRM)

Let the open subsetΩ⊂ℝ2 be a simply connected domain with border ∂Ω=Γc described by a Jordan curve parameterized at the angle θr,θ=Γθ,θ, and we consider a diffeomorphism between open subsets of IR2, as given by the coordinate chart(1)Z:Ω⟶ΩZ⊂IR2,such that for Z=z,ϑ and zr,θ is a linear transformation in r at the polar coordinates r,θ, with ϑr,θ=θand(2)zr,θ=rΓθ.which linksΩ to the unit disk with an easily invertible map for the domain delimited by Γc(Figure 1).Figure 1
Transforming map fromΩ to the unit disk ΩZ with ϑ=θ. The six homothetic curves have Radau’s nodes spacing.The applicationZΩ¯⟶Ω¯Z is a proper map that takes the compact region Ω¯=Ω∪∂Ω to another compact Ω¯Z=ΩZ∪∂ΩZ . Then, if dZ is bijective, the Hadamard–Caccioppoli theorem ensures that Z is a diffeomorphism. The local diffeomorphic map ZP in a point P∈Ω induces, by its differential, another natural linear map dZP:TPΩ⟶TZPΩZ between the corresponding tangent spaces TPΩ and TZPΩZ, provided by ordered bases and vectors ∂r,∂θ and ∂z,∂ϑ, respectively, as defined by the coordinate chart Z. Therefore, a function F, analytic in a Sobolev subspace, as described in (2), is diffeomorphically mapped by Z as Fr,θ⟶ψz,ϑ. Under these conditions, the notation ∂z/∂r≡∂rz:=zr and the differential bases vectors in the following form:∂r∂θ≡J∂z∂ϑ , with JT=∂z,ϑ/∂r,θ=∂rz∂θz∂rϑ∂θϑ⇒(3)Jr,θ=zr0zθ1.because∂rϑ=0 and ∂θϑ=1.Therefore, the first-order derivatives are(4)∂r≡zr∂z,∂θ≡∂ϑ+zθ∂z.Straightforward, one getszr:=∂z/∂r=1/Γ and zθ:=∂z/∂θ=−rΓθ/Γ2⇒zθ=−zΓθ/Γ.Noticing thatdetJ≠0 and dZP is a bijection, ZP is a local diffeomorphism.To determine the partial derivatives of the higher-order, it is sufficient to apply (4) successively.Considering the Helmholtz equationΔF=−κ2F and the Laplacian in polar coordinates multiplied by r2 to avoid divisions near zero, we obtain(5)ΔF⋅r2=z⋅Γ2∂2F∂r2+z⋅Γ⋅∂F∂r+∂2F∂θ2.By substituting the derivatives for the obtained expressions and grouping the terms, we obtain(6)ΔF⋅r2=A1⋅∂2ψ∂z2+A2⋅∂2ψ∂θ2+A3⋅∂2ψ∂z∂θ+A4⋅∂ψ∂z,where(7)A1=r2Γ2+zθ2=z2+zθ2.A2=1,A3=2zθ, A4=r/Γ+zθθ+zθ⋅zθz=z+zθθ+zθ⋅zθzThe right-hand side becomes−z2Γ2⋅κ2ψ.These are all transformations necessary to change the domain to a unit circle, considering the changes in the equation. With these equations, we write the system of equationsS1 with dimension n⋅2m+1×n+1⋅2m+1, as described in Section 4.


### 2.2. Normal Vector and Neumann Condition

This section presents the explicit expressions used to impose the Neumann condition, considering the diffeomorphic map. Starting with the original PDE, the gradient in polar coordinates is expressed as follows:∇⟶F=∂F/∂re^r+1/r∂F/∂θe^θ, where e^θ denotes the usual unit vectors. Considering the equations from section 2.1, we get∂F/∂r=1/Γ⋅∂ψ/∂z and ∂F/∂θ=∂ψ/∂θ+zθ⋅∂ψ/∂z.Asv⟶=rθe^r+re^θ is a tangent vector to the curve r=rθ described in polar coordinates, the unitary normal vector outside-oriented is n^=−rθ⋅e^θ+r⋅e^r/v⟶.Hence, the Neumann condition results in(8)∂F∂n^≔∇⟶F⋅n^=∂F/∂rr−rθ/r∂F/∂θv⟶.Using the expressions in 2.1, with a little algebra and observing that in∂Ωz=1 and r=Γ, we obtain(9)∂F∂n^=Γ2+Γθ2/Γ2∂ψ1,θ/∂z−Γθ/Γ∂ψ1,θ/∂θv⟶,where v⟶=Γ2+Γθ21/2. In this way, it is possible to build a system S2, with dimension 2m+1×2m+1⋅n+1, as described in Section 4, completing a system S1, and solving the eigenvalue problem.


## 2.1. Radial Reciprocal Map (RRM)

Let the open subsetΩ⊂ℝ2 be a simply connected domain with border ∂Ω=Γc described by a Jordan curve parameterized at the angle θr,θ=Γθ,θ, and we consider a diffeomorphism between open subsets of IR2, as given by the coordinate chart(1)Z:Ω⟶ΩZ⊂IR2,such that for Z=z,ϑ and zr,θ is a linear transformation in r at the polar coordinates r,θ, with ϑr,θ=θand(2)zr,θ=rΓθ.which linksΩ to the unit disk with an easily invertible map for the domain delimited by Γc(Figure 1).Figure 1
Transforming map fromΩ to the unit disk ΩZ with ϑ=θ. The six homothetic curves have Radau’s nodes spacing.The applicationZΩ¯⟶Ω¯Z is a proper map that takes the compact region Ω¯=Ω∪∂Ω to another compact Ω¯Z=ΩZ∪∂ΩZ . Then, if dZ is bijective, the Hadamard–Caccioppoli theorem ensures that Z is a diffeomorphism. The local diffeomorphic map ZP in a point P∈Ω induces, by its differential, another natural linear map dZP:TPΩ⟶TZPΩZ between the corresponding tangent spaces TPΩ and TZPΩZ, provided by ordered bases and vectors ∂r,∂θ and ∂z,∂ϑ, respectively, as defined by the coordinate chart Z. Therefore, a function F, analytic in a Sobolev subspace, as described in (2), is diffeomorphically mapped by Z as Fr,θ⟶ψz,ϑ. Under these conditions, the notation ∂z/∂r≡∂rz:=zr and the differential bases vectors in the following form:∂r∂θ≡J∂z∂ϑ , with JT=∂z,ϑ/∂r,θ=∂rz∂θz∂rϑ∂θϑ⇒(3)Jr,θ=zr0zθ1.because∂rϑ=0 and ∂θϑ=1.Therefore, the first-order derivatives are(4)∂r≡zr∂z,∂θ≡∂ϑ+zθ∂z.Straightforward, one getszr:=∂z/∂r=1/Γ and zθ:=∂z/∂θ=−rΓθ/Γ2⇒zθ=−zΓθ/Γ.Noticing thatdetJ≠0 and dZP is a bijection, ZP is a local diffeomorphism.To determine the partial derivatives of the higher-order, it is sufficient to apply (4) successively.Considering the Helmholtz equationΔF=−κ2F and the Laplacian in polar coordinates multiplied by r2 to avoid divisions near zero, we obtain(5)ΔF⋅r2=z⋅Γ2∂2F∂r2+z⋅Γ⋅∂F∂r+∂2F∂θ2.By substituting the derivatives for the obtained expressions and grouping the terms, we obtain(6)ΔF⋅r2=A1⋅∂2ψ∂z2+A2⋅∂2ψ∂θ2+A3⋅∂2ψ∂z∂θ+A4⋅∂ψ∂z,where(7)A1=r2Γ2+zθ2=z2+zθ2.A2=1,A3=2zθ, A4=r/Γ+zθθ+zθ⋅zθz=z+zθθ+zθ⋅zθzThe right-hand side becomes−z2Γ2⋅κ2ψ.These are all transformations necessary to change the domain to a unit circle, considering the changes in the equation. With these equations, we write the system of equationsS1 with dimension n⋅2m+1×n+1⋅2m+1, as described in Section 4.


## 2.2. Normal Vector and Neumann Condition

This section presents the explicit expressions used to impose the Neumann condition, considering the diffeomorphic map. Starting with the original PDE, the gradient in polar coordinates is expressed as follows:∇⟶F=∂F/∂re^r+1/r∂F/∂θe^θ, where e^θ denotes the usual unit vectors. Considering the equations from section 2.1, we get∂F/∂r=1/Γ⋅∂ψ/∂z and ∂F/∂θ=∂ψ/∂θ+zθ⋅∂ψ/∂z.Asv⟶=rθe^r+re^θ is a tangent vector to the curve r=rθ described in polar coordinates, the unitary normal vector outside-oriented is n^=−rθ⋅e^θ+r⋅e^r/v⟶.Hence, the Neumann condition results in(8)∂F∂n^≔∇⟶F⋅n^=∂F/∂rr−rθ/r∂F/∂θv⟶.Using the expressions in 2.1, with a little algebra and observing that in∂Ωz=1 and r=Γ, we obtain(9)∂F∂n^=Γ2+Γθ2/Γ2∂ψ1,θ/∂z−Γθ/Γ∂ψ1,θ/∂θv⟶,where v⟶=Γ2+Γθ21/2. In this way, it is possible to build a system S2, with dimension 2m+1×2m+1⋅n+1, as described in Section 4, completing a system S1, and solving the eigenvalue problem.


## 3. Bases and Operational Matrices

### 3.1. Chebyshev Polynomials as a Basis Set

The Chebyshev polynomialTnx is a polynomial in x of degree n, defined by the following relation:Tncosθ=cosnθ, θ∈0,π⊂IR, and x=cosθ⇒x∈−1,1⊂IR.Computationally, a quick way to get the polynomials up to ordern is to use the recurrence relationship as follows:Tk+1x=2xTkx−Tk−1x,k=1,2,…, where T0x=1 and T1x=x.The first six Chebyshev polynomials are(10)T0x=1T1x=xT2x=2x2−1T3x=4x3−3xT4x=8x4−8x2+1T5x=16x5−20x3+5xT6x=32x6−48x4+18x2−1.The relation of orthogonality with respect to weightwx=1−x2−1/2 is(11)∫−11TnxTmxdx1−x2=0:n≠mπ:n=m=0π2:n=m≠0..As a basis set, a continuous and bounded function can be approximated by(12)fx≃fnx=∑k=0nckTkx,limn⟶∞fx−fnx=0with geometric convergence for functions of classC∞ uniformly [18, 19].For continuous functions, the coefficients of the expansion are(13)ck=2π∫−11fx⋅Tkx1−x2dx,fork=1:n and half of this for k=0.Transforming to the physical spaceyk=∑j=0ncj⋅Tjxk, and considering Bk,j=Tjxk, we have the matrixial discretization Y=BC, where C=c0,c1,…,cnT, given by (9). For discrete approximation, it is useful to use the quadrature rules, as follows:The quadrature of Radau–Chebyshev for the range [–1,1) has the following weights and abscissas [19, 20]:(14)xj=−cos2πj2n+1,0≤j≤n,ω0=π2n+1,ωj=2π2n+1,1≤j≤n.For the interval (–1, 1], it is sufficient to consider opposite abscissas. In ascending order, reads as follows:(15)xj=cos2πn−j2n+1,0≤j≤n,ωN=π2n+1,ωj=2π2n+1,0≤j≤n−1.For a grid withn+1 points, this quadrature is exact until it reaches a polynomial of degree 2n [18]. In the expansion of Chebyshev, we can obtain the coefficients using the weights for (9), making the quadrature of Chebyshev–Radau by(16)ck=2π⋅∑j=0nωjTkxj⋅fxj=2π⋅∑j=0nα⋅π2n+1Tkxj⋅fxjck=22n+1⋅∑j=0nαjTkxj⋅fxj,where αn=1,αj=2j<n.ConsideringC=c0,c1,…,cnTY=fx0,fx1…,fxnT, the discrete transformation between physical space for spectral readsC=B−1Y, where Bj,k−1=2/2n+1⋅∑j=0nαjTkxj.The choice of nameB−1 becomes clear because B−1⋅B=B⋅B−1=I.With this transformation,u=x+1/2, the domain x∈−1,1⊂IR is mapped to u∈0,1⊂IR, the inverse transformation is x=2u−1, forming the complete polynomial basis known as Shifted Chebyshev [21]. This basis is discussed in detail in Bhrawy & Alofi [22], and it is applied in fractional equations [23]. In addition, this basis allows fast computation via FFT.In the frequency space, the matrix of differentiation for a shift Chebyshev read [24]DT=0206…0080…00012…⋮⋮⋮⋮⋱, such that, using Dirac notation, if the numerical approximation of the function is fnx=∑k=0nckTkx:=B|Cn, then, dfnx/dx=B|DT|Cn with the same procedure as for the second derivative, observing that d2fnx/dx2=B|DT2|Cn,


### 3.2. Fourier Series

For a periodic function, the natural choice of basis to numerical approximation is the most well-known discrete Fourier series [25].In complex form, consideringθ∈0,2π⊂IR, the discrete expansion of a periodic function is(17)fθ≃fmθ=∑k=−mmck⋅eikθ≔BF|CF2m+1.The coefficients can be obtained byck=1/2π∫02πfθ⋅e−ikθdθ, or using a discrete approach,(18)ck=1/2m+1∑j=02m+1fθj⋅e−ikθj.with−m≤k≤m⊂Z [18].Thus, the matrixBk,jinv=1/2m+1e−ik⋅θj, where θj=2π/2m+1⋅j denotes the inverse of BF, since Binv⋅BF=BF⋅Binv=I2m+1, an identity matrix.In the frequency space, using the matrix form, one obtains(19)dfmθdθ=B|DF|CF2m+1,DF=diagk.For the second derivative, we note thatDF2=DF2 and so forth to higher derivatives.Similar to the Chebyshev series, the Fourier series also presents geometrical convergence for bounded functions of classC∞ [13].


## 3.1. Chebyshev Polynomials as a Basis Set

The Chebyshev polynomialTnx is a polynomial in x of degree n, defined by the following relation:Tncosθ=cosnθ, θ∈0,π⊂IR, and x=cosθ⇒x∈−1,1⊂IR.Computationally, a quick way to get the polynomials up to ordern is to use the recurrence relationship as follows:Tk+1x=2xTkx−Tk−1x,k=1,2,…, where T0x=1 and T1x=x.The first six Chebyshev polynomials are(10)T0x=1T1x=xT2x=2x2−1T3x=4x3−3xT4x=8x4−8x2+1T5x=16x5−20x3+5xT6x=32x6−48x4+18x2−1.The relation of orthogonality with respect to weightwx=1−x2−1/2 is(11)∫−11TnxTmxdx1−x2=0:n≠mπ:n=m=0π2:n=m≠0..As a basis set, a continuous and bounded function can be approximated by(12)fx≃fnx=∑k=0nckTkx,limn⟶∞fx−fnx=0with geometric convergence for functions of classC∞ uniformly [18, 19].For continuous functions, the coefficients of the expansion are(13)ck=2π∫−11fx⋅Tkx1−x2dx,fork=1:n and half of this for k=0.Transforming to the physical spaceyk=∑j=0ncj⋅Tjxk, and considering Bk,j=Tjxk, we have the matrixial discretization Y=BC, where C=c0,c1,…,cnT, given by (9). For discrete approximation, it is useful to use the quadrature rules, as follows:The quadrature of Radau–Chebyshev for the range [–1,1) has the following weights and abscissas [19, 20]:(14)xj=−cos2πj2n+1,0≤j≤n,ω0=π2n+1,ωj=2π2n+1,1≤j≤n.For the interval (–1, 1], it is sufficient to consider opposite abscissas. In ascending order, reads as follows:(15)xj=cos2πn−j2n+1,0≤j≤n,ωN=π2n+1,ωj=2π2n+1,0≤j≤n−1.For a grid withn+1 points, this quadrature is exact until it reaches a polynomial of degree 2n [18]. In the expansion of Chebyshev, we can obtain the coefficients using the weights for (9), making the quadrature of Chebyshev–Radau by(16)ck=2π⋅∑j=0nωjTkxj⋅fxj=2π⋅∑j=0nα⋅π2n+1Tkxj⋅fxjck=22n+1⋅∑j=0nαjTkxj⋅fxj,where αn=1,αj=2j<n.ConsideringC=c0,c1,…,cnTY=fx0,fx1…,fxnT, the discrete transformation between physical space for spectral readsC=B−1Y, where Bj,k−1=2/2n+1⋅∑j=0nαjTkxj.The choice of nameB−1 becomes clear because B−1⋅B=B⋅B−1=I.With this transformation,u=x+1/2, the domain x∈−1,1⊂IR is mapped to u∈0,1⊂IR, the inverse transformation is x=2u−1, forming the complete polynomial basis known as Shifted Chebyshev [21]. This basis is discussed in detail in Bhrawy & Alofi [22], and it is applied in fractional equations [23]. In addition, this basis allows fast computation via FFT.In the frequency space, the matrix of differentiation for a shift Chebyshev read [24]DT=0206…0080…00012…⋮⋮⋮⋮⋱, such that, using Dirac notation, if the numerical approximation of the function is fnx=∑k=0nckTkx:=B|Cn, then, dfnx/dx=B|DT|Cn with the same procedure as for the second derivative, observing that d2fnx/dx2=B|DT2|Cn,


## 3.2. Fourier Series

For a periodic function, the natural choice of basis to numerical approximation is the most well-known discrete Fourier series [25].In complex form, consideringθ∈0,2π⊂IR, the discrete expansion of a periodic function is(17)fθ≃fmθ=∑k=−mmck⋅eikθ≔BF|CF2m+1.The coefficients can be obtained byck=1/2π∫02πfθ⋅e−ikθdθ, or using a discrete approach,(18)ck=1/2m+1∑j=02m+1fθj⋅e−ikθj.with−m≤k≤m⊂Z [18].Thus, the matrixBk,jinv=1/2m+1e−ik⋅θj, where θj=2π/2m+1⋅j denotes the inverse of BF, since Binv⋅BF=BF⋅Binv=I2m+1, an identity matrix.In the frequency space, using the matrix form, one obtains(19)dfmθdθ=B|DF|CF2m+1,DF=diagk.For the second derivative, we note thatDF2=DF2 and so forth to higher derivatives.Similar to the Chebyshev series, the Fourier series also presents geometrical convergence for bounded functions of classC∞ [13].


## 4. Solving the PDE in the Unit Circle

This section presents an overview of spectral methods to define the fundamental concepts used in this article. For polynomial bases, operational matrices constitute a similarity class [24]. Thus, the numerical approach can also be made in physical space, a pure pseudospectral method, using the MATLAB differentiating suite [26], with a slightly lower computational cost. To calculate the physical values of the functions at arbitrary points, barycentric interpolation is available for trigonometric [27] and polynomial bases [28]. Our choice of the spectral solution was the facility for the integration and evaluation of error estimates.The general form of a linear second-order two-dimensional PDE in polar coordinates for the unit circle is(20)A1∂2F∂z2+A2∂2F∂θ2+A3∂∂z∂F∂θ+A4∂F∂z+A5∂F∂θ+A6F=fr,θ.The domain isΩ¯≡Ω∪∂Ω, where Ω is the open ball z,θϵ0,1×0,2π⊂ℝ2 and ∂Ω is the circumference of the unit radius.Section3 shows the method to transform the original PDE to this domain.Observing thatψ:Ω¯⊂ℝ2 belongs to a Sobolev space [29] wherein a complete basis exists, a bivariate basis can approximate smooth functions with a small number of terms. A natural choice is a Fourier series for the azimuthal coordinates and a shift-Jacobi for the radial variable, usually Chebyshev or Legendre. Here, we choose the shift-Chebyshev.Hence, the bivariate basis is a Cartesian product of the shift-Chebyshev basis for the radial coordinate and Fourier basis for the azimuthal basis, both described in Section3, such that(21)Bz,θ=BT×BF.is the serial expansion of the numerical approximation, given by(22)ψ˜zp,θq=∑k=0n∑j=−mmmk,j⋅Tkzp⋅eijθq,p=0:n,q=0:2m,or in matrix formψ˜z,θ=Bz⋅M⋅BθFT, where M is the matrix of coefficients mk,j.To isolateM, the previous equation can be considered a Sylvester equation [30].Using the Kronecker product properties with operation vecM, which successively stacks the columns of a matrix, into one column [31], the last equation reads(23)vecψn+1,2m+1=B2m+1θ⊗Bn+1z⋅vecM.


### 4.1. Variable Coefficients

To deal with the variable coefficients, assumingU and V two generic matrices with the same dimension and U°V=R, where the symbol ° denotes the Hadamard product, then, using the “vec” operator, the following property is valid [32]:(24)vecU°V=vecU°vecV=vecR.BecauseC is the Hadamard product of the two-column vectors P and Q, we can write ci=pi⋅qi=∑j=1Npi,j⋅qi, where pi,j=0ifj≠i , corresponding to the ordinary matrix product C=diagP⋅Q.We consider matrixV as the discretization of ψ (or any of its derivatives) and matrix U as the discretization of the bivariable function by which the function is multiplied.For instance, the variable coefficientsA6z,θ in matrix form become A6,matrix=diagvecA6n+1,2m+1, and, finally, the term of (14) A6⋅ψ in matrix form is(25)A6,matrix⋅B2m+1θ⊗Bn+1z⋅vecM.For the derivatives, we proceed in the same manner. For instance, the conversion to the matrix form of the first term in (14) leads to A1∂2ψ/∂z2⟶A1,matrix⋅B2m+1θ⊗Bn+1z⋅Iθ⊗Dz2vecM, where Iθ is the identity of order 2m+1.


### 4.2. Internal Points Equation and Boundary Conditions

Letzn=z0,z1,…,znT the Radau–Chebyshev abscissas and θj=2π/2m+1j,withj=0:2m the azimuthal coordinates, the grid z0:n−1×θ0:2m determine the internal points of the domain. All left-hand side terms of (2) will be calculated in these points, resulting in a system S1 of n⋅2m+1×n+1⋅2m+1 equations, having in the right-hand side vecfzk,θj with n⋅2m+1 elements.In∂Ω, we have z=1, so to the Dirichlet condition, we write as(26)Bz1⋅M⋅B˜θ⟶Bθ⊗Bz1⋅vecM=02m+1,1.resulting in another system of equationsS2.For Neumann’s condition, the systemS2 is built with (8).Gathering the systemsS1 and S2, results in a square matrix, obeying(27)S1S2⋅vecM=−κ2B2m+1θ⊗Bn+1z0⋮⋱0⋅vecM,which is a generalized eigenvalue problem of typeS⋅v=−κ2⋅Q⋅v, where S and Q are square matrices with the same dimension, dense and semi-dense, and asymmetrical?Under the homogeneous condition of Dirichlet or pure Neumann’s condition∂F/∂n^=0, Q is the same and singular. Under Neumann’s pure condition, S is singular.A natural choice to solve this eigenvalue problem numerically is the QZ algorithm [33], which is a resident function of MATLAB. The computational cost is not small; for example, using n=30, which results in 2m+1=153, and using equation (28) to be derived in the next section, the dimensions of S and Q are 4743×4743. The run time was approximately 25 min on a standard laptop with 8 GB of RAM.The computational cost can be drastically reduced to less than 1 min if we restrict the search for eigenvalues to a subspace using the shift-and-invert technique [29, 34, 35], which requires an initial assumption about the eigenvalue range. Accuracy is slightly better because it allows for the use of a relatively large n, but sometimes eigenvalues and the corresponding eigenfunctions are skipped and do not appear in the result.This problem can be solved if we process it twice, with different ranges of initial settings, and then join the two sets.In a “tour de force,” using a preliminary result of the QZ algorithm with low order, one can build the chain of eigenvalues and eigenfunctions by varying the search range manually and making the set of eigenvalues/eigenfunctions.Again, let us note that the problematic pointx,y≡0,0 that would be mapped to z,θ≡0,anyθ does not belong to the grid, and the solution for this point will emerge naturally by the continuity condition implicitly in Radau’s nodes choice. Naturally, each eigenvector vec(M) must be reshaped to build matrix Mn+1,2m+1.It remains to be seen the manner in which we chose the order for expansion in thez coordinate, in which we should choose the expansion order for coordinate θ to have the same resolution. The following short theorem answers this question.


### 4.3. Theorem about Optimal Grid

As the range for azimuthal coordinates is approximately six times greater than the radial one, at first glance, it is expected that this would be the ratio between the number of points needed to have the same resolution in both coordinates, but the bases are different.Withn being the order of the shift-Chebyshev series to solve the PDE in the unit circle, the order of the Fourier series for the angular coordinate to have the same wave resolution is(28)2m+1=8sinπ/2n.Proof.
According to the Shannon–Nyquist theorem, [36] it requires 2 points per wavelength to solve an oscillating system with a trigonometric basis. For the Chebyshev basis in the standard domain, this density must be π points per wavelength [37]. For the shift-Chebyshev, half of the standard interval, this number becomes π/2.
The coordinates of the points in the expansion of shift-Chebyshev, in ascending order arexj=−cosjπ/n+1/2,j=0:n; hence, the distance between consecutive points is(29)Δxj=xj+1−xj=−cosj+1π/n+cosjπ/n2.
Using trigonometric identities(30)Δxj=sin2j+1π2nsinπ2n.
To maximize this distance, we consider the variational problemdΔxj/dj=0, which results in j=n−1/2, that has meaning for n odd, since j is an integer.
This is a well-known result and is expected because the maximum distance occurs between two central points.
Substitutingj into (30), we obtain Δxmax=sinπ/2n.
The maximum arc between two points at the azimuthal coordinate occurs at the periphery of the unit circle and is(31)Δℓ=2π2m+1.
Considering the conditions expressed in the first paragraph of this section, whereλ is the wavelength, we have λ/π/2=sinπ/2n for the polynomial basis and λ/2=2π/2m+1 for the trigonometric model. Therefore, π/2sinπ/2n=4π/2m+1.
This completes the proof.Corollary 1.
For a largen, the asymptotic behavior is sinπ/2n≃π/2n=8/2m+1, therefore m⟶8/πn as n⟶∞.Remark 1.
Occasionally, for boundary value problems, the expected solution is highly oscillatory in one variable, but not in the other, and different relations between the gridsz and j could save computational cost, but to ensure the same resolution in both variables, (30) remains mandatory.


## 4.1. Variable Coefficients

To deal with the variable coefficients, assumingU and V two generic matrices with the same dimension and U°V=R, where the symbol ° denotes the Hadamard product, then, using the “vec” operator, the following property is valid [32]:(24)vecU°V=vecU°vecV=vecR.BecauseC is the Hadamard product of the two-column vectors P and Q, we can write ci=pi⋅qi=∑j=1Npi,j⋅qi, where pi,j=0ifj≠i , corresponding to the ordinary matrix product C=diagP⋅Q.We consider matrixV as the discretization of ψ (or any of its derivatives) and matrix U as the discretization of the bivariable function by which the function is multiplied.For instance, the variable coefficientsA6z,θ in matrix form become A6,matrix=diagvecA6n+1,2m+1, and, finally, the term of (14) A6⋅ψ in matrix form is(25)A6,matrix⋅B2m+1θ⊗Bn+1z⋅vecM.For the derivatives, we proceed in the same manner. For instance, the conversion to the matrix form of the first term in (14) leads to A1∂2ψ/∂z2⟶A1,matrix⋅B2m+1θ⊗Bn+1z⋅Iθ⊗Dz2vecM, where Iθ is the identity of order 2m+1.


## 4.2. Internal Points Equation and Boundary Conditions

Letzn=z0,z1,…,znT the Radau–Chebyshev abscissas and θj=2π/2m+1j,withj=0:2m the azimuthal coordinates, the grid z0:n−1×θ0:2m determine the internal points of the domain. All left-hand side terms of (2) will be calculated in these points, resulting in a system S1 of n⋅2m+1×n+1⋅2m+1 equations, having in the right-hand side vecfzk,θj with n⋅2m+1 elements.In∂Ω, we have z=1, so to the Dirichlet condition, we write as(26)Bz1⋅M⋅B˜θ⟶Bθ⊗Bz1⋅vecM=02m+1,1.resulting in another system of equationsS2.For Neumann’s condition, the systemS2 is built with (8).Gathering the systemsS1 and S2, results in a square matrix, obeying(27)S1S2⋅vecM=−κ2B2m+1θ⊗Bn+1z0⋮⋱0⋅vecM,which is a generalized eigenvalue problem of typeS⋅v=−κ2⋅Q⋅v, where S and Q are square matrices with the same dimension, dense and semi-dense, and asymmetrical?Under the homogeneous condition of Dirichlet or pure Neumann’s condition∂F/∂n^=0, Q is the same and singular. Under Neumann’s pure condition, S is singular.A natural choice to solve this eigenvalue problem numerically is the QZ algorithm [33], which is a resident function of MATLAB. The computational cost is not small; for example, using n=30, which results in 2m+1=153, and using equation (28) to be derived in the next section, the dimensions of S and Q are 4743×4743. The run time was approximately 25 min on a standard laptop with 8 GB of RAM.The computational cost can be drastically reduced to less than 1 min if we restrict the search for eigenvalues to a subspace using the shift-and-invert technique [29, 34, 35], which requires an initial assumption about the eigenvalue range. Accuracy is slightly better because it allows for the use of a relatively large n, but sometimes eigenvalues and the corresponding eigenfunctions are skipped and do not appear in the result.This problem can be solved if we process it twice, with different ranges of initial settings, and then join the two sets.In a “tour de force,” using a preliminary result of the QZ algorithm with low order, one can build the chain of eigenvalues and eigenfunctions by varying the search range manually and making the set of eigenvalues/eigenfunctions.Again, let us note that the problematic pointx,y≡0,0 that would be mapped to z,θ≡0,anyθ does not belong to the grid, and the solution for this point will emerge naturally by the continuity condition implicitly in Radau’s nodes choice. Naturally, each eigenvector vec(M) must be reshaped to build matrix Mn+1,2m+1.It remains to be seen the manner in which we chose the order for expansion in thez coordinate, in which we should choose the expansion order for coordinate θ to have the same resolution. The following short theorem answers this question.


## 4.3. Theorem about Optimal Grid

As the range for azimuthal coordinates is approximately six times greater than the radial one, at first glance, it is expected that this would be the ratio between the number of points needed to have the same resolution in both coordinates, but the bases are different.Withn being the order of the shift-Chebyshev series to solve the PDE in the unit circle, the order of the Fourier series for the angular coordinate to have the same wave resolution is(28)2m+1=8sinπ/2n.Proof.
According to the Shannon–Nyquist theorem, [36] it requires 2 points per wavelength to solve an oscillating system with a trigonometric basis. For the Chebyshev basis in the standard domain, this density must be π points per wavelength [37]. For the shift-Chebyshev, half of the standard interval, this number becomes π/2.
The coordinates of the points in the expansion of shift-Chebyshev, in ascending order arexj=−cosjπ/n+1/2,j=0:n; hence, the distance between consecutive points is(29)Δxj=xj+1−xj=−cosj+1π/n+cosjπ/n2.
Using trigonometric identities(30)Δxj=sin2j+1π2nsinπ2n.
To maximize this distance, we consider the variational problemdΔxj/dj=0, which results in j=n−1/2, that has meaning for n odd, since j is an integer.
This is a well-known result and is expected because the maximum distance occurs between two central points.
Substitutingj into (30), we obtain Δxmax=sinπ/2n.
The maximum arc between two points at the azimuthal coordinate occurs at the periphery of the unit circle and is(31)Δℓ=2π2m+1.
Considering the conditions expressed in the first paragraph of this section, whereλ is the wavelength, we have λ/π/2=sinπ/2n for the polynomial basis and λ/2=2π/2m+1 for the trigonometric model. Therefore, π/2sinπ/2n=4π/2m+1.
This completes the proof.Corollary 1.
For a largen, the asymptotic behavior is sinπ/2n≃π/2n=8/2m+1, therefore m⟶8/πn as n⟶∞.Remark 1.
Occasionally, for boundary value problems, the expected solution is highly oscillatory in one variable, but not in the other, and different relations between the gridsz and j could save computational cost, but to ensure the same resolution in both variables, (30) remains mandatory.


## 5. Error Estimates

Let the internal product be defined as usualu|v:=∫Ωu⋅vdA, and(32)u≔u|u1/2≔∫Ωu2dA1/2.Assumingψ˜ is a numerical approximation of the eigenfunction of order n satisfying (7), we can obtain the residual, or defect, by(33)Δψ˜+κ˜2ψ˜=Resθ,z.According to Still [16], Theorem 1, being ε=Res/ψ˜ , there exists an eigenvalue κ2 such that κ2−κ˜2≤ε.In terms of relative errors,σκ=ε/κ∗2 , where κ∗2=minκ2,κ˜2.The analysis of the examples showed that this error bound is very conservative; therefore, this section shows an estimate of the errors considering the geometric convergence and other relations.


### 5.1. Eigenvalues Error Estimates via Rayleigh Quotient

In the Neumann or Dirichlet condition, given a coupleκ˜D2,ψ˜ candidate for the solution of the Helmholtz equation, a weak criterion for analyzing the error in the eigenvalue is to analyze the internal products as follows, considering the error concentrated in the eigenvalue.Sinceψ˜ is a candidate eigenfunction, it must correspond to an eigenvalue κ˜I2, when applied to the Helmholtz equation, such thatΔψ˜=−κ˜I2ψ˜ , so Δψ|ψ=−κ˜I2ψ|ψ. Using this integration approach, we obtain the value of κ˜I2 by the ratio −κ˜I2=Δψ|ψ/ψ|ψ, which is known as the Rayleigh quotient [38].Ifψ˜ is an analytical solution, we get κ˜I2=κ˜D2, of course, but for the numerical solutions the difference between the two values provides a way to get a standard deviation by(34)σκ=κ˜D2−κI2minκ˜D2,κ˜I2.The quadratures were obtained by padding zeros to matrixM to ensure the analysis of the function at many different points and accurate integration. Thus, we use at least n=36 for the quadrature process n=36.


### 5.2. Error Estimates for the Eigenfunction

To sharpen the error estimate, the value ofε, which is different from Still [16], was taken from the Rayleigh quotient in the previous section.As pointed out by Jones [8], as the approximation order is increased, the numerical eigenvalue oscillates around the exact value. This empirical process, through inspection in sequence, almost matches the Rayleigh quotient’s criteria. Defining dκ˜2=minκp2≠κ˜2κp2−κ˜2 [16], and ε=κ2⋅σκ⋅ψ|ψ1/2, the relative error estimate for the eigenfunction is(35)‖ψ−ψ˜‖‖ψ‖∞≤εdκ˜2‖ψ‖∞1+ε2dκ˜21/2.


### 5.3. Inequalities on Eigenvalues

A crisscross was made in the examples presented, and others agreed with some of the inequalities obeyed by the eigenvalues and eigenfunctions. For the condition of Dirichlet, some of them are as follows:(i)
Faber–Krahn’s inequality [38, 39] establishes that the first eigenvalue κ12 of the domain Ω with area A obeys the inequality(36)κ12⋅Aπ≥J0,12,whereJ0,1 is the first root of the Bessel function of order 0 andJ0,12=5.783185962946773 corresponds to the first eigenvalue of the unit circle, therefore(37)CFK=J0,12⋅πκ12⋅A≤1.(ii)
The eigenfunctionψ1 associated with first eigenvalue κ1 does not change its signal in Ω and obeys CPayne=4π∫Ωψ12dA/κ12∫Ωψ1dA2≤1, with the equality held only in the circular domain [40].(iii)
Polya et al. [41] derived that the ratio between two consecutive eigenvalues must obey κn+12/κn2≤3, where it was conjectured that κ22/κ12 attains its maximum for the circular regions i.e.,κ22/κ12≤J1,1/J0,12=2.5387…. To follow the same form as the previous two constants, we write CPolya=κ12/κ22J1,1/J0,12≤1 and the equality is valid only for the unit circle. The conjecture became a theorem, with the proof given in [42].Other inequalities, equality, and properties can be found in [43], which provides a detailed review of the subject’s history.


#### 5.3.1. Asymptotic Behavior

The asymptotic behavior of the eigenvalues in the Dirichlet condition is conditioned by Weyl Law and enhanced by Baltes [44]as follows:(38)κn2∼4πnA+4πnLA−3/2,where n is the order of the eigenvalue counting the duplicates, L is the perimeter, and A is the flat domain area delimited by Γθ.The Weyl law can be used as a guide to inspect the limit when the numerical process is no longer defined to resolve high harmonics.


## 5.1. Eigenvalues Error Estimates via Rayleigh Quotient

In the Neumann or Dirichlet condition, given a coupleκ˜D2,ψ˜ candidate for the solution of the Helmholtz equation, a weak criterion for analyzing the error in the eigenvalue is to analyze the internal products as follows, considering the error concentrated in the eigenvalue.Sinceψ˜ is a candidate eigenfunction, it must correspond to an eigenvalue κ˜I2, when applied to the Helmholtz equation, such thatΔψ˜=−κ˜I2ψ˜ , so Δψ|ψ=−κ˜I2ψ|ψ. Using this integration approach, we obtain the value of κ˜I2 by the ratio −κ˜I2=Δψ|ψ/ψ|ψ, which is known as the Rayleigh quotient [38].Ifψ˜ is an analytical solution, we get κ˜I2=κ˜D2, of course, but for the numerical solutions the difference between the two values provides a way to get a standard deviation by(34)σκ=κ˜D2−κI2minκ˜D2,κ˜I2.The quadratures were obtained by padding zeros to matrixM to ensure the analysis of the function at many different points and accurate integration. Thus, we use at least n=36 for the quadrature process n=36.


## 5.2. Error Estimates for the Eigenfunction

To sharpen the error estimate, the value ofε, which is different from Still [16], was taken from the Rayleigh quotient in the previous section.As pointed out by Jones [8], as the approximation order is increased, the numerical eigenvalue oscillates around the exact value. This empirical process, through inspection in sequence, almost matches the Rayleigh quotient’s criteria. Defining dκ˜2=minκp2≠κ˜2κp2−κ˜2 [16], and ε=κ2⋅σκ⋅ψ|ψ1/2, the relative error estimate for the eigenfunction is(35)‖ψ−ψ˜‖‖ψ‖∞≤εdκ˜2‖ψ‖∞1+ε2dκ˜21/2.


## 5.3. Inequalities on Eigenvalues

A crisscross was made in the examples presented, and others agreed with some of the inequalities obeyed by the eigenvalues and eigenfunctions. For the condition of Dirichlet, some of them are as follows:(i)
Faber–Krahn’s inequality [38, 39] establishes that the first eigenvalue κ12 of the domain Ω with area A obeys the inequality(36)κ12⋅Aπ≥J0,12,whereJ0,1 is the first root of the Bessel function of order 0 andJ0,12=5.783185962946773 corresponds to the first eigenvalue of the unit circle, therefore(37)CFK=J0,12⋅πκ12⋅A≤1.(ii)
The eigenfunctionψ1 associated with first eigenvalue κ1 does not change its signal in Ω and obeys CPayne=4π∫Ωψ12dA/κ12∫Ωψ1dA2≤1, with the equality held only in the circular domain [40].(iii)
Polya et al. [41] derived that the ratio between two consecutive eigenvalues must obey κn+12/κn2≤3, where it was conjectured that κ22/κ12 attains its maximum for the circular regions i.e.,κ22/κ12≤J1,1/J0,12=2.5387…. To follow the same form as the previous two constants, we write CPolya=κ12/κ22J1,1/J0,12≤1 and the equality is valid only for the unit circle. The conjecture became a theorem, with the proof given in [42].Other inequalities, equality, and properties can be found in [43], which provides a detailed review of the subject’s history.


### 5.3.1. Asymptotic Behavior

The asymptotic behavior of the eigenvalues in the Dirichlet condition is conditioned by Weyl Law and enhanced by Baltes [44]as follows:(38)κn2∼4πnA+4πnLA−3/2,where n is the order of the eigenvalue counting the duplicates, L is the perimeter, and A is the flat domain area delimited by Γθ.The Weyl law can be used as a guide to inspect the limit when the numerical process is no longer defined to resolve high harmonics.


## 5.3.1. Asymptotic Behavior

The asymptotic behavior of the eigenvalues in the Dirichlet condition is conditioned by Weyl Law and enhanced by Baltes [44]as follows:(38)κn2∼4πnA+4πnLA−3/2,where n is the order of the eigenvalue counting the duplicates, L is the perimeter, and A is the flat domain area delimited by Γθ.The Weyl law can be used as a guide to inspect the limit when the numerical process is no longer defined to resolve high harmonics.


## 6. Examples for Helmholtz Problems

This section presents problems related to the eigenvalues of the Helmholtz equation, as it is favorable for comparison with literature results.However, for the presented reasoning, the proposed diffeomorphism can be used in other eigenvalue problems such as Hamiltonians. In the examples, the first 2 graphs will show the following:(i)
some modes of vibration, seen by nodal lines;(ii)
one of the eigenfunctionsψk.A third chart displays the following error estimates for the obtained results as follows:σκ: for eigenvalues, according to Rayleigh’s quotient, compared to that obtained by the matrix equation; σψ: error of eigenfunctions according (35).This error analysis shows the spectral accuracy for hundreds of eigenvalues and eigenfunctions obtained in one shot.


### 6.1. Elliptical Domain-Dirichlet

In 2007, Wilson and Scharstein [17] published a detailed study of eigenvalues for elliptical domains, comparing the results of the Galerkin method with the method they proposed using Mathieu functions. The eigenvalues were presented with nine digits for the Dirichlet and Neumann conditions.With the precision developed here, there are at least ten digits of precision in the first 50 eigenvalues, and we identify some more undetected there.The ellipse that delimits the domain is described by the equation(39)Γθ=13sin2θ+1,θ∈0,2π⊂IR.Regarding the vibration modes, among the first 40 published by Wilson and Scharstein, we identified five more that did not emerge from the numerical process using Mathieu’s functions in Dirichlet’s homogeneous condition. Figure2 illustrates the vibration modes from 26 to 50, including those not detected by the Mathieu functions’ based algorithm and Figure 3 shows the vibration mode 37 and Figure 4 the errors for eigenvalues and eigenfunctions.Figure 2
Vibration modes from26th to 50th in ascending order of eigenvalues. The modes 36, 37, 40, 41, and 43 were not detected in [17] so mode 45 appears as 40th there.Figure 3
37th eigenfunction with κ372=3.4121083712669e+02.Figure 4
Error estimates for eigenvalues and eigenfunctions.


### 6.2. Elliptical Domain-Neumann

For the Dirichlet condition, let us now compare the results considering the Neumann boundary condition:∇⟶F⋅n^=0on∂Ω for the same elliptic curve.Figure5 illustrates the vibration modes from 26 to 50, Figure 6 shows the vibration mode 37, and Figure 7 shows the errors for eigenvalues and eigenfunctions.Figure 5
Vibration modes from the26th to 50th in ascending order of eigenvalues. The modes 36, 39, 40, and 43 were not detected in [17] so mode 44 appears as 40th there.Figure 6
36th eigenfunction κ362=2.3316072649495e+02 in Neumann condition.Figure 7
Error estimates for eigenvalues and eigenfunctions.Although some vibration modes have escaped the algorithm’s perception with Mathieu’s functions, all eigenvalues presented are in complete agreement with those we obtained, with a slight difference only in the last digit in two of the eigenvalues shown.


### 6.3. Trefoil-Dirichlet

This domain is delimited by the Jordan curve as follows:(40)∂Ω=Γθ≔1+0.3⋅cos3θ+0.2cosθ,θ∈0,2π⊂R.Barnet & Hessel [45] carried out a deep study of high eigenvalues in this case, using a sophisticated transformation map with fundamental solutions and dividing the space to search for solutions and obtain spectral precision. As mentioned in the introduction, RRM diffeomorphism can obtain hundreds of eigenvalues with spectral accuracy in one shot with a low computational cost. Moreover, it is also possible to choose the range using the shift and invert technique [35] for the matrix equation to improve the performance run time by almost 30 times, although it is restricted to grid resolution, as discussed in Section 4.1. With the flexible finite element method, the error in the first 20 eigenvalues is less than 0.5%.Compared with the results of Barnet and Hessel [45], who kindly made available the software MPSpack for free download, there is at least ten-digit agreement until the 200th harmonic.Figure8 illustrates the vibration modes from 76 to 100.Figure 8
Vibrations modes from 76 to 100.Figure9 shows the vibration mode 93 and Figure 10 shows the errors for eigenvalues and eigenfunctions.Figure 9
93rd eigenfunction in Dirichlet condition, whose eigenvalue is κ93=19.94995891589. All digits match with [39].Figure 10
Error estimates for eigenvalues and eigenfunctions.


### 6.4. Limacon-Dirichlet

The limaconΓθ=2.5+cosθ+2sinθ,withθ∈0,2π⊂R, was studied by Atkinson [46] using “ridge polynomials” and other bases in a Galerkin approach. This article shows the first two eigenvalues with five correct digits. As in other examples, we obtain at least 10 digits for the first hundred harmonics. A similar Pascal limacon was used by Kuttler and Sigillito [47] to model the outline of the human mitral valve for noninvasive diagnostics and waveguides, studying its vibration modes.Figure11 illustrates the vibration modes from 73 to 100, Figure 12 shows the vibration mode 93, and Figure 13 shows the errors for eigenvalues and eigenfunctions.Figure 11
Vibrations modes from 76 to 100.Figure 12
93rd eigenfunction in Dirichlet condition.Figure 13
Error estimates for eigenvalues and eigenfunctions.


### 6.5. Pentafoil-Neumann Condition

This problem was also addressed in [45] under the homogeneous condition of Dirichlet for very high frequencies. Here, we analyzed Neumann’s pure condition, compared it with the finite element method to validate the first four digits, and found a difference of less than 0.5% in the first 32 eigenvalues.For this example, the Jordan curve is(41)∂Ω=Γθ=1+0.3⋅sin5t,θ∈0,2π⊂R.Figure14 illustrates the vibration modes from 53 to 100, Figure 15 shows the vibration mode 80, and Figure 16 shows the errors for eigenvalues and eigenfunctions.Figure 14
Vibrations modes from 76 to 100.Figure 15
80th eigenfunction, counting duplicates, in Neumann condition, whose eigenvalue is κ80=15.895457.Figure 16
Error estimates for eigenvalues and eigenfunctions.


### 6.6. Epitrochoid-Dirichlet Condition

This curve was studied in a boundary value problem by Liu [48] under the Dirichlet condition. Here, we search for eigenvalues for the domain enlaced by this curve, whose form in polar coordinates is Γθ:=26−10⋅cos4θ, where θ∈0,2π⊂R.To the best of our knowledge, there have been no published results on these eigenvalues. Therefore, we used the FE method to validate and compare the results. For the first 50 eigenvalues, the difference is less than 0.5%.Figure17 illustrates the vibration modes from 1 to 25, Figure 18 shows the vibration mode 20, and Figure 19 shows the errors for eigenvalues and eigenfunctions.Figure 17
Vibrations modes from 1 to 25.Figure 18
20th eigenfunction, counting duplicates, in Dirichlet condition, whose eigenvalue is κ50=1.969489165159.Figure 19
Error estimates for eigenvalues and eigenfunctions.Using the shift-invert technique [29, 34, 35] and limiting the search of eigenvalues to short-range is feasible to reach the machine limit in terms of precision, increasing the expansion series until N=40, as we did in this example, for the range between the 129 and 249 eigenvalue orders with less than 30 s runtime (see Figure 20).Figure 20
Error estimates for high eigenvalues limiting the range of search.Although complicated, this applies to all the other examples. It is worth mentioning that among these 120 eigenvalues, in 16, the relative error estimate was less thaneps=2,204−16, the resolution limit of MATLAB, including four with zero error, for which we assigned the proper eps to appear in the semilog chart.


## 6.1. Elliptical Domain-Dirichlet

In 2007, Wilson and Scharstein [17] published a detailed study of eigenvalues for elliptical domains, comparing the results of the Galerkin method with the method they proposed using Mathieu functions. The eigenvalues were presented with nine digits for the Dirichlet and Neumann conditions.With the precision developed here, there are at least ten digits of precision in the first 50 eigenvalues, and we identify some more undetected there.The ellipse that delimits the domain is described by the equation(39)Γθ=13sin2θ+1,θ∈0,2π⊂IR.Regarding the vibration modes, among the first 40 published by Wilson and Scharstein, we identified five more that did not emerge from the numerical process using Mathieu’s functions in Dirichlet’s homogeneous condition. Figure2 illustrates the vibration modes from 26 to 50, including those not detected by the Mathieu functions’ based algorithm and Figure 3 shows the vibration mode 37 and Figure 4 the errors for eigenvalues and eigenfunctions.Figure 2
Vibration modes from26th to 50th in ascending order of eigenvalues. The modes 36, 37, 40, 41, and 43 were not detected in [17] so mode 45 appears as 40th there.Figure 3
37th eigenfunction with κ372=3.4121083712669e+02.Figure 4
Error estimates for eigenvalues and eigenfunctions.


## 6.2. Elliptical Domain-Neumann

For the Dirichlet condition, let us now compare the results considering the Neumann boundary condition:∇⟶F⋅n^=0on∂Ω for the same elliptic curve.Figure5 illustrates the vibration modes from 26 to 50, Figure 6 shows the vibration mode 37, and Figure 7 shows the errors for eigenvalues and eigenfunctions.Figure 5
Vibration modes from the26th to 50th in ascending order of eigenvalues. The modes 36, 39, 40, and 43 were not detected in [17] so mode 44 appears as 40th there.Figure 6
36th eigenfunction κ362=2.3316072649495e+02 in Neumann condition.Figure 7
Error estimates for eigenvalues and eigenfunctions.Although some vibration modes have escaped the algorithm’s perception with Mathieu’s functions, all eigenvalues presented are in complete agreement with those we obtained, with a slight difference only in the last digit in two of the eigenvalues shown.


## 6.3. Trefoil-Dirichlet

This domain is delimited by the Jordan curve as follows:(40)∂Ω=Γθ≔1+0.3⋅cos3θ+0.2cosθ,θ∈0,2π⊂R.Barnet & Hessel [45] carried out a deep study of high eigenvalues in this case, using a sophisticated transformation map with fundamental solutions and dividing the space to search for solutions and obtain spectral precision. As mentioned in the introduction, RRM diffeomorphism can obtain hundreds of eigenvalues with spectral accuracy in one shot with a low computational cost. Moreover, it is also possible to choose the range using the shift and invert technique [35] for the matrix equation to improve the performance run time by almost 30 times, although it is restricted to grid resolution, as discussed in Section 4.1. With the flexible finite element method, the error in the first 20 eigenvalues is less than 0.5%.Compared with the results of Barnet and Hessel [45], who kindly made available the software MPSpack for free download, there is at least ten-digit agreement until the 200th harmonic.Figure8 illustrates the vibration modes from 76 to 100.Figure 8
Vibrations modes from 76 to 100.Figure9 shows the vibration mode 93 and Figure 10 shows the errors for eigenvalues and eigenfunctions.Figure 9
93rd eigenfunction in Dirichlet condition, whose eigenvalue is κ93=19.94995891589. All digits match with [39].Figure 10
Error estimates for eigenvalues and eigenfunctions.


## 6.4. Limacon-Dirichlet

The limaconΓθ=2.5+cosθ+2sinθ,withθ∈0,2π⊂R, was studied by Atkinson [46] using “ridge polynomials” and other bases in a Galerkin approach. This article shows the first two eigenvalues with five correct digits. As in other examples, we obtain at least 10 digits for the first hundred harmonics. A similar Pascal limacon was used by Kuttler and Sigillito [47] to model the outline of the human mitral valve for noninvasive diagnostics and waveguides, studying its vibration modes.Figure11 illustrates the vibration modes from 73 to 100, Figure 12 shows the vibration mode 93, and Figure 13 shows the errors for eigenvalues and eigenfunctions.Figure 11
Vibrations modes from 76 to 100.Figure 12
93rd eigenfunction in Dirichlet condition.Figure 13
Error estimates for eigenvalues and eigenfunctions.


## 6.5. Pentafoil-Neumann Condition

This problem was also addressed in [45] under the homogeneous condition of Dirichlet for very high frequencies. Here, we analyzed Neumann’s pure condition, compared it with the finite element method to validate the first four digits, and found a difference of less than 0.5% in the first 32 eigenvalues.For this example, the Jordan curve is(41)∂Ω=Γθ=1+0.3⋅sin5t,θ∈0,2π⊂R.Figure14 illustrates the vibration modes from 53 to 100, Figure 15 shows the vibration mode 80, and Figure 16 shows the errors for eigenvalues and eigenfunctions.Figure 14
Vibrations modes from 76 to 100.Figure 15
80th eigenfunction, counting duplicates, in Neumann condition, whose eigenvalue is κ80=15.895457.Figure 16
Error estimates for eigenvalues and eigenfunctions.


## 6.6. Epitrochoid-Dirichlet Condition

This curve was studied in a boundary value problem by Liu [48] under the Dirichlet condition. Here, we search for eigenvalues for the domain enlaced by this curve, whose form in polar coordinates is Γθ:=26−10⋅cos4θ, where θ∈0,2π⊂R.To the best of our knowledge, there have been no published results on these eigenvalues. Therefore, we used the FE method to validate and compare the results. For the first 50 eigenvalues, the difference is less than 0.5%.Figure17 illustrates the vibration modes from 1 to 25, Figure 18 shows the vibration mode 20, and Figure 19 shows the errors for eigenvalues and eigenfunctions.Figure 17
Vibrations modes from 1 to 25.Figure 18
20th eigenfunction, counting duplicates, in Dirichlet condition, whose eigenvalue is κ50=1.969489165159.Figure 19
Error estimates for eigenvalues and eigenfunctions.Using the shift-invert technique [29, 34, 35] and limiting the search of eigenvalues to short-range is feasible to reach the machine limit in terms of precision, increasing the expansion series until N=40, as we did in this example, for the range between the 129 and 249 eigenvalue orders with less than 30 s runtime (see Figure 20).Figure 20
Error estimates for high eigenvalues limiting the range of search.Although complicated, this applies to all the other examples. It is worth mentioning that among these 120 eigenvalues, in 16, the relative error estimate was less thaneps=2,204−16, the resolution limit of MATLAB, including four with zero error, for which we assigned the proper eps to appear in the semilog chart.


## 7. Asymptotic Behavior-Dirichlet Condition

The inequalities referred to in Section5.3, with the three constants CFK, CPayne, and CPolya, as well as the asymptotic behavior referred to in (24), are necessary conditions to be satisfied by the solutions. These three constants must be less than unity. The divergence of the asymptotic behavior reveals the extent to which the numerical method has sufficient grid resolution to provide reliable results. All verifications were carried out with the grid n+1×2m+1≡33×163, which corresponds to 5379 degrees of freedom. Let kappan,A2 be κn,A2 the asymptotic order N eigenvalue, given by (38), and κ˜n2 be the numerical result for this eigenvalue, and we define ξ:=κn,A2−κ˜n2/κn,A2. As κn2⟶κn,A2 for n⟶∞, as long as the numerical approach is good, ξ⟶0, for large N, with oscillatory behavior as N increases. Due to the limitation of points by the wavelength of the grid and floating-point [4], the eigenvalues obtained numerically will be reliable only up to a specific order. The behavior of ξ with respect to the harmonic order can reveal the reliability limit. Graphs of the behavior and constants related to some theoretical inequalities are presented hereafter.


### 7.1. Ellipse from Example 6.1.1

For this domain, we gotCFK=0.8107, CPayne=0.8420, and CPolya=0.6930, as shown in Figure 21.Figure 21
The behavior ofξ with respect to the numerical eigenvalue order. The ξ divergence starts around N  = 600.


### 7.2. Trefoil from Example 6.3

For this domain we gotCFK=0.75887, CPayne=0.8305, and CPolya=0.7807, as shown in Figure 22.Figure 22
The behavior ofξ with respect to the numerical eigenvalue order.


### 7.3. Limacon from Example 6.4

For this domain, the three constants areCFK=0.9945, CPayne=0.9957, and CPolya=0.9537, as shown in Figure 23.Figure 23
The behavior ofξ with respect to the numerical eigenvalue order.The proximity of the units in these constants shows that the curve is very similar to that of the circle.


### 7.4. Pentafoil-6.4-Dirchlet Condition

For this domain, the inequality constants areCFK=0.6570, CPayne=0.8472, and CPolya=0.9365, as shown in Figure 24.Figure 24
The behavior ofξ with respect to the numerical eigenvalue order.


### 7.5. Epitrochoid-6.6-Dirchlet Condition

For this domain, the inequality constants areCFK=0.8444, CPayne=0.9154, and CPolya=0.9511, as shown in Figure 25.Figure 25
The behavior ofξ with respect to the numerical eigenvalue order.


## 7.1. Ellipse from Example 6.1.1

For this domain, we gotCFK=0.8107, CPayne=0.8420, and CPolya=0.6930, as shown in Figure 21.Figure 21
The behavior ofξ with respect to the numerical eigenvalue order. The ξ divergence starts around N  = 600.


## 7.2. Trefoil from Example 6.3

For this domain we gotCFK=0.75887, CPayne=0.8305, and CPolya=0.7807, as shown in Figure 22.Figure 22
The behavior ofξ with respect to the numerical eigenvalue order.


## 7.3. Limacon from Example 6.4

For this domain, the three constants areCFK=0.9945, CPayne=0.9957, and CPolya=0.9537, as shown in Figure 23.Figure 23
The behavior ofξ with respect to the numerical eigenvalue order.The proximity of the units in these constants shows that the curve is very similar to that of the circle.


## 7.4. Pentafoil-6.4-Dirchlet Condition

For this domain, the inequality constants areCFK=0.6570, CPayne=0.8472, and CPolya=0.9365, as shown in Figure 24.Figure 24
The behavior ofξ with respect to the numerical eigenvalue order.


## 7.5. Epitrochoid-6.6-Dirchlet Condition

For this domain, the inequality constants areCFK=0.8444, CPayne=0.9154, and CPolya=0.9511, as shown in Figure 25.Figure 25
The behavior ofξ with respect to the numerical eigenvalue order.


## 8. Conclusion

In this study, a diffeomorphic RRM is easily invertible, with all transformations resulting in PDE equations in two-dimensional planar domains.The application of the RRM combined with the choice of Radau nodes allows the construction of complete and orthogonal bivariate bases, leading to a well-conditioned linear system and the resolution of eigenvalue problems. In one shot, we obtain hundreds of eigenvalues and eigenfunctions with spectral precision and low computational cost, either in the Dirichlet or Neumann condition, with no one left behind until the resolution limit.There is no need fora priori symmetry considerations because they emerge spontaneously in the numerical solution process once the optimal grid presented here is chosen, which imposes the same wavelength resolution on a bivariable basis. The error estimates in the eigenvalues and the eigenfunctions were also discussed.Finally, using the theoretical expression of asymptotic behavior, we were able to show that the spectral accuracy extends to harmonics of the order of 4 hundred or more for smooth Jordan curves.In future work, we will analyze the application of a similar transformation for three-dimensional problems and the adaptation for cases in which the equation that describes the Jordan curve does not have continuous derivatives.



---

*Source: 7379766-2022-04-11.xml*