A Generalized Technique in Numerical Integration

Integration by parts is one of the most popular techniques in the analysis of integrals and is one of the simplest methods to generate asymptotic expansions of integral representations. The product of the technique is usually a divergent series formed from evaluating boundary terms; however, sometimes the remaining integral is also evaluated. Due to the successive differentiation and anti-differentiation required to form the series or the remaining integral, the technique is difficult to apply to problems more complicated than the simplest. In this contribution, we explore a generalized and formalized integration by parts to create equivalent representations to some challenging integrals. As a demonstrative archetype, we examine Bessel integrals, Fresnel integrals and Airy functions.


Introduction
It is well known that in applied mathematics and in the numerical treatment of scientific problems, slowly convergent integrals occur abundantly. They are produced by approximation procedures depending on a parameter, iterative methods, quadrature schemes, perturbation techniques and reliable evaluation of functions defined by integrals. These slowly convergent integrals present severe numerical and computational difficulties. Traditional quadrature rules and summation techniques have failed to provide accurate approximations to such integrals. Numerous methods and techniques were developed for improving convergence of these challenging integrals [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] and extremely efficient methods were introduced such as numerical steepest descent, Filon-type and Levin-type methods. Unfortunately, their application to complicated integrals is extremely challenging. Examples of such integrals are the twisted tail, Bessel integrals, Sommerfeld integrals and the so-called molecular multicenter integrals involved in molecular structure calculations.
The technique of integration by parts has frequently been used to create divergent asymptotic series representations of integrals. As an example, we consider the Euler series arising from integrating the Euler integral by parts: e-mail: hsafouhi@ualberta.ca In this example, the remaining integral after n integrations by parts: is often discarded and the divergent series has either been summed straightforwardly or summed through the use of sequence transformations. Another example of the use of integration by parts in numerical integration arises in [19], applied to the oscillatory spherical Bessel integral function involved in molecular integrals. Through a reformalized integration by parts with respect to x dx, we transformed the initial spherical Bessel integral into an integral involving the simple sine function: From a numerical point of view, the transformed sine integral function is more favorable than the initial Bessel integral function. The above integral transformation, which we caled the S transformation, was successfully applied to all molecular integrals leading to an unprecedented accuracy and efficiency [19,20]. However, this transformation requires the boundary terms to vanish at both limits of integration and was only applied to spherical Bessel integrals [21]. In the present work, we generalize this approach to a larger class of integrals. This generalization does not require the boundary terms to vanish at both limits of integration. The contribution that will be introduced is a reformalized integration by parts with respect to w(x) dx for some weight function w(x) = x µ for µ ∈ C whose choice depends on the integrand. For the successive differentiation d x µ dx k for k = 0, 1, 2, . . . required by this generalization, we will apply the Slevinsky-Safouhi formulae (SSF) [22].
We applied successfully the general reformalization of integration by parts to three well known challenging integrals, namely Bessel integrals, Fresnel integrals and Airy functions. The integral representations obtained for these integrals have better convergence properties and more favourable asymptotic behaviours.
For the computation of the integral representations, we propose a robust algorithm which will be referred to as the staircase algorithm. The algorithm uses both the boundary terms and the transformed integrals and progressively descends the integrand to an asymptotically favourable situation by applying the reformalized integration by parts at each iteration.

Slevinsky-Safouhi formulae (SSF) for differentiation
For more details on SSF and their applications, we refer the reader to [22,23]. Applications of SSF in numerical integration of challenging integrals are presented in [24,25].

well-defined and easy
to compute for m, n ∈ C. For µ, ν ∈ C, the term d Moreover, for m −1, these coefficients have the explicit expression: where (x) n,k is the Pochhammer k-symbol [26] and can be computed as The SSF 2 corresponds to the case (µ, ν, m, n) = (0, 0, 1, 0) which is the most popular.

Reformalized integration by parts
Definition 3.1 Let the following operator be known as integration of a function f (x) by w(x)dx, for some weight function w(x): Let us begin with the general integral f (x)dx for which the general antiderivative is F(x). By performing the substitution (x) ↔ x, we arrive at f ( (x))d (x). The differential element d (x) can be written as (x)dx and by denoting w(x) = (x), we assert that this becomes integration by w(x)dx, with a weight function w(x).
The integral f (x)dx generally does not have a simple antiderivative; however, the integral f ( (x)) (x)dx most certainly has an antiderivative with a simple representation with respect to the first integral, notably F( (x)).
Integration by parts is one of the main methods used in the formation of asymptotic series expansions of a given integral. When integrating by parts, we generally begin with: When integrating by parts by xdx as in [19], we begin with: Our method, which formalizes this process is suitable for integrals of composite functions, like the aforementioned f ( (x)).
In order to expand integrals of composite functions in asymptotic series, we must incorporate this formalism with integration by parts by w(x) with w(x) = x µ for µ ∈ Z. Naturally, then, we begin with: By using the weight function w(x) = x µ , we must continue by integrating by x µ dx. To clarify the following process, we make use of these functionals of g(x) and h(x): and The computation of the functionals G l (x) is obtained using SSF.
With these functionals, we can rewrite (10) as: and therefore with integration by x µ dx, the first term on the right hand side becomes a boundary term and the second term on the right hand side becomes the transformed integral.
and be of the general form: where the functionals G l (x) and H l (x) are defined as in (11). Let us consider the integral: By integrating the right hand side by parts, we obtain: And, by integrating by parts another n − 1 times, we obtain the equivalent representation: Naturally, then, the question arises in how to use the equivalent representations. Indeed, one could simply neglect the remaining integral as n → ∞ and form a series from boundary terms; or, one could evaluate the boundary terms up to m and continue to numerically evaluate the new integral. To these two proposed methods, there exist advantages and disadvantages. The boundary terms, for example, usually form divergent series, which are generally challenging to sum; however, that the integral becomes a simply series can lead to reduced calculation times. The new integral, on the other hand, may be asymptotically favourable; however, the computation of the new integrand may become challenging after a certain order.
The next algorithm we propose is actually a blend of the two extremes outlined above. We numerically integrate the integral to some point x = x 0 , a < x 0 < b. At such a point, we now apply the first iteration of the reformalized integration by parts by x µ dx to the remaining integral. We continue by evaluating the boundary terms at x = b and x = x 0 and we numerically integrate the new integral to another point x = x 1 , x 0 < x 1 < b. We continue in this manner until the desired accuracy is attained. This blended algorithm will be aptly known as the staircase algorithm, as it progressively descends the integrand into an evermore asymptotically favourable situation.
This methodology is concisely summarized in: For l = 1, 2, . . . , n, and for the sequence {x l } n l=1 that satisfies a < x l−1 < x l < b define: and the approximations to the integral b a f (x)dx form the sequence {S l } n l=0 .

Bessel Integrals
The integral that follows appeared in Numerical Recipes as a computational challenge and an example of application of sequence transformations: We begin by noting that the integrand f (x) satisfies f (x) ∈ C ∞ (−∞, +∞). This integral is already in an appropriate form for the reformalized integration by parts, since we can subdivide the integrand as the product of: Starting with these root forms, we can explicitly develop the functionals G l (x) and H l (x) as: through Bessel function properties. The integral then has the equivalent representation: In (22), all the boundary terms vanish. In Numerical Recipes, the value of the integral is determined through Cauchy's residue theorem to be I 1 = K 0 (1). Although the equivalent representation does not give this formula, it follows from (22) that: That is, each and every one of these integrals for parameterized n is also given explicitly by the Bessel function K 0 (1).
As previously alluded to, a numerical evaluation of the integrals I 1 , I 2 (a, v),Ĩ 2 (a, v) and I 3 (a, z) can be achieved through a multiplicity of methods. Historically, the popular candidates were the summation of the generally divergent boundary terms, although the transformed integral has also been used in cases where the boundary terms vanished. In this work, we implemented the Staircase algorithm, as the newest method to evaluate these integrals and serves as an original contribution to the overall methodology we will be developing.
The selection of an appropriate sequence {x l } n l=0 for the integration subintervals in the Staircase algorithm is critical. This sequence should leave the integral subintervals to be numerically integrable -the points cannot be spaced so far apart that a quadrature routine is incapable of providing the solution. Conversely, the sequence should not be so squished that the sum of the boundary terms and the integral subintervals introduces a reduction in accuracy from the addition of a very large number and a very large negative number. Generally speaking, the dependence on l of the sequence {x l } n l=0 resides in the integrand. For the cases at hand, the integrands are oscillatory, and thus if the difference x l − x l−1 is on the same order of magnitude as the difference between successive zeros of the integrand, both of the aforementioned scenarios can be averted.
As can be seen from the numerical table, the algorithm is capable of reaching high accuracy. This is mainly due to the fact that reformalized integration by parts progressively descends the integrand to an asymptotically favourable situation. Table 1. Numerical Results of Staircase Algorithm applied to Bessel integrals I 1 with x l = 2π(l + 1), Fresnel integrals with I 2 (0, 1). x l = √ 2π(l + 1), and Airy functions I 4 (0, 1) with x l = 3 √ 6π(l + 1).