Numerical Steepest Descent Method for Hankel Type of Hypersingular Oscillatory Integrals in Electromagnetic Scattering Problems

We present a fast and accurate numerical scheme for approximating hypersingular integrals with highly oscillatory Hankel kernels. The main idea is to first change the integration path by Cauchy’s theorem, transform the original integral into an integral on 
 
 
 
 a
 ,
 +
 ∞
 
 
 
 , and then use the generalized Gauss Laguerre integral formula to calculate the corresponding integral. This method has the advantages of high-efficiency, fast convergence speed. Numerical examples show the effect of this method.


Introduction
In time-harmonic electromagnetic scattering, the following integral equations arise frequently [1,2]: where H ð1Þ 0 ðxÞ is the Hankel function of order 0 and qðxÞ is the unknown function.
For scatterers with sharp edges or corners, the unknown qðxÞ should be sought in the following form [3]: where ϕðxÞ is smooth on ð−1, 1Þ, which leads to the following integral: In addition, the following integral appears frequently in the fields of physics and engineering [4,5]: which is equivalent to the following form: where wðxÞ = ðx − aÞ α ðb − xÞ β and α > −1, β > −1 and f ðxÞ is a given function. In order to ensure the existence of integrals, we need to assume that the m − 1-th order derivative of f is continuous and f ðmÞ is Hölder continuous on ½a, b.
In the case wðxÞ ≡ 1 and k = 0, many calculation methods have been proposed for (4), such as the Gauss-type method, the (composite) Newton-Cotes method, and others [6][7][8][9][10]. A popular method is using the Taylor formula and eliminates the singularity at c by the following formula [7,8,11]: Although these methods are simple and generally feasible, they also have some disadvantages (e.g., numerical cancellations and computation of the higher derivatives). By using Chebyshev interpolants of f ðxÞ, Hasegawa and Torii presented an efficiently uniform approximation algorithms for following integrals [12,13].
However, it requires OðN 2 Þ operations, where N − 1 is the degree of Chebyshev interpolants. There is a traditional method of calculating these integrals.
In the case wðxÞ ≡ 1, k ≫ 1, the integrand is highly oscillatory. The traditional calculation methods of these integrals have the disadvantages of low efficiency and poor accuracy and will encounter difficulties (4). Xiang et al. presented an efficiently uniform approximation scheme for this case [14]. The principle is to establish the following new approximate formula: where p N ðxÞ = ða 0 /2ÞT 0 ðxÞ + ∑ N−1 j=1 a j T j ðxÞ + ða N /2ÞT N ðxÞ and T j ðxÞ is the Chebyshev polynomial of the first kind.
Moreover, by rewriting ðp N ðxÞ − p N ðcÞÞ/ðx − cÞ in terms of T j ðxÞ as follows: Equation (8) is transferred to where the prime denotes the summation whose first term is halved and M j ðkÞ = Ð 1 −1 T j ðxÞe ikx dx can be calculated using a recursive formula.
However, when wðxÞ = ðx + 1Þ α ð1 − xÞ β and α > −1, β > −1, the recursion formula for M j ðkÞ and d j is complicated. In [15], using the numerical steepest descent method, the following Cauchy principal value integral is calculated: Numerical results show that the calculation effect of the proposed method is better, but the situation of hypersingular is not considered in the article. In [16], using Hermite interpolation and a recurrence formula, Liu and Xiang present a method for calculating integrals (4) in combination with the numerical steepest descent method and gave the error analysis. For more details on singular integrals with oscillatory function, see Ref. [15,[17][18][19][20][21][22] and references therein.
In this paper, we study the direct steepest descent method for a class of hypersingular integral.
where GðkxÞ = e ikx or H ð1Þ v ðkxÞ, H ð1Þ v denotes Hankel functions of the first kind of order v(v ≤ 1).
This paper is organized as follows: in Section 2, we review the basic formula for the steepest descent methods. In Section 3, the performance of the method is demonstrated by numerical examples, which verify the efficiency and accuracy of the algorithm.

Direct Steepest Descent Methods for
Approximating the Integral (12) In this section, we focus on the steepest descent method on the evaluation of the following integral: In the case GðkxÞ = e ikx , the method is similar except for deleting an integral path around 0.  Advances in Mathematical Physics where R is a large number and r is small enough such that D′ ⊂ D (see Figure 1). Theorem 1. Suppose that f ðzÞ is an analytic function in the half-strip of the complex plane, a ≤ RðzÞ ≤ b, and IðzÞ ≥ 0. If there are two constants M and k 0 such that for 0 ≤ k 0 < k then the hypersingular integral with highly oscillatory kernels Ið f , w, H ; c, m, kÞ can be calculated by the following formula: where Proof. Since ðwðxÞH ð1Þ v ðkxÞÞ/ððx − cÞ m+1 Þf ðxÞ is analytic in the region D − D ′ , using the Cauchy's theorem, the following formula can be obtained: Setting it derives that I 1 + I 2 + I 3 + I 4 + I 5 − I 7 − I 9 = I 6 + I 8 + I 10 with all the contours shown in Figure 1.
Let z − a = re iθ , θ ∈ ½0, π/2, then Similarly, we have jI 5 j ⟶ 0 as r ⟶ 0 and For large arguments, the Hankel functions behave like an oscillatory complex exponential with a decaying amplitude: As can be seen from the asymptotic behavior (22), for complex arguments with a positive imaginary part, the Hankel function decays exponentially, which follows that |I 3 | ⟶ 0, as R ⟶ +∞: ð23Þ

Advances in Mathematical Physics
In addition, I 2 and I 4 can be represented as Due to the Hankel function of order zero H ð1Þ 0 ðzÞ having a logarithmic singularity at z = 0 and Hankel functions of higher order having algebraic singularities of the form 1/z v , z ⟶ 0 [23], we achieve

Numerical Examples
In this section, we illustrate the accuracy and efficiency of the method described in this paper with several numerical examples. All the calculations are done on Matlab R2016b, and all the exact values of the integrals are calculated by Maple with 32-digit arithmetic. The number of nodes used in generalized Gauss-Laguerre quadrature rule is set to 32. For the detail of the algorithm, see Appendix.
Example 1. Consider the calculation of the following Hadamard finite part integrals: where "MP" indicates that the error has reached the accuracy of the machine in Matlab. Numerical results are illustrated in Table 1.
It can be seen from Table 1 that the calculated result has an error of order Oð10 −15 Þ, which indicates the accuracy of the proposed method is very high.
Advances in Mathematical Physics

Example 2. Consider the computation of the hypersingular integral
Its exact value can be expressed as follows [14]: It can be seen from Table 2 that the calculation effect of the new method is very good, and as the frequency k increases, the calculation accuracy also increases.
Its exact value can be expressed by the following formula: where 1 F 2 ða ; b, c ; zÞ is the hypergeometric function.