Two formulae with nodes related to zeros of Bessel functions for semi-infinite integrals: extending Gauss–Jacobi-type rules

For approximating integrals ∫0∞xαf(x)dx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\int _0^{\infty }\!\!} \varvec{x}^{\varvec{\alpha }}\varvec{f(x)dx}$$\end{document} (α>-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\alpha >-1}$$\end{document}) over a semi-infinite interval [0,∞)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{[0,\infty )}$$\end{document} with a given function f(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{f(x)}$$\end{document}, two formulae, one of them new and another associated with an existing formula, are presented. They are constructed in a limiting process to a semi-infinite interval [0,∞]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{[0,\infty ]}$$\end{document} with a linear transformation for a well-known approximation method, the Gauss–Jacobi (GJ) rule and its family rules: the Gauss–Jacobi–Radau (GR) and Gauss–Jacobi–Lobatto (GL) rules on a finite interval. This procedure was used in constructing our previous limit Clenshaw–Curtis-type formulae. The limit GJ formula (LGJ) constructed in this way uses as nodes zeros of the Bessel function Jα(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\varvec{J_{\alpha }(x)}}$$\end{document} squared after multiplied by a positive constant a and the limit GR (LGR) (and limit GL (LGL)) formula those with zeros of Jα+1(x)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{J_{\alpha +1}(x)}$$\end{document}. The LGJ formula is also shown to be derived from the formula developed by Frappier and Olivier (Math. Comp. 60:303–316, 1993) for an integral on [0,∞)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{[0,\infty )}$$\end{document}. The LGR and LGL formulae give the same and new formula. We show that for a function f(z) analytic on a domain in the complex plane z and satisfying some appropriate conditions, there exists a constant d>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${d>0}$$\end{document} such that the errors of both formulae are O(e-2d/a)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{O(e^{-2d/a})}$$\end{document} as a→+0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{a\rightarrow +0}$$\end{document}. The average of the LGJ and LGR formulae gives smaller quadrature errors than each formula. Numerical examples confirm these behaviors and show that the LGJ and LGR formulae give asymptotically the same quadrature errors of opposite sign. Consequently, the LGR formula behaves like an anti-LGJ formula in the same way as the Lobatto rule for integrals on [-1,1]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{[-1,1]}$$\end{document} behaving like the anti-Gauss rule.


Introduction
For a given function f (x) defined on a semi-infinite interval [0, ∞) and a constant α > −1, we want to approximate an integral I f over [0, ∞), with the assumption that x α f (x) decays sufficiently fast as x → ∞.Two formulae, one of them being new and another associated with an existing formula, are presented.They are constructed by extending interpolatory quadrature formulae (cf.[1, 2.5]) for integrals with the Jacobi weight function In our previous work [2], the C-C rule and its four family rules (cf.[3]) are extended to the semi-infinite integral I f in (1.1) with a case α = 0, to develop five new formulae.
In the present work, analogous extensions of the Gauss-Jacobi (GJ) (cf.[1, p. 95]) and its family rules, the GJ-Radau and GJ-Lobatto rules (cf.[1, pp.103-104]) with degree of exactness higher than the C-C and its family rules, are attempted in order to develop two approximation schemes to the integral I f in (1.1); the limit GJ (LGJ)(see (1.2)) and the limit GJ-Radau (LGR) (see (1.4)) formulae.It is also shown that the LGJ formula is derived from the formula involving zeros of the Bessel function J α (z) developed by Frappier and Olivier [4].The LGR formula is, on the other hand, a new one.Further, the limit GJ-Lobatto (LGL) formula coincides with the LGR formula.
Throughout the paper, denote by α, β and h real constants such that α > −1, β > −1 and h > 0, and n an integer ≥ 1.Let 0 < j α,1 < j α,2 < . . .and j α,−l = j α,l (l = 1, 2, . . . ) be zeros of the Bessel function J α (z) of the first kind of order α (cf.[5]).An approximation to I f in (1.1), the limit Gauss-Jacobi (LGJ) formula I LGJ h f , is given by where nodes x LGJ h,l and weights w LGJ h,l are given, respectively, by (1.3) Appendix A shows that this LGJ formula is derived from the formula (8) in Frappier and Olivier [4].The summation in the center of (1.2) can be used as a numerical integration formula.
The limit Gauss-Radau-Jacobi formula (LGR) I LGR h f that is a new approximation scheme with not only interior nodes, but also the boundary node at the end point 0 of [0, ∞), for I f in (1.1), is given by where nodes x LGR h,l and weights w LGR h,l are given, respectively, by ) (1.6) with j α+1,l (l ≥ 1) being zeros of the Bessel function J α+1 (z).
Remark 1 Comparison of (1.3), (1.5) and (1.6) reveals a remarkable similarity between the LGJ and LGR formulae.The nodes x LGR h,l (l ≥ 1), except for x LGR h,0 , are x LGJ h,l with α replaced by α + 1 and the weights w LGR h,l (l ≥ 1) are w LGJ h,l with α of J α and j α,l replaced by α + 1.
Lemma 1. 1 The LGL formula coincides with the LGR formula.
Frappier and Olivier [4] show that their formula coincides with the midpoint rule when α = −1/2 and that the trapezoidal rule is induced from their formulae with α = 1/2.For the case α = −1/2, our results are summarized below.Proposition 1.2 When α = −1/2, with the change of variables x = t 2 for I f in (1.1) and the definition g(t) = f (x) = f (t 2 ), the LGJ and LGR formulae coincide with the midpoint and trapezoidal rules, respectively, g(hl), (1.8) where the prime denotes the summation whose first term is halved.
Asymptotic convergence behaviors of the LGJ and LGR formulae are briefly discussed here.To this end, we introduce a function family.Throughout the paper, for a complex number z, we write z = x + iy (x, y ∈ R).Denote by P d the right region of a parabola in the complex plane z, with a number d > 0 (see Fig. 1) and by ∂ P d the contour traversing the boundary of P d in the counterclockwise direction, ) We start by defining a complex function z ν with a real number ν, by According to Abramowitz and Stegun [6, Section 9] and the DLMF [7], for a complex number τ ∈ C, we define the Bessel functions α (τ ) and H (2) α (τ ) (−π < arg τ < π) that are analytic on the τ -plane cut along the negative real axis.With these Bessel functions, we define α (τ ) by H (1)  α (τ ) J α (τ ) (0 < arg τ < π).
(1.14) Next, an inequality is introduced to discuss the errors of LGJ and LGR formulae.For the notation of order O(•) and a function F(•), we define the inequality This means that there exists a function g The error of the LGJ formula is given below.
where 0 < arg z < 2π for z ∈ ∂ P d in (1.10).The error is bounded by Note that the term z α α (π z 1/2 /h) in the integral in (1.15) is analytic for 0 < arg z < 2π except for the positive real axis in the complex plane z.
The error of the LGR formula is given below.

Theorem 1.5 If f ∈ B(α, d), then the error E I
LGR h f of the LGR formula in (1.4) is expressed as a contour integral: (1.17 where 0 < arg z < 2π for z ∈ ∂ P d in (1.10).The error is bounded by Theorems 1.4 and 1.5 show that for a function f ∈ B(α, d), the errors of the LGJ and LGR formulae converge exponentially in O(e −2π d/h ), twice as fast as O(h j e −π d/h ) ( j = 1, . . ., 5) for the limit C-C-type formulae, as step size h → +0 (cf.[2]).This ratio of the convergence rates for two groups of formulae implies that the LGJ and LGR formulae inherit the nice property of the Gauss-Jacobi-type rules that converge asymptotically twice as fast as the C-C-type rules (cf.[8]).
Define the average I LAV h f of the LGJ and LGR formulae by Then, I LAV h f converges faster than the LGJ and LGR formulae as h → 0, as shown in Theorem 1.6.Remark 2 Laurie [9] proposes an anti-Gauss rule to estimate the error of the Gauss rule to the integral 1 −1 f (x)dω (cf.[10,11]).The average of the Gauss and anti-Gauss rules is also shown to give a better approximation.For the n-point Gauss and (n + 1)-point anti-Gauss rules that have the same degree of exactness 2n − 1, the definition of the anti-Gauss rule means that their average has degree of exactness higher by at least two, 2n + 1.Further, a generalized averaged Gaussian formula is proposed and shown to have the maximum degree of precision 2n + 2 at least, and 2n + 3 for the special case where the measure dω is symmetric with respect to the origin; see [11] and references therein.
By adopting these properties of the average rules, namely, degree of exactness higher than the Gauss rule, we define an asymptotic anti-formula for a semi-infinite integration formula.Assume that A and B formulae have the same order of convergence and their average has higher order of convergence.Then, we call A (or B) formula asymptotic anti-B (or anti-A) formula.On the other hand, the LAV formula converges faster by the factor h than the LGJ and LGR formulae; see (1.21) below.
In this sense, the LGR formula is the asymptotic anti-LGJ formula.
Theorem 1.6 Define α (τ ) by ) is expressed as a contour integral: where 0 < arg z < 2π for z ∈ ∂ P d in (1.10).The error is bounded by

.21)
When α = −1/2, the leading term vanishes, and we have These properties are experimentally confirmed for the integral I f in (1.1), where f (x) = e −x /((x + 1) 2 + 1/4) and α = 0 and 5.The left image of Fig. 2 depicts the error behaviors of the LGJ, LGR, their average LAV and the limit C-C (LC-C) formulae for α = 0.This image shows that the LGJ and LGR formulae converge twice as fast as the rate O(h 3 e −π d/h ) of the LC-C formula.The LAV formula gives approximations of higher accuracy by about two significant digits than the LGJ and Fig. 2 (left) Errors of the approximations to the integral I f in (1.1) by the LGJ and LGR formulae, their average (LAV) and the limit C-C (LC-C) formula, where f (x) = e −x /((x + 1) 2 + 1/4) and α = 0. (right) Three error orders are shown for the LGJ, LGR and LAV formulae in black and red curves.The error ratios LGR formulae.The right image of Fig. 2 shows that the LGJ and LGR formulae give the same quadrature errors of opposite sign except for small 1/h.This supports that the LGR (equivalent to the LGL) formula is asymptotically an anti-LGJ formula and inherits the property of the Gauss-Lobatto rule that is the modified anti-Gauss rule (cf.[10]); see also [1, p. 106].Figure 3 depicts the error behaviors of the LGJ, LGR and LAV formulae for the same function f as in Fig. 2 with α = 5.
The computation in this experiment is performed in multiple precision arithmetic of seventy-five significant digits with Mathematica.
The paper is organized as follows.The LGJ formula (1.2) is derived in Section 2 and the LGR formula (1.4) in Section 3. Section 4 shows that the LGL formula coincides with the LGR formula.In Sections 3 and 4, we exploit expressions for the nodes and the weights of the GJ-Radau and GJ-Lobatto rules constructed by Gautschi [12,13], although for the weights of both rules at interior nodes, we develop and use alternative 5 Fig. 3 (left) Errors of the approximations to the integral I f in (1.1) by the LGJ and LGR formulae and their average (LAV), where f (x) = e −x /((x + 1) 2 + 1/4) and α = 5. (right).The error ratios r (h 123 expressions suited to our purpose; see Lemmas 3.1.and 4.1.Theorems 1.4 and 1.5 are proved in Section 5. Theorem 1.6 is proved in Section 6. Section 7 concludes the paper.Appendix A shows that the LGJ formula is also derived from the Frappier-Olivier formula.Proposition 1.2 is proved in Appendix B.

Deriving the limit Gauss-Jacobi (LGJ) formula
In this section, the LGJ formula (1.2) is derived from the Gauss-Jacobi rule.For real numbers a and b and an integer n ≥ 1, we begin by defining γ n (a, b) as follows

The Gauss-Jacobi (GJ) rule and its limit formula
We consider the integral with the Jacobi weight The n-point GJ rule 2) is given by where the nodes 1 > ξ > −1, in decreasing order, are zeros of the Jacobi polynomial are given by ω The LGJ formula is derived from the GJ rule in (2.3) as following.Define a n by and the integral I a n f n by Then, we see that lim Applying the change of variables x = a n (1 − t)/2 in (2.6) gives Using the GJ rule (2.3), where we set β = 0, to the integral in (2.7), we have the approximation where nodes x and weights w are given, respectively, by ), (2.9) (2.10) The formula 2) is given by the limit of the formula The proof is given below.

Proof of Theorem 2.1
We begin with some preliminary results.The key one is Theorem 2.2 below that gives a relation between the Jacobi polynomial P (α,β) n (z) of degree n and the Bessel function J α (z); see Szegö [14, p. 192] and (7) in [4].

Theorem 2.2 (Formula of Mehler-Heine type
(2.12) Let l be an integer l ≥ 1 and ξ be a zero of P (α,β) n (z) and j α,l a zero of ) is a zero of In view of Theorem 2.2 and Hurwitz's Theorem (cf.[15, p. 231]), it follows that lim Proof Differentiating the both terms in (2.12) with respect to z before multiplying 1/z on the both sides gives where the convergence is uniform in a bounded domain.Insert n 2(1 − ξ ) in place of z in the left-hand side of (2.15).Then, recalling (2.13) and inserting a zero j α,l of J α (z) in place of z in the right-hand side of (2.15), we have (2.14).Now, we prove Theorem 2.1.We start by verifying the first relation in (2.11).Since, in view of (2.13), we see that In view of (1.3), (2.5) and (2.16), we have We verify the second relation in (2.11).In view of (2.16), it follows that (2.18) In view of (2.14), we see that (2.20)

Deriving the limit Gauss-Radau-Jacobi (LGR) formula
In this section, the LGR formula (1.4) is derived from the Gauss-Radau-Jacobi rule in the same limiting process as in Section 2. Gautschi [12] gives the (n + 1)-point Gauss-Radau-Jacobi rule by where the nodes, 1 = ξ with ξ being the zeros of the Jacobi polynomial P (α+1,β) n (x) of degree n, and the weights ω ) Note that our weights in (3.3) and (3.4) are of slightly different forms from those in Gautschi [12], where the preassigned boundary node is the lower end of the integration interval, while it is the upper end in our case, as shown in (3.2).Specifically, by noting that the right-hand side of (3.1) is written by < 1, exchanging the roles of α and β in [12] gives our results above.In Section 4, note this type of differences as well.
For the weights ω ), however, we use an alternative and simpler expression given below.

Lemma 3.1 For the weights ω
where ω are the weights for the GJ rule for the integral Then, g is a polynomial of degree 2n.Since the (n + 1)-point GJ-Radau rule is exact for the polynomial g, in view of (3.2) and the fact that g(1) = 0, we have ).
(3.6)Then, the rightmost-hand side of (3.6) is the n-point GJ rule with the weights ω . Hence, we have (3.5).
Using the RJ rule (3.1) to the integral I a n f in (2.7), we obtain the approximation where the nodes x are given by and the weights w Here, we used (3.5) for ω and (2.10) for w

The formula I
LGR h f in (1.4) is given by the limit of the formula Next, we verify the second relation in (3.11) with the two relations in (1.6).It is easy to verify the first relation in (1.6) for l = 0. Indeed, in view of (2.1) and (3.9), and recalling that a n = (2hn/π ) 2 , we see that (n → ∞).
For l ≥ 1, using (2.16) in (3.10) and in view of (2.20), we see that

Limit Gauss-Lobatto-Jacobi (LGL) formula coincides with the LGR formula
In this section, the LGL formula is derived from Gauss-Lobatto-Jacobi rule in the same limiting process as in Section 2 and is shown to coincide with the LGR formula.
Gautschi [13] gives the (n + 2)-point Gauss-Lobatto-Jacobi rule by where the nodes, 1 = ξ n,n+1 = −1 (note that the node sequence is opposite to the conventional one, namely, with ξ being the zeros of the Jacobi polynomial P (α+1,β+1) n (x) of degree n, and the weights ω Note that the weights ω 3) are expressed in a simpler form, However, we use an alternative expression for the weights ω where ω J(α+1,β+1) n,l are the weights for the GJ rule for the integral with the weight function (1 − x) α+1 (1 + x) β+1 , see (2.4).

Proof
The proof is done in a similar way to that for Lemma 3.1 except for the choice of g(x) = (1 − x 2 ) f (x), where f (x) is a polynomial of degree 2n − 1.In view of (4.2) and the fact that g(±1) = 0. we have 1 −1 ).
(4.6)Then, the rightmost-hand side of (4.6) is the n-point GJ rule with the weights ω . Hence, we have (4.5).
Applying the Gauss-Lobatto-Jacobi rule (4.1) to the integral I a n f in (2.7), we obtain the approximation where the nodes x LJ(α,0) n,l are given by and the weights w ) n,l where we used (4.5) for ω LJ(α,0) n,l (1 ≤ l ≤ n) and (2.10) for ω J(α+1,1) n,l .Theorem 4.2 below shows that the Gauss-Lobatto-Jacobi rule I LJ(α,0) n f in (4.7) goes to the LGR formula given by (1.4) as n → ∞.Therefore, the LGL formula coincides with the LGR formula.
with x LGR h,l given by (1.5) and w LGR h,l by (1.6).

Proof of Theorems 1.4 and 1.5
To prove Theorems 1.4 and 1.5, the following Lemma 5.1 is required with some definitions.For an integer parameter λ = 0 and 1, let the approximation Q (λ) f to the integral I f in (1.1) and its error E Q (λ) f be written by where the nodes x (λ) l and weights w (λ) l are given by x (1) x respectively, with the dummy node x (0) 0 = 0 for convenience.Then, in view of (1.3), (1.5) and (1.6), for h = π , we see that The errors are bounded by We prove Theorem 1.4 from Lemma 5.1 (λ = 0) and Theorem 1.5 from Lemma 5.1 (λ = 1) in Section 5.1.Lemma 5.1 is proved in Section 5.2.
Further properties of α (τ ) are shown in the Lemmas 5.2-5.5 below.
The second and third integrals are evaluated as follows.
Lemma 5.6 For A n in (5.38) and C n in (5.39), (5.47) Proof We verify (5.46).Similarly, we can verify (5.47) and its proof is therefore omitted.Let τ = −w n + it ∈ A n .In view of (1.11), for an arbitrary ε > 0, there exists a sufficient large N 1 such that for all n ≥ N 1 and any t , Therefore, in view of (5.15) and Lemma 5.5, for all n ≥ max{N 1 , N }, we have Hence, we have By letting ε → 0. the relation (5.46) is verified.

Conclusion
In this paper, we presented two formulae for approximating the semi-infinite integral (1.1), namely, formulae LGJ and LGR, by extending the Gauss-Jacobi and Gauss-Jacobi-Radau rules for an integral on the finite interval [−1, 1].Developing contour integration representations of the errors of the formulae, we showed that for a function f (z) analytic on the right region of a parabola, and its boundary ∂ P d in the complex plane z, the errors of both formulae converge in O(e −2π d/h ) as step size h → +0.This convergence rate is much faster than the rate O(h j e −π d/h ) ( j = 1, . . ., 5) of the extended formulae of the Clenshaw-Curtis-type rules [2].The average of the LGJ and LGR formulae is shown to converge faster in O(h e −2π d/h ).By numerical examples, we showed that the LGR formula behaves like an anti-LGJ formula and the average of the LGJ and LGR formulae gives better approximations.

. 10 )Definition 1 . 3 (
and we set arg z = θ (0 < θ < 2π ) for z ∈ ∂ P d .function family B(α, d)) For the domain P d given by (1.9) and its closure P d , define by B(α, d) the set of functions f (z) analytic on P d and satisfying the conditions max

Fig. 4
Fig. 4 Contour A n + B n,ε + C n + ∂ D d,n on the closure D d