FAST APPROXIMATION OF ALGEBRAIC AND LOGARITHMIC HYPERSINGULAR TYPE SINGULAR INTEGRALS WITH HIGHLY OSCILLATORY KERNEL

Herein, highly oscillatory integrals with hypersingular type singularities are studied. After transforming the original integral into a sum of line integrals over a positive semi-infinite interval, a Gaussrelated quadrature rule is constructed. The vehicle utilized is the moment’s information. The comparison of two algorithms (Chebyshev and its modified one) to produce the recursion coefficients that satisfy orthogonal polynomial with respect to Gautschi logarithmic weight function, is investigated. Lastly, numerical examples are given to substantiate the effectiveness of the proposed method.


Introduction
The main aspect of our interest in this paper is the computational efficiency of hypersingular integral of the form where α, β > −1, s ∈ {0, 1} , |ω| >> 1, a < µ < b, −∞ < a < b < ∞ and f is an analytic function of x in an appropriately large region containing the interval [a, b] . For simplicity let the above integral be written as where ω (α,β,s) (x) = (x − a) s is the weight function on [a, b] . While analyzing the above integral, we encounter the following cases: 1. When s ∈ {0, 1} and k = 0 the integral becomes algebraic (s = 0) or logarithmic (s = 1) singular highly oscillatory integral, 2. When k = 1 the integral is considered as the Cauchy Principal Value (CPV) integral with a highly oscillatory kernel, 3. When k = 2 the integral is understood in the Hadamard Finite Part (HFP) sense integral with a highly oscillatory kernel, 4. When k ≥ 2 the integral is considered as Hypersingular, highly oscillatory integral.
Typically, it is somewhat difficult to evaluate CPV or Hypersingular as ordinary improper integrals due to the existence of a pole at x = µ. If the order of singularity of the kernel function is higher than the dimension of integral, then the integral itself is called hypersingular. Having said that, hypersingular integrals play a vital role in physics and applied mechanics, where they can be utilized to formulate mixed boundary value problems. Besides these, we may find several applications of finite part integrals also in electrodynamics, aerodynamics, acoustics, and so on, for more about their applications, please see [1][2][3][4][5].
In an attempt to evaluate the Cauchy Hyperbolic problem for the differentiation of several variables, Hadamard [6,7] introduced the concept of the finite-part in hypersingular integrals. Consequently, numerous approximations have been achieved, since then, to obtain the fastest algorithm for hypersingular integrals, particularly in the case where ω (α,β,s) = 1 and the frequency (ω = 0), for example, Gauss related types quadrature [8][9][10][11][12][13]. One may simply rewrite (1.1) in terms of the kth order derivative of the Cauchy principal value integral [11] as In the integral part of the above equality when ω (α,β,s) (x) = 1 and the frequency ω = 0 the integration part becomes the Hilbert transform and the sufficient condition for the existence of Hilbert transform is that f (x) has to satisfy a Lipschitz and Hölder condition in a closed interval [a, b] . In recent decades, many scientists have studied numerous computation methods and applications for numerical evaluation of hypersingular integral equations (for more details see [17][18][19][20][21][22][23] and the references therein). Since such types of integrals contain a broad range of applications, it is of great interest to provide a fast and accurate algorithm for the numerical computation of these types of integrals.
It is essential to mention that our primary goal here is to provide a fast algorithm for the efficient computation of the algebraic and logarithmic hypersingular type integrals with highly oscillatory kernel. Based on the assumption that f is a holomorphic function in an appropriately closed large region consisting [a, b], the original integral is changed into several line integrals in a positive semi-infinite interval. Subsequently, these line integrals are computed using the construction of a Gauss related quadrature rule. Moreover, algorithms to produce the recursion coefficients of the three-term recurrence relation that satisfy an orthogonal polynomial with respect to Gautschi logarithmic weight function over positive half-infinite interval are compared by employing the moment's information approach.
This paper is structured as follows. In the next section, we evaluate the integral (1.1) using the proposed method. In section 3, we give the construction of the Gauss-related quadrature rule. Section 4 comprises numerical experiments to substantiate the efficiency of our proposed algorithms. Concluding remarks are provided in Section 5.

2.
A method for the evaluation of the integral (1.1) In this section, we demonstrate a method based on contour integration in the complex plane for the integral (1.1). For k, s = 0, the application of the steepest descent or Clenshaw-Curtis methods are straightforward [24,25]. However, for k ≥ 0 and s = 1, some modifications will be required inasmuch as some of the integrands will contain multi-valued functions that cannot be computed directly. Here, let z be defined as z = |z| e i(θ+2nπ) , for n = 0, ±1, ..., to get a single-valued function, we define our principal branch as ln |z| + iθ, 0 ≤ θ ≤ 2π. Recall that the principal value of ln z is the value obtained when n = 0. For instance, ln (±i) = ln 1 + 2n ± 1 2 πi, which is equal to ± π 2 i for n = 0. Γ j , such that (2.1) The above theorem (Cauchy-Goursat), which is the focal theorem in this section, is used to prove the following theorem: Theorem 2. Let f be analytic inside and on a simply-connected curve in a complex plane, briefly denoted by , in Fig. 1. For a sufficiently large L, suppose that there exist two nonnegative constants M and ω o such that Γ |f (x + iL)| dx ≤ M e ωoL , for ω o ∈ [0, ω). Then, the hypersingular oscillatory integral (1.1) can be transformed in the following form: Proof: Let's denote the integrand φ (z) as and let be the region in the domain D = {z ∈ C : a ≤ Re (z) ≤ b, 0 ≤ Im (z) ≤ L} , bounded by Γ as defined in (2.1). Since φ (z) is analytic, by employing the Cauchy-Goursat integral theorem, we have with the direction being taken in the positive sense (counterclockwise), as depicted in Fig. 1. We can show, with ease, that the integral over quarter circles HA and DE tend to zero as epsilon goes to zero. More precisely, let Γ 1 : z = a + εe ix , x ∈ 0, π 2 , this yields It can be easily shown that the function G (t, θ) is continuous on t ∈ [0, ε] and θ ∈ 0, π 2 . Thus, taking the limit on both sides as ε → 0, we get Γ1 φ (z) dz → 0. The same technique can be applied to obtain wherein the above expressions Note that for any nonnegative constants M 1 , ω 1 such that ω > ω 0 + ω 1 and for sufficiently large L, we can Also, on Γ 6 , letting Γ 6 : z = b + ix for x ∈ [ε, L] and applying the change of variable where ωx = y, we get similarly on Γ8, we have Integral over Γ3 : The above right-hand side integral part gives as ε → 0, the integral over Γ3 becomes Taking the limits as ε → 0, and L → ∞ of the above results and substituting all in (2.5) completes the proof.

Computation of (2.2) by Gaussian quadrature rule
In this section, we construct a Gaussian rule for efficiently computing the integral (2.2). The technique applied in this paper is based on the recent development of Walter Gautschi in his paper [28], followed the paper by Ball J.S. et al. [29]. For simplification let us rewrite (2.2) as ... Drawing from the Gautschi analysis, we can transform the complex logarithmic multi-valued function that appeared in (3.1) as Substituting the above equations (3.2) and (3.3) into (3.1) leads us to the new form of (3.1) which can be denoted as Moreover, analyzing (3.4) we came across the following related (nonnegative) weight functions: Here w (α,LG) and w (β,LG) should be recognized as logarithmic Gautschi-Laguerre weight functions over a positive half-infinite interval whereas the rest are known to be the conventional Generalized Gauss Laguerre weight functions y → w (y; m) = y m e −y for (m = α, β > −1) .
By combining logarithmic Gautschi-Laguerre with the traditional Generalized Laguerre weight functions, a Gaussrelated quadrature rule for the integral (3.4) can be formulated as follows In the above summations (3.6) y Provided that u is orthogonal to v if (u, v) = 0 and the norm of u is obtained if u = v which yields The monic orthogonal polynomials with respect to (3.7) are denoted by πi (.) = πi (., w (y)) , i = 0, 1, and they satisfy the three-term recurrence relation Furthermore, Golub  yp. Furthermore, in the determination of the eigenvalues of the matrix J, the QR method proposed by Francis [32] is deemed the best method available since it converges very rapidly for symmetric matrices.
There exist several methods for generating the n recursion coefficients. As stated previously, we employ the method based on the moment's information.
It is well known that the first n recursion coefficients αi, βi, i = 0, 1, .., n − 1 can be evaluated by the first 2n moments µr = b a y r w (y) dy, r = 0, 1, ...., 2n − 1. Consequently, the moments associated with the weight functions in (3.5) (taking α as an example) can be calculated exactly in an efficient way as where Γ (z) = ∞ 0 y z−1 e −y dy is the Gamma function and ψ (z) is the logarithmic derivative of the Gamma function ψ (z) = Γ (z) /Γ (z) known as Digamma function [33], and it can be implemented in Mathematica as PolyGamma[z].
Several formulas exist, for instance, by expressing the coefficients in terms of Hankel determinants in these moments, we obtain the n recursive coefficients. However, the sensitivity of the formula to small errors increases with n due to the condition number that grows with n if an extended precision is not employed [34]. For constructing recursive coefficients that satisfy three-term recurrence relation for orthogonal polynomials with respect to the given weight functions w (y; m) = y m e −y and w (y; m) = (y − 1 − ln y) y m e −y for m = α, β > −1, the moments obtained in (3.11) are sufficient.  Another alternative is to utilize the modified moments where instead of evaluating moments by the power y r in In this case, letL m i (y) be the monic generalized Laguerre polynomials [35] and suppose that L m i are monic polynomials satisfying a three-term recurrence relation similar to (3.9) In the above recurrence relation, the coefficients ai ∈ R, bi ≥ 0 can be determined in the following way: We already know that Equating (3.12) and (3.14) we observe that the coefficients ai ∈ R, bi ≥ 0 are (3.15) ai = 2i + m + 1 and bi = i (i + m) , b0 = µ0 = Γ (m + 1) .
Thus, the modified moments related to the weight function w (m,LG) (y) = y m e −y (y − 1 − ln y) , over [0, ∞) for (m = α, β > −1) is expressed in the Gamma and its logarithmic derivative function as Thanks to the modified Chebyshev algorithm [ [36], p. 76], which directly takes 2n modified moments and the in (3.12) to produce the coefficients αi (w) , βi (w) , i = 0, 1, ..., n − 1 that we strongly desired. In order to generate coefficients using the modified moments' approach, we provide an algorithm that can be implemented in any mathematical software. Our implementation was done using Mathematica version 9.0.
10: Evaluate βi = λi,i / λi−1,i−1 ; 11: Return αi; βi for i = 0, 1, ..., n − 1; Coefficients obtained from Algorithm 2 by using modified moments are in full agreement with the ones provided by Algorithm 1 using ordinary moments. Although both algorithms lead to the same results, the first algorithm requires the employment of extended precision compared to the second algorithm, resultantly, a slight computation time difference between the two algorithms may be observed. This problem is due to the ill-conditioning that grows exponentially with n.
for α = −0.3, β = 0.6, µ = 0.5. Table 3 shows the approximation values obtained using the proposed method with different values of n and ω fixed. Moreover, Table 3 exhibits that as n increases, the computation accuracy also increases.
for α = −0.31, β = 0.69, µ = 0 using the two algorithms. Table 5 shows the errors and CPU time in seconds using the proposed algorithms with ω = 350 (frequency) fixed. In Table 5, precision was fixed and equal to 40 for each algorithm. Table 6 exhibited that the accuracy of the two algorithms improved with n, although it can be seen that as n increases, the computation time in Algorithm 1 increased proportionally. Moreover, all the results were found to be entirely satisfactory.
Example 5. Consider the following algebraic singular Hadamard finite part: (4.5) Table 6 shows Errors and CPU time in seconds using the proposed method with ω = 100 fixed and µ = 0. Additionally, Table 7 shows the approximation values for µ = 0.5 using the method proposed.  Taking n = 20 as a reference, we can easily see that as n increases, the computational accuracy also increases.
Example 6. We consider the computation of  Table 8 exhibits the accuracy of the proposed method for n = 4 fixed. Furthermore, Table 8 shows the absolute errors and the executed time in seconds for both functions f1 and f2 with ω = 10 2 ,10 3 ,10 4 ,10 5 ,10 6 . We observe that the proposed method can guarantee at least |errors| ≤ 10 −15 for both functions f1 and f2. It can be also seen through Table 8 that the proposed method is more efficient for large values of frequency.

Conclusions
Herein, the computational efficiency of algebraic and logarithmic hypersingular type singular integrals with highly oscillatory kernel was discussed. The approach employed proved that it can easily yield better approximation accuracy as n increases. Furthermore, the used approach also exhibited that taking n as fixed, increasing the frequency resultantly increases approximation accuracy. Also, two algorithms to produce the recursion coefficients that satisfy orthogonal polynomial with respect to Gautschi logarithmic weight function were given and compared.
Moreover, both algorithms showed that they can give satisfactory results when one among them is used. The tables and convergence rate figures given hereinbefore provide evidence in support of our analysis.

Conflicts of Interest:
The author(s) declare that there are no conflicts of interest regarding the publication of this paper.