Weighted quadrature formulas
for semi-infinite range integrals\(^\ast \)

Gradimir V. Milovanović\(^\S \)

November 6, 2015.

\(^\ast \) This paper was supported by the Serbian Ministry of Education, Science and Technological Development (No. #OI 174015).

\(^\S \) Serbian Academy of Sciences and Arts, Beograd, Serbia & State University of Novi Pazar, Serbia, e-mail: gvm@mi.sanu.ac.rs.

Dedicated to prof. Ion Păvăloiu on the occasion of his 75th anniversary

Weighted quadrature formulas on the half line \((a,+\infty )\), \(a{\gt}0\), for non-exponentially decreasing integrands are developed. Such \(n\)-point quadrature rules are exact for all functions of the form \(x\mapsto x^{-2}P(x^{-1})\), where \(P\) is an arbitrary algebraic polynomial of degree at most \(2n-1\). In particular, quadrature formulas with respect to the weight function \(x\mapsto w(x)=x^\beta \log ^m x\) (\(0\le \beta {\lt}1\), \(m\in \mathbb {N}_0\)) are considered and several numerical examples are included.

MSC. 65D30, 65D32

Keywords. Gaussian quadrature rules, nodes, Christoffel numbers, non-exponentially decreasing integrands.

1 Introduction

In this paper we consider weighted quadrature formulae on the half line \((a,+\infty )\),

\begin{equation} \label{GQCBS} \int _a^{+\infty }w(x) f(x){\, \mathrm{d}}x=\sum _{k=1}^n A_k f(x_k)+R_n(f), \end{equation}
1.1

where \(a\) is a finite real number and \(x\mapsto w(x)\) is a given weight function. Such a quadrature formula for \(a=0\) and \(w(x)=x^\alpha {\mathrm{e}}^{-x}\), \(\alpha {\gt}-1\), is the well known generalized Gauss-Laguerre quadrature rule (cf. [ 10 , p. 325 ] ), which is exact for all algebraic polynomials of degree at most \(2n-1\), i.e., when \(f\in {\EuScript P}_{2n-1}\).

Error analysis and convergence of such Gaussian formulas on unbounded intervals (with the classical measures of Laguerre and Hermite) was given in 1928 by Uspensky [ 17 ] . Otherwise, the corresponding problems for quadrature rules on finite intervals was studied much earlier by [ 15 ] , Markov [ 8 ] , Stieltjes [ 16 ] , etc. On some new results in this directions see books [ 5 ] and [ 10 ] , including the so-called truncated quadrature rules obtained by ignoring the last part of its nodes (see Mastroianni and Monegato [ 9 ] ).

Very recently Gautschi [ 6 ] has constructed a special logarithmically weighted quadrature formula on \((0,+\infty )\), when \(x\mapsto (x-1-\log x)\mathrm{e}^{-x}\). Also, Xu and Milovanović [ 18 ] have developed generalized Gaussian quadrature rules of the form \((\ref{GQCBS})\), with \(x\mapsto w(x)=\mathrm{e}^{-x}\) on \((0,+\infty )\), which are exact on the set of basis functions \(\{ 1,\log x,x,x\log x,\ldots ,x^{n-1},x^{n-1}\log x\} \). In the other words, these rules are exact for each \(f(x)=p(x)+q(x)\log x\), where \(p,q\in {\EuScript P}_{n-1}\), so that they can calculate integrals with a sufficient accuracy, regardless of whether their integrands contain a logarithmic singularity, or they do not. For a similar approach for integrals on the finite intervals see [ 12 ] and [ 14 ] .

On the other side, a large number of integrals of the form \(\int _a^{+\infty }F(x){\, \mathrm{d}}x\) which appear in applications do not have exponentially decreasing integrands \(F(x)\), and in such cases Gauss-Laguerre quadrature rules are notoriously poor (see Evans [ 2 ] ). As a starting simple example, Evans [ 2 ] has considered \(F(x)=1/(x^2+0.25)\), where the convergence of the corresponding integral depends on the \(1/x^2\) term for large \(x\). He has proposed a quadrature method based on the set of basis functions \(\{ 1/x^k\} \) and demonstrated its effectiveness on a series of numerical examples.

In this paper we develop a general approach for constructing a class of \(n\)-point generalized quadrature rules \((\ref{GQCBS})\) of Gaussian type on \((a,+\infty )\), \(a{\gt}0\), which are exact for all functions of the form \(x\mapsto x^{-2}P(x^{-1})\), where \(P\) is an arbitrary algebraic polynomial of degree at most \(2n-1\). In particular, we consider quadrature formulas with respect to the weight function \(x\mapsto w(x)=x^\beta \log ^m x\) (\(0\le \beta {\lt}1\), \(m\in \mathbb {N}_0\)), which reduces to the constant weight for \(\beta =0\) and \(m=0\). In order to show the efficiency of the obtained quadrature rules we present a few numerical examples.

2 Generalized Weighted Gaussian Rules

Suppose \(a{\gt}0\), as well as that the weight function \(x\mapsto w(x)\) on \((a,+\infty )\) is such that

\begin{equation} \label{uslovi} 0<\int _a^{+\infty }\frac{w(x)}{x^2}{\, \mathrm{d}}x<+\infty . \end{equation}
2.1

Following Evans [ 2 ] , we develop a general approach for constructing generalized Gaussian quadrature formulas of the form \((\ref{GQCBS})\). In the cases of integrals on \((\alpha , +\infty )\), when \(\alpha {\lt} a\), we simply take

\[ \int _\alpha ^{+\infty }w(x)f(x){\, \mathrm{d}}x=\int _\alpha ^{a}w(x)f(x){\, \mathrm{d}}x+ \int _a^{+\infty }w(x)f(x){\, \mathrm{d}}x \]

and apply to first integral on the right hand side some of rules for calculating integrals on the finite intervals. Also, we mention here that a faster convergence of the corresponding quadrature process can be achieved by taking a greater value of \(a\).

Thus, the basic idea is to construct a quadrature formula of the form \((\ref{GQCBS})\), which is exact for all functions of the form

\[ x\mapsto \frac1{x^2}P_m\Big(\frac1{x}\Big),\quad m=0,1,\ldots ,2n-1, \]

where \(P_m(t)\) are arbitrary selected algebraic polynomials in \(t\) of degree \(m\), i.e.,

\begin{equation} \label{SystemEq} \int _a^{+\infty }w(x) \frac1{x^2}P_m\Big(\frac1{x}\Big){\, \mathrm{d}}x=\sum _{k=1}^n \frac{A_k}{x_k^2}P_m\Big(\frac1{x_k}\Big), \quad m=0,1,\ldots ,2n-1. \end{equation}
2.2

Remark

Because of linearity, it is easy to see that this system of \(2n\) nonlinear equations in \(x_k\) and \(A_k\), \(k=1,\ldots ,n\), is equivalent to the corresponding system with monomials, i.e., when \(P_m(x)=x^m\), \(m=0,1,\ldots ,2n-1\).â–¡

On the other side we consider the Gauss-Christoffel quadrature formula with respect to the weight function \(t\mapsto w(1/t)\) on \((0,1/a)\), i.e.,

\begin{equation} \label{IntegFin} \int _0^{1/a}w\Big(\frac1{t}\Big)g(t){\, \mathrm{d}}t= \sum _{k=1}^n B_k g(\tau _k)+R_n^G(g), \end{equation}
2.3

where \(\tau _k\) and \(B_k\) are its nodes and Christoffel numbers, respectively, and \(R_n^G(g)\) is the corresponding remainder term. According to \((\ref{uslovi})\), such quadrature formulas exist uniquely, because the all moments \(\mu _k=\int _0^{1/a} w(1/t)t^k{\, \mathrm{d}}t\), \(k\ge 0\), exist and \(\mu _0{\gt}0\).

It is known that the nodes \(\tau _k\) in \((\ref{IntegFin})\) are eigenvalues of the following symmetric tridiagonal Jacobi matrix (cf. [ 10 , pp. 325–328 ] )

\begin{equation} \label{Jac-matr} J_n=\left[\begin{array}{ccccc} \alpha _0& \sqrt{\beta _1}& & & {\mathbf{O}}\\[2mm] \sqrt{\beta _1}& \alpha _1& \sqrt{\beta _2}\\[2mm]& \sqrt{\beta _2}& \alpha _2& \ddots \\[2mm]& & \ddots & \ddots & \sqrt{\beta _{n-1}}\\[2mm] {\mathbf{O}}& & & \sqrt{\beta _{n-1}}& \alpha _{n-1} \end{array}\right], \end{equation}
2.4

where \(\alpha _k\) and \(\beta _k\) are coefficients in the three-term recurrence relation

\begin{eqnarray} \label{TTRR} \pi _{k+1}(t)& \! \! \! =\! \! \! & (t-\alpha _k)\pi _k(t)-\beta _k\pi _{k-1}(t),\quad k=0,1,\ldots ,\\[2mm] \pi _0(t)& \! \! \! =\! \! \! & 1,\ \ \pi _{-1}(t)=0,\nonumber \end{eqnarray}

for the (monic) polynomials \(\{ \pi _k\} _{k\in \mathbb {N}_0}\) orthogonal with respect to the inner product

\begin{equation} \label{inner} (p,q)=\int _0^{1/a}w\Big(\frac1{t}\Big)p(t)q(t){\, \mathrm{d}}t. \end{equation}
2.6

In fact, \(\pi _n(t)=(t-\tau _1)\cdots (t-\tau _n)\).

The weight coefficients \(B_k\) in (2.3) are given by

\[ B_k=\beta _0 v_{k,1}^2,\quad k=1,\ldots ,n, \]

where \(v_{k,1}\) is the first component of the eigenvector \({\mathbf{v}}_k\ (=[v_{k,1}\ \ldots \ v_{k,n}]^{\mathrm T})\) corresponding to the eigenvalue \(\tau _k\) and normalized such that \({\mathbf{v}}_k^{\mathrm T}{\mathbf{v}}_k=1\), and \(\beta _0=\mu _0=\int _0^{1/a}w(1/t){\, \mathrm{d}}t\).

The most popular method for solving this eigenvalue problem is the Golub-Welsch procedure, obtained by a simplification of QR algorithm [ 7 ] . This procedure is implemented in several packages including the most known ORTPOL given by Gautschi [ 4 ] .

As we can see from \((\ref{Jac-matr})\), for constructing Gauss–Christoffel quadratures \((\ref{IntegFin})\) for any number of nodes less than or equal to \(n\), we need the first \(n\) recursion coefficients \(\alpha _k\) and \(\beta _k\), \(k=0,1,\ldots ,n-1\), in \((\ref{TTRR})\).

In general, the recursion coefficients are known explicitly only for some narrow classes of orthogonal polynomialsc(e.g. for the classical orthogonal polynomials). In the case of the so-called strongly non-classical polynomials, these recursion coefficients must be constructed numerically (cf. [ 3 ] , [ 5 ] , [ 10 , pp. 159–166 ] ). However, recent progress in symbolic computation and variable-precision arithmetic today makes it possible to generate the recursive coefficients in \((\ref{TTRR})\) directly by using the original Chebyshev method of moments. Respectively symbolic/variable-precision software for orthogonal polynomials and Gaussian (and similar) quadratures is available. Our Mathematica package OrthogonalPolynomials (see [ 1 ] and [ 13 ] ), is downloadable from the web site http://www.mi.sanu.ac.rs/~gvm/. Also, there is Gautschi’s software in Matlab (packages OPQ and SOPQ).

Now, we can give our main result:

Theorem 2.1

Let \(x\mapsto w(x)\) be a weight function on \((a,+\infty )\), \(a{\gt}0\), such that the condition \((\ref{uslovi})\) holds. Assume also that \(\tau _k\) and \(B_k\), \(k=1,\ldots ,n\), are nodes and Christoffel numbers of the Gaussian quadrature formula \((\ref{IntegFin})\), respectively. Then there exists the generalized Gaussian quadrature formula

\begin{equation} \label{novaQ} \int _a^{+\infty }w(x) f(x){\, \mathrm{d}}x=\sum _{k=1}^n A_k f(x_k)+R_n(f), \end{equation}
2.7

with

\begin{equation} \label{PARAM} x_k=\frac{1}{\tau _k}, \quad A_k=\frac{B_k}{\tau _k^2}>0,\quad k=1,\ldots ,n, \end{equation}
2.8

which is exact for all functions of the form \(f(x)=x^{-2}P(x^{-1})\), where \(P\in {\EuScript P}_{2n-1}\).

The remainder term in this quadrature rule can be expressed in the following form \(R_n(f)=R_n^G(g)\), where \(g(t)=t^{-2}f(t^{-1})\).

Proof â–¼
We start with the system of \(2n\) nonlinear equations \((\ref{SystemEq})\), whose solution determines the parameters of the quadrature formula \((\ref{novaQ})\). Our aim is to prove that this solution uniquely exists.

First, we take the sequence of orthogonal polynomials \(\{ \pi _m\} _{m=0}^{2n-1}\) in the system \((\ref{SystemEq})\) and then by a simple change of variables \(x=1/t\) in the integral on the left hand side we obtain

\begin{eqnarray*} \int _a^{+\infty }w(x) \frac1{x^2}\pi _m\Big(\frac1{x}\Big){\, \mathrm{d}}x& =& \int _0^{1/a}w\Big(\frac1{t}\Big)\pi _m(t){\, \mathrm{d}}t\\[1mm]& =& (\pi _0,\pi _m)\\[1mm]& =& \mu _0\delta _{0,m}, \end{eqnarray*}

where the inner product is defined by \((\ref{inner})\) and \(\delta _{k,m}\) is Kronecker’s delta.

Evidently, this leads to the system of equations

\begin{equation} \label{sys1} \sum _{k=1}^n A_k \frac{1}{x_k^2}\pi _m\Big(\frac1{x_k}\Big)= \mu _0\delta _{0,m},\quad m=0,1,\ldots ,2n-1, \end{equation}
2.9

but, by an application of the Gaussian rule \((\ref{IntegFin})\), it gives also another system of equations

\begin{equation} \label{sys2} \sum _{k=1}^n B_k \pi _m(\tau _k)=\mu _0\delta _{0,m},\quad m=0,1,\ldots ,2n-1, \end{equation}
2.10

because \(R_n^G(g)=0\) for each \(g\in {\EuScript P}_{2n-1}\). The last system has the unique solution, and it represents the parameters \(\tau _k\) and \(B_k\), \(k=1,\ldots ,n\), of the Gaussian quadrature \((\ref{IntegFin})\).

Since the systems of equations \((\ref{sys1})\) and \((\ref{sys2})\) are equivalent, the statement of this theorem follows directly.

Proof â–¼

3 Special cases and numerical examples

In this section we consider special cases of quadrature formulas with respect to the weight function \(x\mapsto w(x)=x^\beta \log ^m x\), where \(0\le \beta {\lt}1\) and \(m\in \mathbb {N}_0\). For \(\beta =0\) and \(m=0\), it reduces to the constant weight \(w(x)=1\). In order to show the efficiency of the obtained quadrature formulas we present a few numerical examples.

We start this section with the weight function \(x\mapsto w(x)=x^\beta \), \(0\le \beta {\lt}1\).

The condition \((\ref{uslovi})\) is satisfied, because

\[ \int _a^{+\infty } \frac{w(x)}{x^2}{\, \mathrm{d}}x=\frac{a^{\beta -1}}{1-\beta }. \]

Here we consider only the case \(\beta =0\), i.e., when \(w(x)=1\). Since

\[ \int _a^{+\infty }f(x){\, \mathrm{d}}x=a\int _1^{+\infty }f(ax){\, \mathrm{d}}x, \]

we see that for this important case the following statement holds.

Proposition 3.1

Let

\begin{equation} \label{const} \int _a^{+\infty }f(x){\, \mathrm{d}}x=\sum _{k=1}^n A_k(a)f(x_k(a))+R_n(f;a),\quad a>0, \end{equation}
3.1

be a generalized Gaussian quadrature \((\ref{novaQ})\) \((\)with the constant weight function \(w(x)=1)\). Then

\[ A_k(a)=aA_k(1)\quad \mbox{and}\quad x_k(a)=ax_k(1),\ \ k=1,\ldots ,n. \]

This means that it is enough to know only quadrature parameters for \(a=1\). These parameters can be obtained directly using \((\ref{PARAM})\) and Gauss-Legendre parameters \(\tau _k\) and \(B_k\) for transformed interval \((0,1)\).

Recursive coefficients in \((\ref{TTRR})\), in this case for translated monic Legendre polynomials, are

\[ \alpha _k=\frac12,\ \ k\ge 0,\quad \beta _0=1,\ \ \beta _k=\frac{k^2}{4(4k^2-1)},\ \ k\ge 1. \]

Otherwise, it can be obtained using our Mathematica Package OrthogonalPolynomials in symbolic form (see [ 1 ] and [ 13 ] ). For example, if we need the first forty recurrence coefficients, then we start with the first eighty moments \(\mu _k=1/(k+1)\), \(k=0,1,\ldots ,79\), and then we use the standard Chebyshev algorithm (cf. [ 10 , 160–162 ] : 

  << orthogonalPolynomials`
  mom=Table[1/(k+1), {k,0,79}];
  {al,be} = aChebyshevAlgorithm[mom, Algorithm -> Symbolic]
These recursive coefficients enable us to construct quadrature formulas \((\ref{IntegFin})\) for any number of nodes up to \(40\).

However, in this Legendre case (translated to \((0,1)\)) we can directly use aGaussianNodesWeights routine to construct nodes and weights in the Gaussian quadrature formula \((\ref{IntegFin})\), as well as ones in the quadrature formula \((\ref{novaQ})\): 

  << orthogonalPolynomials`
  transLeg[n_] := (aGaussianNodesWeights[n, {aLegendre},
     WorkingPrecision -> 70, Precision -> 65] + {1,0})/2;
     parQF = Table[transLeg[n], {n,2,40,2}];
     For[m = 1, m < 21, m++,
        parQF[[m]][[2]] = parQF[[m]][[2]]/parQF[[m]][[1]]^2;
        parQF[[m]][[1]] = 1/parQF[[m]][[1]];]
Thus, in this way for \(a=1\), we obtain quadrature parameters \(x_k\) and \(A_k\) for each \(n=2(2)40\).

Example 3.2

In order to show the efficiency of our quadrature rule \((\ref{novaQ})\) we apply it to the integral

\[ J(a;c)=\int _a^{+\infty }\frac{1}{(x-2)^2+c^2}{\, \mathrm{d}}x =\frac{1}{2 c}\left[\pi -2 \arctan \big(\frac{a-2}{c}\big)\right], \]

for different values of \(a{\gt}0\) and \(c{\gt}0\). In Figure 3.1 we present graphics of the function

\[ x\mapsto f(x;c)=\frac1{(x-2)^2+c^2} \]

for \(c=\frac18,\frac14,\frac12\), and \(1\), as well as the corresponding graphics of the exact values of this integral \(J(a;c)\) (right).

In order to test the quadrature formula \((\ref{const})\), we apply it to \(J(a;1)\) for \(a=\frac12\), \(1\), \(2\), \(3\), \(4\), and \(8\), when \(n=2(2)40\).

\[ \includegraphics[width=0.50\textwidth ]{GraffunBOJA.pdf}\includegraphics[width=0.50\textwidth ]{TVRintBOJA.pdf} \]
Figure 3.1 The function \(x\mapsto f(x;c)\) (left) and the integral \(a\mapsto J(a;c)\) (right) for \(c=\frac18\) (blue line), \(c=\frac14\) (black line), \(c=\frac12\) (brown line), and \(c=1\) (red line).

Relative errors in the quadrature sums

\[ Q_n(f(\cdot ;c);a)=\sum _{k=1}^n A_k(a)f(x_k(a);c), \]

defined by

\[ \operatorname {err}_n(f(\cdot ;c);a)=\Bigl|\frac{Q_n(f(\cdot ;c);a)-J(a;c)}{J(a;c)}\Bigr|, \]

are displayed in Figure 3.2 in a log-scale.

\[ \includegraphics[width=0.80\textwidth ]{err_graf1.pdf} \]
Figure 3.2 Relative errors \(\operatorname {err}_n(f(\cdot ;1);a)\) in quadrature sums \(Q_n(f(\cdot ;1);a)\).

Numerical results show that the convergence is much faster if the parameter \(a\) is larger. For example, if \(a = 2\), then for \(n=10(10)40\), the relative errors are \(1.71\times 10^{-7}\), \(1.83\times 10^{-14}\), \(1.91\times 10^{-21}\), \(1.94\times 10^{-28}\), respectively, while the corresponding errors for \(a=4\) are \(5.52\times 10^{-15}\), \(1.21\times 10^{-29}\), \(1.40\times 10^{-44}\), \(1.44\times 10^{-59}\).

Otherwise, this integrand \(f(x;c)\) has poles at the points \(2\pm {\mathrm{i}}c\), which are approaching the real line when \(c\) tends to zero. In this case, for small values of \(a\) (near \(2\) or less than \(2\)), the convergence of the quadrature process slows down considerably, because of a strong influence of these singularities. This effect can be seen from Table 3.1, where quadrature approximations and corresponding relative errors are presented for \((a,c)=(1,\frac14)\), \((\frac{21}{10},10^{-6})\), and \((4,10^{-6})\). In order to save space, in last case only relative errors are given. Digits in error are underlined, and numbers in parenthesis indicate the decimal exponents.

Notice that the integral \(J(a;0)\) for \(a\le 2\) does not exist.

Finally, the last column shows that for \((a,c)=(4,10^{-6})\), the convergence of the quadrature rule \((\ref{const})\) is very fast.

\(n\)

\((a,c)=(1,1/4)\)

\((a,c)=(21/10,10^{-6})\)

\((a,c)=(4,10^{-6})\)

\(2\)

\(\underline{2.83088}\ \)

\(7.56(-1)\)

\(\underline{4.21706255691703}\)

\(5.78(-1)\)

\(5.92(-3)\ \, \)

\(4\)

\(\underline{5.38719}\ \)

\(5.35(-1)\)

\(\underline{8.01223217799471}\)

\(1.99(-1)\)

\(9.70(-6)\ \, \)

\(6\)

\(\underline{7.41379}\ \)

\(3.60(-1)\)

\(9.\underline{47887835712778}\)

\(5.21(-2)\)

\(1.24(-8)\ \, \)

\(8\)

\(\underline{8.88711}\ \)

\(2.33(-1)\)

\(9.\underline{88043864297441}\)

\(1.20(-2)\)

\(1.42(-11)\)

\(10\)

\(\underline{9.89102}\ \)

\(1.46(-1)\)

\(9.9\underline{7447558340612}\)

\(2.55(-3)\)

\(1.53(-14)\)

\(20\)

\(11.\underline{45438}\ \)

\(1.14(-2)\)

\(9.99999\underline{276505451}\)

\(7.23(-7)\)

\(1.47(-29)\)

\(30\)

\(11.5\underline{7808}\ \)

\(7.23(-4)\)

\(9.99999999\underline{813998}\)

\(1.53(-10)\)

\(1.08(-44)\)

\(40\)

\(11.586\underline{06}\ \)

\(3.41(-5)\)

\(9.999999999666\underline{38}\)

\(2.86(-14)\)

\(6.99(-60)\)

Table 3.1 Quadrature sums \(Q_n(f(\cdot ;c);a)\) and their relative errors         
\(err_n(f(\cdot ;c);a)\) for integrals \(J(a;c)\).

In the sequel we consider quadrature rules with respect to the weight function \(x\mapsto w(x)=x^\beta \log x\), \(0\le \beta {\lt}1\). Here we suppose that \(a\ge 1\). The condition \((\ref{uslovi})\) is satisfied, because

\begin{equation} \label{U1} 0<\int _a^{+\infty }\frac{w(x)}{x^2}{\, \mathrm{d}}x= \frac{a^{\beta -1}}{(1-\beta )^2}\left[1+(1-\beta )\log a\right]. \end{equation}
3.2

In this case, the moments

\[ \mu _k=\int _0^{1/a}w(1/t)t^k{\, \mathrm{d}}t=\int _0^{1/a}t^{k-\beta }\log \frac1t{\, \mathrm{d}}t \]

can be expressed in the form

\begin{equation} \label{Mom1} \mu _k=\frac{a^{\beta -k-1} [(k+1-\beta )\log a+1]}{(k+1-\beta )^2},\quad k\ge 0. \end{equation}
3.3

Taking the first one hundred moments (mom) (e.g. for \(a=1\) and \(\beta =1/4\)) and using Mathematica Package OrthogonalPolynomials, we can get the first fifty recurrence coefficients \(\alpha _k\) and \(\beta _k\) (denoted by {al1,be1}) in the three-term recurrence relation (2.5) in a symbolic form 

  << orthogonalPolynomials`
  mom=Table[(a^(-1+b-k)(1+(1-b+k)Log[a]))/(1-b+k)^2, {k,0,99}];
  mom1=mom/. {a->1, b->1/4}
  {al1,be1} = aChebyshevAlgorithm[mom, Algorithm -> Symbolic]
For example, first four coefficients are
\[ \alpha _0=\frac{9}{49},\, \alpha _1=\frac{209897}{452025},\, \alpha _2=\frac{6582284926939}{13538179995075},\, \alpha _3=\frac{7618613698603068100869609}{15464687102113919816429449} \]
and
\[ \beta _0=\frac{16}{9},\ \, \beta _1=\frac{11808}{290521},\ \, \beta _2=\frac{213147564896}{3717280400625},\ \, \beta _3=\frac{421267942813254097088}{6997413354065613077481}. \]

Example 3.3

As a test example we consider the function

\[ x\mapsto f(x)=\frac1{(x+1)^{2}} \]

and integral (see [ 2 ] )

\begin{align*} I(f;a)=& \int _a^{+\infty }\frac{x^{1/4}\log x}{(x+1)^2}{\, \mathrm{d}}x\\ =& \frac{1}{36 a^{\frac74} (a+1)} \left\{ -9 (a+1)\Phi \Bigl(-\frac{1}{a},2,\frac{7}{4}\Bigr)\right.\qquad \quad \\[2mm]& \qquad \left.+4 a \left[3 (a+1)(\log a+4) \, _2F_1\Bigl(\frac{3}{4},1;\frac{7}{4};-\frac{1}{a}\Bigr) +4 a+9 a \log a+4\right]\right\} , \end{align*}

where \(\Phi \) and \(_2F_1\) are the Lerch transcendent and Gauss hypergeometric function, defined by

\[ \Phi (z,s,a)=\sum _{k=0}^{+\infty } \frac{z^k}{(k+a)^s}\quad \mbox{and}\quad _2F_1(a,b;c;z)=\sum _{k=0}^{+\infty }\frac{(a)_k(b)_k}{(c)_k}\frac{z^k}{k!}, \]

respectively, and \((a)_k=a(a+1)\ldots (a+k-1)\) is the Pochhammer symbol.

We consider this integral for two values of the lower bound: \(a=1\) and \(a={\mathrm{e}}\), i.e.,

\[ I(f;1)={1.35974328097600895}\ldots \ \ \mbox{and}\ I(f;\mathrm{e})={1.22897618668037255}\ldots \ . \]
Applying Gauss-Laguerre rule to \(I(f;1)\) (translated from \((1,+\infty )\) to \((1,+\infty )\)) gives poor results. Relative errors in the corresponding Gauss-Laguerre quadrature sums are presented in Table 3.2. As we can see only two two decimal digits are true in quadrature sum with \(2048\) nodes!

\(n=2\)

\(n=8\)

\(n=32\)

\(n=128\)

\(n=512\)

\(n=2048\)

\(6.72(-1)\)

\(3.60(-1)\)

\(1.64(-1)\)

\(7.00(-2)\)

\(2.90(-2)\)

\(1.18(-2)\)

Table 3.2 Relative errors in Gauss-Laguerre quadrature sums         
with \(n=2,8,32,128,512\) and \(2048\) nodes.

Now, we apply our quadrature formula \((\ref{novaQ})\), with parameters given by \((\ref{PARAM})\), to \(I(f;1)\) and \(I(f;\mathrm{e})\), with only \(n=2(2)12\) nodes. The relative errors in the quadrature sums \(Q_n(f;a)=\sum _{k=1}^n A_kf(x_k)\),

\[ \operatorname {err}_n(f;a)=\Bigl|\frac{Q_{n}(f;a)-I(f;a)}{I(f;a)}\Bigr|, \]

are presented in Table 3.3.

\(a\)

\(n=2\)

\(n=4\)

\(n=6\)

\(n=8\)

\(n=10\)

\(n=12\)

\(1\)

\(2.94(-3)\)

\(4.24(-6)\ \, \)

\(5.15(-9)\ \)

\(5.72(-12)\)

\(4.74(-13)\)

\(7.07(-13)\)

\({\mathrm{e}}\)

\(2.40(-4)\)

\(1.64(-8)\ \, \)

\(\ 8.91(-13)\)

\(8.83(-14)\)

\(5.31(-14)\)

\(3.80(-14)\)

\({\mathrm{e}}^2\)

\(7.18(-6)\)

\(1.28(-11)\)

\(\ 3.10(-14)\)

     
Table 3.3 Relative errors \(\operatorname {err}_n(f;a)\) in quadrature sums \(Q_n(f;a)\)         
for different number of nodes \(n\) and three values of \(a\) \((=1,{\mathrm{e}},\ \mbox{and}\ {\mathrm{e}}^2)\).

As we can see, the convergence is faster when \(a\) is bigger. In the third line of the same table we also present the corresponding relative errors when \(a={\mathrm{e}}^2\) and \(n=2,4\), and \(6\).

Finally, we mention that this approach can be applied also in the case of the weight functions

\[ w(x)=w_m(x)=x^\beta \log ^m x,\quad 0\le \beta {\lt}1,\ m=2,3,\ldots \ , \]

on the interval \((a,+\infty )\), with \(a\ge 1\).

The condition \((\ref{uslovi})\) for \(U_m=\int _a^{+\infty }x^{-2}w_m(x){\, \mathrm{d}}x\) is also satisfied, because

\[ U_m=\frac1{1-\beta }\left[mU_{m-1}+a^{\beta -1}\log ^m a\right],\quad m=2,3,\ldots \, , \]

where \(U_1\) is given in \((\ref{U1})\). The corresponding moments

\[ \mu _k^{[m]}=\int _0^{1/a}t^{k-\beta }\log ^m\frac1{t}{\, \mathrm{d}}t,\quad k=0,1,\ldots \, , \]

can be expressed recursively in terms of the moments \(\mu _k^{[m-1]}\),

\[ \mu _k^{[m]}=\frac{1}{k+1-\beta }\left(m \mu _k^{[m-1]}+a^{\beta -k-1}\log ^m a\right), \quad m=2,3,\ldots \, , \]

where the moments \(\mu _k^{[1]}\ (\equiv \mu _k)\) are given by \((\ref{Mom1})\). For example, for \(m=2\) we get

\[ \mu _k^{[2]}=\frac{a^{\beta -k-1} \left[(k+1-\beta )^2\log ^2a +2(k+1-\beta ) \log a+2\right]}{(k+1-\beta )^3},\quad k\ge 0. \]

The corresponding recursive coefficients, for example for \(\beta =0\) and \(a=1\), are

\[ \alpha _0=\frac{1}{8},\, \alpha _1=\frac{115}{296},\, \alpha _2=\frac{28200187}{62721512},\, \alpha _3=\frac{28003451041760695}{59414538084233528}, \ \ldots \]
and
\[ \beta _0=2,\, \beta _1=\frac{37}{1728},\ \, \beta _2=\frac{211897}{4620375},\ \, \beta _3=\frac{945381680572419}{17600932734728000},\ \ldots . \]

Example 3.4

We consider the function \(x\mapsto 1/(1+x^2)\) and the corresponding weighted integral over \((a,+\infty )\)

\[ I(f;a)=\int _a^{+\infty }\frac{\log ^2x}{1+x^2}{\, \mathrm{d}}x, \]

for two values of \(a\) \((=1\) and \(={\mathrm{e}})\), for which

\[ I(f;1)={1.93789229251873876096726969169}\ldots \]
and
\[ \ \ I(f;{\mathrm{e}})={1.80988687939786942602016447246}\ldots \ . \]

Applying our quadrature formula \((\ref{novaQ})\), with parameters given by \((\ref{PARAM})\), to \(I(f;1)\) and \(I(f;\mathrm{e})\), with \(n=2(2)12\) nodes, we get the corresponding quadrature approximations \(Q_n(f;a)\) with the relative errors \(\operatorname {err}_n(f;a)\) presented in Table 3.4.

\(a\)

\(n=2\)

\(n=4\)

\(n=6\)

\(n=8\)

\(n=10\)

\(n=12\)

\(1\)

\(1.66(-4)\)

\(1.31(-6)\ \ \)

\(1.98(-10)\, \)

\(5.73(-12)\)

\(2.08(-15)\)

\(2.56(-17)\)

\({\mathrm{e}}\)

\(5.33(-5)\)

\(5.04(-10)\)

\(1.86(-13)\)

\(2.05(-17)\)

\(1.22(-21)\)

\(3.30(-26)\)

Table 3.4 Relative errors \(\operatorname {err}_n(f;a)\) in quadrature sums \(Q_n(f;a)\)       
for different number of nodes \(n\) and two values of \(a\) \((=1\) and \(={\mathrm{e}})\).

Here also we can note a faster convergence when \(a\) has a larger value.

Bibliography

1

A.S. Cvetković and G.V. Milovanović, The Mathematica Package “OrthogonalPolynomials”, Facta Univ. Ser. Math. Inform., 19, pp. 17–36, 2004.

2

G.A. Evans, Some new thoughts on Gauss-Laguerre quadrature, Int. J. Comput. Math. 82, pp. 721–730, 2005. \includegraphics[scale=0.1]{ext-link.png}

3

W. Gautschi, On generating orthogonal polynomials, SIAM J. Sci. Statist. Comput., 3, pp. 289–317, 1982. \includegraphics[scale=0.1]{ext-link.png}

4

W. Gautschi, Algorithm 726: ORTHPOL – A package of routines for generating orthogonal polynomials and Gauss-type quadrature rules, ACM Trans. Math. Software, 20, pp. 21–62, 1994. \includegraphics[scale=0.1]{ext-link.png}

5

W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Clarendon Press, Oxford, 2004.

6

W. Gautschi, Gauss quadrature routines for two classes of logarithmic weight functions, Numer. Algor. 55, pp. 265–277, 2010. \includegraphics[scale=0.1]{ext-link.png}

7

G. Golub and J.H. Welsch, Calculation of Gauss quadrature rules, Math. Comp. 23, pp. 221–230, 1969. \includegraphics[scale=0.1]{ext-link.png}

8

A. Markov, Sur la méthode de Gauss pour le calcul approché des intégrales, Math. Ann., 25, pp. 427–432, 1885. \includegraphics[scale=0.1]{ext-link.png}

9

G. Mastroianni and G. Monegato, Truncated quadrature rules over \((0, \infty )\) and Nyström-type methods, SIAM J. Numer. Anal. 41, pp. 1870–1892, 2003. \includegraphics[scale=0.1]{ext-link.png}

10

G. Mastroianni and G.V. Milovanović, Interpolation Processes – Basic Theory and Applications, Springer-Verlag, Berlin – Heidelberg – New York, 2008.

11

G.V. Milovanović, Construction and applications of Gaussian quadratures with nonclassical and exotic weight function, Stud. Univ. Babeş-Bolyai Math., 60, pp. 211–233, 2015.

12

G.V. Milovanović, Generalized Gaussian quadratures for integrals with logarithmic singularity, FILOMAT (to appear).

13

G.V. Milovanović and A.S. Cvetković, Special classes of orthogonal polynomials and corresponding quadratures of Gaussian type, Math. Balkanica, 26, pp. 169–184, 2012.

14

G.V. Milovanović, T.S. Igić and D. Turnić, Generalized quadrature rules of Gaussian type for numerical evaluation of singular integrals, J. Comput. Appl. Math. 278, pp. 306–25, 2015. \includegraphics[scale=0.1]{ext-link.png}

15

C. Posse, Sur les quadratures, Nouv. Ann. Math. (2) 14, pp. 49–62, 1875.

16

T.J. Stieltjes, Quelques recherches sur la théorie des quadratures dites mécaniques, Ann. Sci. Éc. Norm. Paris, Sér. 2, 1, pp. 409–426, 1884.

17

J.V. Uspensky, On the convergence of quadrature formulas related to an infinite interval, Trans. Amer. Math. Soc., 30, pp. 542–559, 1928. \includegraphics[scale=0.1]{ext-link.png}

18

Zhenhua Xu and G.V. Milovanović, Efficient method for the computation of oscillatory Bessel transform and Bessel Hilbert transform (submitted).