Efficient Evaluation of Slowly Converging Integrals Arising from MAP Application to a Spectral-Domain Integral Equation

: In this paper, we devised an analytical technique to efficiently evaluate the improper integrals of oscillating and slowly decaying functions arising from the application of the method of analytical preconditioning (MAP) to a spectral-domain integral equation. The reasoning behind the method’s application may consistently remain the same, but such a procedure can significantly differ from problem to problem. An exhaustive and understandable description of such a technique is provided in this paper, where we applied MAP for the first time to analysis of electromagnetic scattering from a zero-thickness perfectly electrically conducting (PEC) disk in a planarly layered medium. Our problem was formulated in the vector Hankel transform domain and discretized via the Galerkin method, with expansion functions reconstructing the physical behavior of the surface current density. This ensured fast convergence in terms of the truncation order, but involved numerical evaluation of slowly converging integrals to fill in the coefficient matrix. To overcome this problem, appropriate contributions were pulled out of the kernels of the integrals, which led to integrands transforming into exponentially decaying functions. Subsequently, integrals of the extracted contributions were expressed as linear combinations of fast-converging integrals via the Cauchy integral theorem. As shown in the numerical results section, the proposed technique drastically outperformed the classical analytical asymptotic-acceleration technique.


Introduction
The classical statement of a general electromagnetic propagation, radiation, and scattering problem requires that fields, being solutions of Maxwell equations, satisfy boundary, edge (i.e., local power boundedness conditions), and radiation conditions [1]. Among the techniques devised to search for solutions to these kinds of problems, a special place is occupied by integral-equation approaches, i.e., integral-equation formulations associated with discretization techniques. First, this is because the radiation condition can be automatically satisfied through a suitable choice of the kernel of the integral equation; second, because integral equations and unknowns to be discretized can usually be defined on finite supports [2].
It is well known that the existence of a solution of an arbitrary integral equation cannot be mathematically stated and, if such a solution exists, that there are no theorems proving the convergence of an arbitrary discretization scheme [3,4]. This is what happens when dealing, for example, with open surfaces or nonsmooth boundaries for which problem formulation can lead to weakly singular integral equations of the first kind or integral equations of the second kind with strongly singular kernels.
Conversely, Fredholm theory [5,6] can be applied to an integral equation if the operator is the superposition of a continuously invertible operator and a completely continuous operator [7], meaning that, for these kinds of integral equations, the existence of an exact solution and the convergence of a general discretization scheme can be proven. For this reason, the scientific community has been seriously engaged in formulating integral equations of the mentioned type and discretization schemes generating well-conditioned matrix equations [8][9][10].
Methods of analytical regularization are aimed at converting first-kind integral equations and strongly singular second-kind integral equations to integral or matrix equations, for which the Fredholm theory is valid [11]. We achieved this by analytically inverting the most singular part of the integral operator and obtaining a Fredholm second-kind integral equation that could be solved via any direct discretization that keeps Fredholm's nature [12,13]. On the other hand, an analytically regularized matrix equation can be obtained in a single step through the suitable choice of a discretization scheme [14][15][16][17]. This is what happened when we used a complete set of orthogonal eigenfunctions of a suitable operator containing the most singular part of the original integral operator as the expansion basis of a Galerkin scheme. Such a method is called analytical preconditioning, as the Galerkin projection technique acts as a perfect analytical preconditioner for the considered integral equation.
We observed that the operator containing the most singular part of the integral operator to be diagonalized could be selected in different ways; hence, the rate of convergence could be different from scheme to scheme. In any case, the achievable numerical convergence of the method was limited at best by machine precision, or at least by the accuracy of the approximation in the numerical evaluation of matrix coefficients. Thus, if performed poorly, the latter point compromised the results and was generally an important factor to consider when building an efficient algorithm.
In the literature devoted to applications of the method of analytical preconditioning (MAP) to integral equation formulations, it has been widely shown that the choice of expansion functions reconstructing the physical behavior of unknowns guarantees fast convergence, i.e., fewer expansion functions are needed to achieve highly accurate results provided that a suitable accurate evaluation of matrix coefficients has been made [14][15][16][17][18][19][20][21][22][23][24]. Moreover, spectral-domain formulations were preferred in the quoted papers, as convolution integrals resulting from the Galerkin projection technique could be reduced to algebraic products by selecting expansion functions with closed-form spectral-domain counterparts.
The bottleneck of such a technique is that the obtained matrix coefficients are improper integrals of oscillating and, in the worst cases, slowly decaying functions to be numerically evaluated. A classical way to accelerate the convergence of these kinds of integrals consists of extracting asymptotic behavior from the kernels while expressing slowly converging integrals of the extracted parts in closed form [25][26][27]. Unfortunately, the integrands of the accelerated integrals are asymptotically oscillating functions with an algebraic decay, meaning that the choice of integration limits strongly affects the accuracy of that integration. Paradoxically, despite guaranteed convergence with respect to the matrix-truncation number, the algorithm becomes increasingly less efficient in terms of computation time as the required accuracy for a solution increases.
In a series of papers, the problem of the accurate and efficient numerical evaluation of these kinds of integrals was addressed under a different perspective [28][29][30][31][32][33][34][35][36]: (1) If needed, suitable contributions were pulled out of the kernels, making the integrands exponentially decaying functions; (2) an analytical technique was devised on the basis of the Cauchy integral theorem to express the integrals of extracted contributions as a linear combination of fast-converging integrals, or series. In this way, this method drastically outperformed the classical analytical asymptotic acceleration technique (CAAAT). Interestingly enough, despite using the same line of reasoning, the procedures devised in all the relevant papers were, in general, different from problem to problem.
In this paper, the technique detailed above is successfully applied for the first time to the electromagnetic scattering from a zero-thickness perfectly electrically conducting (PEC) disk in a planarly layered medium. The revolution symmetry of the problem allowed us to expand all the involved functions in terms of their series of orthogonal cylindrical harmonics. Hence, the problem was conveniently cast to a set of integral equations in the vector Hankel transform domain. In order to achieve preconditioning and, hence, fast convergence, each of the obtained integral equations was discretized via the Galerkin method, with the expansion functions reconstructing the physical behavior of the corresponding cylindrical harmonic of surface current density, and forming a set of orthogonal eigenfunctions on the basis of static parts of the associated operators. The obtained matrix coefficients, which are improper integrals of oscillating and slowly decaying functions, were rewritten as linear combinations of fast-converging integrals via the technique detailed above.
This paper is organized as follows. In Section 2, for the sake of brevity, we only provide a brief description of the problem formulation and the discretization of the general integral equation, referring the reader to the quoted papers for more details. An exhaustive and understandable description of the devised analytical technique for the accurate and efficient evaluation of the scattering matrix coefficients is presented in Section 3. Section 4 illustrates how the proposed technique drastically outperformed the CAAAT. Our conclusions are summarized in Section 5.

Spectral-Domain Integral Equation
In Figure 1, a zero-thickness PEC disk was located at the q-th interface of a planarly layered medium of L + 1 lossless, homogeneous, and isotropic layers of dielectric permittivity  The revolution symmetry of the problem allowed to expand all involved functions in series of orthogonal cylindrical harmonics Therefore, the problem could be equivalently reduced to an infinite set of independent onedimensional equations, obtained by imposing the n-th harmonic of the total electric field to be vanishing on the disk surface. In the vector Hankel transform domain, the integral equation for arbitrary index n was written as [37][38][39] for a ρ ≤ , where w is the Hankel mate of ρ, with respect to the argument [40]. The spectral domain Green function was given by [41,42]

Discretization Procedure
In order to discretize the obtained integral equations, the Galerkin method was used. Unknown were expanded in a complete series of Bessel functions [43]: where ( ) Following the line of reasoning in [39], it was possible to demonstrate that: (1) the obtained matrix equation was a Fredholm second-kind equation for which convergence of the approximate solution of the truncated matrix equation to the exact solution of the problem as the truncation order tends to infinity can be stated; (2) the behavior of the n-th harmonic of surface current density at the edge and around the center of the disk was correctly reconstructed, leading to fast convergence; (3) convolution integrals that resulted from the Galerkin projection technique being applied were automatically reduced to algebraic products, i.e., the matrix coefficients were expressed as linear combinations of one-dimensional integrals of the kind to be numerically evaluated. Unfortunately, the integrands of the integrals in Formula (7) were asymptotically oscillating and slowly decaying functions. Hence, the numerical evaluation of these kinds of integrals becomes increasingly less efficient in terms of computation time as the required accuracy for the solution increases. In order to accelerate the asymptotic decay of the integrands in Formula (7), suitable asymptotic contributions could be pulled out of the kernels so that the integrals of the extracted contributions could be expressed as linear combinations of Weber-Schafheitlin discontinuous integrals [25][26][27]44]: where 1 2 T P n k h p ≤ < + + + + and , , q T p g are asymptotic expansion coefficients.
However, the convergence rate was still strongly related to the accuracy required for the solution due to the asymptotic oscillating nature and algebraic decay of the integrands of the accelerated integrals. This problem was completely overcome through the application of the analytical technique presented in the next section.  (4), the Green function of two half-spaces was simply obtained:

Analytical Technique for Accurate and Efficient Evaluation of Scattering Matrix Coefficients
Starting from Equations (4j) and (4k), it was simple to verify that where  TM TE TM TE  TM TE TM TE  q  q  q  q  TM TE  q  TM TE TM TE  q  q  q   Z  R  R  R  S R R e were exponentially decaying functions. Hence, the integrals in Formulas (7) could be rewritten as the summation of two contributions where ( ) Formula (10) allows us to conclude that the integrals in Formula (13b) were fast convergent. On the other hand, the integrands of integrals in Formula (13a) were still asymptotically oscillating functions with an algebraic asymptotic decay. An analytical technique to fast evaluate these kinds of integrals, taking advantage of the nonoscillating nature of the kernels, is shown in the following.
Utilizing algebraic manipulations and the recurrence formula for Bessel functions [40] ( ) ( ) ( ) it was simple to obtain and the integration path in complex plane z w jw = + for integrals in Formula (16a) was the one in     Utilizing the Cauchy integral theorem [45], it was possible to write Starting from the behavior for the small and large argument of the Bessel functions of the first kind and Hankel functions [40] ( ) ( ) ( ) it was simple to conclude that functions in Formula (17) were bounded around 0 z = and decayed asymptotically as Therefore, Jordan lemma [45] allowed us to rewrite Formula (18) for 1 l = and 2 l = , respectively, as follows: where, due to relations [40] ( ) ( ) ( ) ( ) ( ) , the following expression could be established: Hence, by taking the difference between Formulas (20a) and (20b), we obtained where {} ℑ ⋅ denotes the imaginary part of a complex number.
By summing and subtracting integrals I ′′ , i.e., by considering the case in which the disk was located at the interface between two half-spaces.
For the sake of completeness, Table 1 shows the fast convergence of the presented technique. The value of the radial component of the surface current density on a disk of normalized radius with the x axis in the xy plane, impinged onto the scatterer surface. As clearly shown, convergence was exponential. We observed that: (1) values in Table 1 were independent of the required accuracy in the numerical evaluation of the integrals due to the smoothness of the integrands and (2) only 8 seconds were needed to reconstruct the solution with machine precision.   For the same example, values of ( ) 3, 0 J a ρ , calculated by using ten expansion functions via the CAAAT, can be seen in Table 2 for different values of the relative accuracy in the numerical evaluation of the integrals (RA). The extraction of only the first-order asymptotic behavior of the kernels of the integrals to be numerically evaluated was unthinkable in terms of computation time. For this reason, the first-and, whenever possible, second-order asymptotic behaviors were extracted from the kernels. It was clear that: (1) different results were obtained by changing the RA and (2) these values tended to be the convergent value obtained by means of the proposed technique as RA decreased. This last behavior could be better appreciated as shown in Figure 4: assuming the convergence value of ( ) 3, 0 J a ρ , obtained by applying the proposed method as reference value, the relative error in reconstructing the solution using the CAAAT can be seen plotted as a function of RA ( ( ) err RA ).
In summary, the accuracy of the solution obtained using the CAAAT increased as RA decreased. As outlined in the Introduction, this behavior was due to the algebraic decay and oscillating nature of the integrands of the accelerated integrals. Indeed, for these kinds of integrals, the choice of integration limits strongly affects integration accuracy.

RA
Jρ (a/3,0  A comparison between the computation times needed to reconstruct the solution, i.e., to fill in the coefficient matrix, via the CAAAT by extracting the first-and, whenever possible, second-order asymptotic behavior of the kernels of the integrals to be numerically evaluated, and by using the presented technique can be seen in Figure 5a for the same case as examined above. The CPU time ratio was plotted for different values of N and for three different values of RA. As expected, the CPU time ratio rapidly increased if higher accuracy was required for the solution, i.e., if a larger number of expansion functions were used, and if higher relative accuracy was required in the numerical evaluation of the integral coefficient matrix. In order to better emphasize the effectiveness of the proposed method, as can be seen in Figure 5b,c, the CPU time ratio obtained for two other cases, i.e., for 15deg.

Conclusions
In this paper, we presented an analytical technique for the efficient evaluation of improper integrals of oscillating and slowly decaying functions arising from MAP application to a spectraldomain integral equation. This is the first time this technique was applied to analysis of electromagnetic scattering from a zero-thickness PEC disk in a planarly layered medium. The proposed technique, which showed the numerical evaluation of improper integrals of exponentially decaying functions and proper integrals of bounded continuous functions, was very effective.
Indeed, it drastically outperformed the analytical asymptotic acceleration technique, which is widely used to compute improper integrals of oscillating and slowly decaying functions, independent of characteristics of the involved media and disk size.