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

Authors: Qinghua Wu; Mengjun Sun
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

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 ona,+∞, 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.


---

## Body

## 1. Introduction

In time-harmonic electromagnetic scattering, the following integral equations arise frequently [1, 2]:
(1)usρ=−ik4∫Cqρ′H01kρ−ρ′dl′,where H01x is the Hankel function of order 0 and qx is the unknown function.For scatterers with sharp edges or corners, the unknownqx should be sought in the following form [3]:
(2)qx=wxϕx,wx=1−xα1+xβ,where ϕx is smooth on −1,1, which leads to the following integral:
(3)∫−11wxH01kφxϕxdx.In addition, the following integral appears frequently in the fields of physics and engineering [4, 5]:
(4)which is equivalent to the following form:
(5)where wx=x−aαb−xβ and α>−1, β>−1 and fx 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 fm is Hölder continuous on a,b.In the casewx≡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–10]. A popular method is using the Taylor formula and eliminates the singularity at c by the following formula [7, 8, 11]:
(6)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 offx, Hasegawa and Torii presented an efficiently uniform approximation algorithms for following integrals [12, 13].
(7)However, it requiresON2 operations, where N−1 is the degree of Chebyshev interpolants. There is a traditional method of calculating these integrals.In the casewx≡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:
(8)where pNx=a0/2T0x+∑j=1N−1ajTjx+aN/2TNx and Tjx is the Chebyshev polynomial of the first kind.Moreover, by rewritingpNx−pNc/x−c in terms of Tjx as follows:
(9)pNx−pNcx−c=∑j=0N−1′djTjx.Equation (8) is transferred to
(10)where the prime denotes the summation whose first term is halved and Mjk=∫−11Tjxeikxdx can be calculated using a recursive formula.However, whenwx=x+1α1−xβ and α>−1, β>−1, the recursion formula for Mjk and dj is complicated. In [15], using the numerical steepest descent method, the following Cauchy principal value integral is calculated:
(11)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–22] and references therein.In this paper, we study the direct steepest descent method for a class of hypersingular integral.(12)where Gkx=eikx or Hv1kx,Hv1 denotes Hankel functions of the first kind of order v(v≤1).This paper is organized as follows: in Section2, 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.


## 2. 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:(13)In the caseGkx=eikx, the method is similar except for deleting an integral path around 0.SetD=z∈ℂ∣a≤Rz≤b,0≤Iz≤R, D1=z∈ℂ∣z−a≤r,0≤argz−a≤π/2, D2=z∈ℂ∣z−b≤r,π/2≤argz−a≤π, D3=z∈ℂ∣z−c≤r,0≤argz−c≤π, D4=z∈ℂ∣z≤r,0≤argz≤π, and D′=D1⋃D2⋃D3⋃D4, where R is a large number and r is small enough such that D′⊂D (see Figure 1).Theorem 1.
Suppose thatfz is an analytic function in the half-strip of the complex plane, a≤Rz≤b, and Iz≥0. If there are two constants M and k0 such that for 0≤k0<k(14)∫abfx+iRdx≤Mek0R,then the hypersingular integral with highly oscillatory kernels If,w,H;c,m,k can be calculated by the following formula:
(15)If,w,H;c,m,k=Saf,w;c,m,k+Sbf,w;c,m,k−1m!dmdcm−iπwcfcHv1kc,0≤v<1,Saf,w;c,m,k+Sbf,w;c,m,k−1m!dmdcm−iπwcfcHv1kc,v=1,+ikπw0f0limr⟶0Hv1krerθkrerθv,where
(16)Saf,w;c,m,k=iα+1∫0∞xαb−a−ixβa+ix−sm+1Hv1ka+ikxfa+ixdx,Sbf,w;c,m,k=−iβ+1∫0∞xβb−a+ixαb+ix−sm+1Hv1kb+ikxfb+ixdx.Proof.
SincewxHv1kx/x−cm+1fx is analytic in the region D−D′, using the Cauchy’s theorem, the following formula can be obtained:
(17)∫Γ1+Γ2+Γ3+Γ4+Γ5−Γ6−Γ7−Γ8−Γ9−Γ10wzHv1kxz−cm+1fzdz=0.
Setting(18)Ij=∫ΓjwzHv1kxz−cm+1fzdz,j=1,2,⋯,10,it derives that
(19)I1+I2+I3+I4+I5−I7−I9=I6+I8+I10,with all the contours shown in Figure 1.
Letz−a=reiθ,θ∈0,π/2, then
(20)I1=ir∫0π/2reiθαb−a−reiθβa+reiθ−cm+1Hv1ka+reiθfa+reiθeiθdθ≤r1+α∫0π/2b−a−reiθβa+reiθ−cm+1Hv1ka+reiθfa+reiθdθ⟶0,asr⟶0.
Similarly, we haveI5⟶0 as r⟶0 and
(21)I3=∫abiR+x−aαb−iR−xβiR+x−cm+1Hv1kiR+xfiR+xdx.
For large arguments, the Hankel functions behave like an oscillatory complex exponential with a decaying amplitude:(22)Hv1z~2πzeiz−1/2vπ−1/4π,−π<argz<π,z⟶∞.
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
(23)∣I3∣⟶0,asR⟶+∞.
In addition,I2 and I4 can be represented as
(24)I2=iα+1∫rRxαb−a−ixβa+ix−sm+1Hv1ka+ikxfa+ixdx,I4=−iβ+1∫rRxβb−a+ixαb+ix−sm+1Hv1kb+ikxfb+ixdx,I9=1m!dmdcm∫Γ9z−aαb−zβz−cfzHv1kzdz=1m!dmdcm−i∫0πc+1+reiθαb−c−reiθβfc+reiθHv1kc+reikθdθ⟶1m!dmdcm−iπwcfcHv1kc,asr⟶0.
Due to the Hankel function of order zeroH01z having a logarithmic singularity at z=0 and Hankel functions of higher order having algebraic singularities of the form 1/zv, z⟶0 [23], we achieve
(25)I7=∫Γ7z−aαb−zβkzvz−cm+1fzHv1kzkzvdz=−ir∫0πreiθ−aαb−reiθβkreiθvreiθ−sm+1freiθHv1kreikθkreiθveiθdθ⟶0,0≤v<1,−ikπw0f0limr⟶0Hv1krerθkrerθv,v=1,asr⟶0.
Thus, we obtain(26)If,w;c,m,k=limr⟶0,R⟶∞I1+I2+I3+I4+I5−I7−I9=Saf,w;c,m,k+Sbf,w;c,m,k−1m!dmdcm−iπwcfcHv1kc,0≤v<1,Saf,w;c,m,k+Sbf,w;c,m,k−1m!dmdcm−iπwcfcHv1kc,v=1,+ikπw0f0limr⟶0Hv1krerθkrerθv.
Formula (16) shows the integrals defined on 0,∞, and the integrand function decays exponentially. It can be calculated by the generalized Gauss-Laguerre quadrature rule. Suppose that xlα,wlαl=1N and xlβ,wlβl=1N are the nodes and corresponding weights of generalized Gauss-Laguerre quadrature rule, where the generalized weight functions are defined as xαe−x and xβe−x, respectively. Then, integrals in (16) can be approximated by
(27)Saf,w;c,m,k≈ikα+1∑l=1Nwlαb−a−ixlα/kβa+ixlα/k−cm+1fa+ixlαkHv1ka+ixlαexlα,Sbf,w;c,m,k≈−ikβ+1∑l=1Nwlβb−a+ixlβ/kαb+ixlβ/k−cm+1fb+ixlβkHv1kb+ixlαexlα☐Figure 1
The integration paths.


## 3. 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:(28)
where “MP” indicates that the error has reached the accuracy of the machine in Matlab. Numerical results are illustrated in Table1.
It can be seen from Table1 that the calculated result has an error of order O10−15, which indicates the accuracy of the proposed method is very high.Example 2.
Consider the computation of the hypersingular integral(29)
Its exact value can be expressed as follows [14]:
(30)I=i4dmdcmes1+ikE11−c1+ik−E1−1+s1+ik+ln−1−ik−ln1+ik+e−c−1+ikE11+c−1+ik−E11−c1−ik+ln1−ik−ln−1+ik.
It can be seen from Table2 that the calculation effect of the new method is very good, and as the frequency k increases, the calculation accuracy also increases.Example 3.
Consider the computation of the Cauchy principle integral(31)
Its exact value can be expressed by the following formula:(32)I=πk1F212;1,32;−k24,where F12a;b,c;z is the hypergeometric function.Example 4.
Consider the computation of the Cauchy principle integral(33)
For more general cases, Tables3–5 show that the accuracy of the algorithm is very high, and it is very effective.Table 1
Relative errors of the direct steepest descent method for.
Relative errorsk=10k=102k=103k=1057.1×10−15MPMPMPExact values-31.44003099-314.16506211-3141.5917622-314159.2653(Maple)497279211660231825375859263Table 2
Relative errors of the direct steepest descent method for.
Relative errorsk=10k=102k=103k=1052.8×10−14MPMP1.2×10−16Exact values-156.3459075115707.459596-1570795.8045-15707963267.4(Maple)00816749072314625343Table 3
Relative errors of the direct steepest descent method for.
Relative errorsk=10k=50k=100k=150MP4.4×10−168.9×10−164.4×10−16Exact values3.3521148732.83186970212.89862991062.9369583887(Maple)80775598748843122875Table 4
Relative errors of the direct steepest descent method for.
Relative errorsk=10k=20k=30k=401.02×10−153.09×10−162.62×10−165.88×10−16Exact values-0.8208250958210744-1.134059952837425-0.21469090499367690.5217520845428135(Maple)1.062845269936712i0.0389775446769420i0.8319561170161881i0.1919435601126078iTable 5
Relative errors of the direct steepest descent method for.
Relative errorsk=10k=20k=30k=409.99×10−163.21×10−161.57×10−165.8×10−16Exact values-0.8021025230103994-1.1913562s12289248-0.23155251704244590.5759125632671623(Maple)-1.0945266295060129i0.02163694829164146i0.8911805394195018i0.1976701245353936i



---

*Source: 8021050-2021-07-22.xml*