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

,is 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. ,e domain must be enclosed by a Jordan curve that is univocally described in polar coordinates by a C2 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. ,e 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.


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. is package was used to validate the results. For the Helmholtz equation in particular, even in nontensor domains, some researchers [6][7][8][9][10][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. e 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 C 2 ) because we restricted the study to 2 nd order PDEs.
Since the derivative coefficients are considered to be variable, the method applies to several classes of PDEs (elliptical, hyperbolic, or parabolic). Section 2.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 Section 3, 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. e resolution of 2 nd 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. e asymptotical behavior is commented upon for the Dirichlet condition, in addition to some inequalities for the eigenvalues.
Finally, in section 6, 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.

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

Radial Reciprocal Map (RRM).
Let the open subset Ω ⊂ R 2 be a simply connected domain with border zΩ � Γ c described by a Jordan curve parameterized at the angle θ(r, θ) � (Γ(θ), θ), and we consider a diffeomorphism between open subsets of IR 2 , as given by the coordinate chart such that for Z � (z, ϑ) and z(r, θ) is a linear transformation in r at the polar coordinates (r, θ), with ϑ(r, θ) � θ and which links Ω to the unit disk with an easily invertible map for the domain delimited by Γ c (Figure 1). e application Z(Ω) ⟶ Ω Z is a proper map that takes the compact region Ω � Ω ∪ zΩ to another compact en, if dZ is bijective, the Hadamard-Caccioppoli theorem ensures that Z is a diffeomorphism. e local diffeomorphic map Z P in a point P ∈ Ω induces, by its differential, another natural linear map dZ P : T P Ω ⟶ T Z(P) Ω Z between the corresponding tangent spaces T P Ω and T Z(P) Ω Z , provided by ordered bases and vectors (z r , z θ ) and (z z , z ϑ ), respectively, as defined by the coordinate chart Z. erefore, a function F, analytic in a Sobolev subspace, as described in (2), is diffeomorphically mapped by Z as F(r, θ) ⟶ ψ(z, ϑ). Under these condiequation. With these equations, we write the system of equations S 1 with dimension n · (2m + 1) × (n + 1) · (2m + 1), as described in Section 4.

Normal Vector and Neumann Condition.
is 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: where e θ denotes the usual unit vectors. Considering the equations from section 2.1, we get zF/zr � 1/Γ · zψ/zz and zF/zθ � zψ/zθ + z θ · zψ/zz.
As v → � 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 Using the expressions in 2.1, with a little algebra and observing that in zΩz � 1 and r � Γ, we obtain where ‖ v → ‖ � (Γ 2 + Γ 2 θ ) 1/2 . In this way, it is possible to build a system S 2 , with dimension (2m + 1) × [(2m + 1) · (n + 1)], as described in Section 4, completing a system S 1 , and solving the eigenvalue problem.

Chebyshev Polynomials as a Basis Set.
e Chebyshev polynomial T n (x) is a polynomial in x of degree n, defined by the following relation: Computationally, a quick way to get the polynomials up to order n is to use the recurrence relationship as follows: e first six Chebyshev polynomials are e relation of orthogonality with respect to weight As a basis set, a continuous and bounded function can be approximated by with geometric convergence for functions of class C ∞ uniformly [18,19].
For continuous functions, the coefficients of the expansion are for k � 1: n and half of this for k � 0. Transforming to the physical space y k � n j�0 c j · T j (x k ), and considering B k,j � T j (x k ), we have the matrixial discretization Y � BC, where C � [c 0 , c 1 , . . . , c n ] T , given by (9). For discrete approximation, it is useful to use the quadrature rules, as follows: e quadrature of Radau-Chebyshev for the range [-1,1) has the following weights and abscissas [19,20]: , For the interval (-1, 1], it is sufficient to consider opposite abscissas. In ascending order, reads as follows: Figure 1: Transforming map from Ω to the unit disk Ω Z with ϑ � θ.
e six homothetic curves have Radau's nodes spacing.
x j � cos , For a grid with n + 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 where the discrete transformation between physical space for spectral reads forming the complete polynomial basis known as Shifted Chebyshev [21]. is 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]

Fourier Series.
For a periodic function, the natural choice of basis to numerical approximation is the most wellknown discrete Fourier series [25].
In complex form, considering θ ∈ [0, 2π) ⊂ IR, the discrete expansion of a periodic function is e coefficients can be obtained by c k � 1/2π 2π 0 f(θ) · e −ikθ dθ, or using a discrete approach, In the frequency space, using the matrix form, one obtains For the second derivative, we note that D 2 F � (D F ) 2 and so forth to higher derivatives.
Similar to the Chebyshev series, the Fourier series also presents geometrical convergence for bounded functions of class C ∞ [13].

Solving the PDE in the Unit Circle
is 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]. us, 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. e general form of a linear second-order two-dimensional PDE in polar coordinates for the unit circle is and zΩ is the circumference of the unit radius.
Section 3 shows the method to transform the original PDE to this domain.
Observing that ψ: Ω ⊂ R 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 Section 3, such that is the serial expansion of the numerical approximation, given by To isolate M, the previous equation can be considered a Sylvester equation [30].
Using the Kronecker product properties with operation vec (M), which successively stacks the columns of a matrix, into one column [31], the last equation reads

Variable Coefficients.
To deal with the variable coefficients, assuming U 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]: Because C is the Hadamard product of the two-column vectors P and Q, we can write We consider matrix V 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 coefficients A 6 (z, θ) in matrix form become A 6,matrix � diag(vec(A 6(n+1,2m+1) )), and, finally, the term of (14) A 6 · ψ in matrix form is 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 where I θ is the identity of order 2m + 1.
In zΩ, we have z � 1, so to the Dirichlet condition, we write as resulting in another system of equations S 2 . For Neumann's condition, the system S 2 is built with (8). Gathering the systems S 1 and S 2 , results in a square matrix, obeying which is a generalized eigenvalue problem of type S · 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 (zF/zn � 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. e 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. e run time was approximately 25 min on a standard laptop with 8 GB of RAM.
e 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.
is 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 point (x, 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 M n+1,2m+1 .
It remains to be seen the manner in which we chose the order for expansion in the z coordinate, in which we should choose the expansion order for coordinate θ to have the same resolution. e following short theorem answers this question.

eorem 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 Mathematical Problems in Engineering between the number of points needed to have the same resolution in both coordinates, but the bases are different.
With n 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 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. e coordinates of the points in the expansion of shift-Chebyshev, in ascending order are x j � [−cos(jπ/n) + 1]/2, j � 0: n; hence, the distance between consecutive points is Using trigonometric identities To maximize this distance, we consider the variational problem d(Δx j )/dj � 0, which results in j � (n − 1)/2, that has meaning for n odd, since j is an integer.
is is a well-known result and is expected because the maximum distance occurs between two central points.

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 grids z and j could save computational cost, but to ensure the same resolution in both variables, (30) remains mandatory.
In terms of relative errors, σ κ � ε/κ 2 * , where κ 2 * � min(κ 2 , κ 2 ). e 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.

Eigenvalues Error Estimates via Rayleigh Quotient.
In the Neumann or Dirichlet condition, given a couple (κ 2 D , ψ) 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.
If ψ is an analytical solution, we get κ 2 I � κ 2 D , of course, but for the numerical solutions the difference between the two values provides a way to get a standard deviation by e quadratures were obtained by padding zeros to matrix M to ensure the analysis of the function at many different points and accurate integration. us, we use at least n � 36 for the quadrature process n � 36.

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. is empirical process, through inspection in sequence, almost matches the Rayleigh quotient's criteria.

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 κ 2 1 of the domain Ω with area A obeys the inequality where J 0,1 is the first root of the Bessel function of order 0 and J 2 0,1 � 5.783185962946773 corresponds to the first eigenvalue of the unit circle, therefore (ii) e eigenfunction ψ 1 associated with first eigenvalue κ 1 does not change its signal in Ω and obeys C Payne � 4π Ω ψ 2 1 dA/κ 2 1 ( Ω ψ 1 dA) 2 ≤ 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 κ 2 n+1 /κ 2 n ≤ 3, where it was conjectured that κ 2 2 /κ 2 1 attains its maximum for the circular regions i.e., κ 2 2 /κ 2 1 ≤ (J 1,1 /J 0,1 ) 2 � 2.5387 . . .. To follow the same form as the previous two constants, we write C Polya � κ 2 1 /κ 2 2 (J 1,1 /J 0,1 ) 2 ≤ 1 and the equality is valid only for the unit circle. e 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.

Asymptotic Behavior.
e asymptotic behavior of the eigenvalues in the Dirichlet condition is conditioned by Weyl Law and enhanced by Baltes [44]as follows: where n is the order of the eigenvalue counting the duplicates, L is the perimeter, and A is the flat domain area delimited by Γ(θ). e Weyl law can be used as a guide to inspect the limit when the numerical process is no longer defined to resolve high harmonics.

Examples for Helmholtz Problems
is 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).
is error analysis shows the spectral accuracy for hundreds of eigenvalues and eigenfunctions obtained in one shot.

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. e 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. e ellipse that delimits the domain is described by the equation 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. Figure 2 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.

Elliptical Domain-Neumann.
For the Dirichlet condition, let us now compare the results considering the Neumann boundary condition: ∇ → F · n � 0onzΩ for the same elliptic curve. Figure 5 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.
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.

Trefoil-Dirichlet.
is domain is delimited by the Jordan curve as follows: 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 200 th harmonic. Figure 8 illustrates the vibration modes from 76 to 100. Figure 9 shows the vibration mode 93 and Figure 10 shows the errors for eigenvalues and eigenfunctions.   e limacon Γ(θ) � 2.5 + cos θ+ 2 sin θ, with θ ∈ [0, 2π) ⊂ R, was studied by Atkinson [46] using "ridge polynomials" and other bases in a Galerkin approach.
is 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. Figure 11 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.

Pentafoil-Neumann Condition.
is 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 Figure 14 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.

Epitrochoid-Dirichlet Condition.
is curve was studied in a boundary value problem by Liu [48] under the Dirichlet  Mode 93 Figure 9: 93 rd eigenfunction in Dirichlet condition, whose eigenvalue is κ 93 � 19.94995891589. All digits match with [39].
condition. Here, we search for eigenvalues for the domain enlaced by this curve, whose form in polar coordinates is To the best of our knowledge, there have been no published results on these eigenvalues. erefore, we used the FE method to validate and compare the results. For the first 50 eigenvalues, the difference is less than 0.5%. Figure 17 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.
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).
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 than eps � 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.

Asymptotic Behavior-Dirichlet Condition
e inequalities referred to in Section 5.3, with the three constants C FK , C Payne , and C Polya , as well as the asymptotic behavior referred to in (24), are necessary conditions to be satisfied by the solutions. ese three constants must be less than unity. e 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 kappa 2 n,A be κ 2 n,A the asymptotic order N eigenvalue, given by (38), and κ 2 n be the numerical result for this eigenvalue, and we define ξ: 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. e behavior of ξ with respect to the harmonic order can reveal the reliability limit. Graphs of the behavior and constants

Conclusion
In this study, a diffeomorphic RRM is easily invertible, with all transformations resulting in PDE equations in two-dimensional planar domains. e 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.
ere is no need for a 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. e 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.

Data Availability
All the data used in this paper are available on the Internet by following the links given in the references.