\documentclass[11pt]{amsart}        %! the draft option shows (by filled rectangles)
                                          %! margin overflows. It disables normal output of figures
                                          %! One should suppress this option after the manuscript
                                          %! is in final form (no overflows occur)

\usepackage{amsmath,amsfonts,amssymb}
\usepackage{jnaat-v4,hyperref}
\usepackage{latexsym,euscript}
\usepackage{graphicx}                     %! used only when dealing with figures
\usepackage{subfigure}                    %! used only when dealing with subfigures
\usepackage{epstopdf}
\usepackage{numprint}                     %! used only when dealing with figures
\npdecimalsign{.}
\usepackage[normalem]{ulem}
\setcounter{MaxMatrixCols}{30}%


%\numberwithin{equation}{section}          %! you may use these two commands
%\newtheorem{theorem}{Theorem}[section]    %! to obtain equations numbered within sections
                                           %! This is recommended for long manuscripts
                                           %! Variant 1 (when not using the numinsec package)
\newtheorem{theorem}{Theorem}
                                           %! End of variant 1

                                           %! Variant 2 (when using the numinsec package)
%\newtheorem{theorem}{Theorem}[section]    %! If the case, please activate this variant, and delete the previous
                                           %! This is recommended for long manuscripts, with numerous sections and results
                                           %! End of variant 2

\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\newtheorem*{unremark}{Remark}             %! unnumbered remark
\newtheorem{example}[theorem]{Example}
\newtheorem*{acknowledgement}{Acknowledgement} %! unnumbered acknowledgement
%...                                       %! other commands may be defined, such as algorithm, etc.


\newcommand{\re}{\mathop{\mathrm{Re}}}
\newcommand{\im}{\mathop{\mathrm{Im}}}
\newcommand{\DS}{\displaystyle}
 \def\I{\mathrm{i}}

\def\Dj{\hbox{D\kern-.73em\raise.30ex\hbox{-}
\raise-.30ex\hbox{}}}
\def\dj{\hbox{d\kern-.33em\raise.80ex\hbox{-}
\raise-.80ex\hbox{\kern-.40em}}}

\def\Cale#1{\EuScript #1}
\newcommand{\PP}{{\Cale P}}
\newcommand{\op}{{\Cale A}}

\newcommand{\RES}{\mathop{\mathrm{Res}}\limits}
\newcommand{\loc}{\mathop{\mathrm{loc}}}
\newcommand{\lineal}{\mathop{\mathrm{span}}}
\newcommand{\Lip}{\mathop{\mathrm{Lip}}}

\def\RR{\mathbb{R}}
\def\NN{\mathbb{N}}
\def\CC{\mathbb{C}}
\def\ZZ{\mathbb{Z}}

\def\Res{{\mathrm Res}}

\newcommand{\E}{\mathrm{e}}
\newcommand{\D}{\,\mathrm{d}}



\def\bin#1#2{\Bigl({#1\atop#2}\Bigr)}

 \newcommand{\DD}{{\Cale{D}}}

\def\al{\alpha}
\def\be{\beta}
\def\de{\delta}
\def\ga{\gamma}
\def\la{\lambda}
\def\om{\omega}
\def\vr{\varrho}
\def\eps{\varepsilon}
\def\VGG{\mit\Gamma}
\def\GG{\Gamma}
\def\si{\sigma}
\def\supp{\mbox{\rm supp\,}}
\def\norm#1{{\Vert{#1}\Vert}}
\def\bnorm#1{{\bigm\Vert{#1}\bigm\Vert}}
\def\Bnorm#1{{\Bigm\Vert{#1}\Bigm\Vert}}
\def\dfrac#1#2{{\displaystyle {{#1}}\over{\displaystyle{#2}}}}
\def\dg{\mbox{\rm dg\,}}


\def\navod#1{$\,{#1}^\circ$}




\def\eqref#1{$(\ref{#1})$}
\def\bin#1#2{\biggl({#1\atop#2}\biggr)}
\def\dbin#1#2{{\displaystyle\biggl({#1\atop#2}\biggr)}}


\renewcommand{\theproposition}{\thesection.\arabic{proposition}}
\renewcommand{\thetheorem}{\thesection.\arabic{theorem}}
\renewcommand{\thelemma}{\thesection.\arabic{lemma}}
\renewcommand{\thecorollary}{\thesection.\arabic{corollary}}
\renewcommand{\thetable}{\thesection.\arabic{table}}
\renewcommand{\thefigure}{\thesection.\arabic{figure}}
\renewcommand{\thedefinition}{\thesection.\arabic{definition}}
\renewcommand{\theremark}{\thesection.\arabic{remark}}
\renewcommand{\theequation}{\thesection.\arabic{equation}}



\renewcommand{\currentvolume}{44}          %! to be inserted by the editors
\renewcommand{\currentissue}{1}            %! to be inserted by the editors
\renewcommand{\currentyear}{2015}          %! to be inserted by the editors


\pagespan{0} {0}                         %! to be inserted by the editors

\begin{document}

\setcounter{firstpg}{1}                   %! to be inserted by the editors


            %!               VARIANT 1: single author  (please delete if not the case)
\title[Weighted quadrature formulas  for semi-infinite range integrals]{Weighted quadrature formulas  \\for semi-infinite range integrals$^\ast$}   %as a rule, replace <text> with your own text

\author[Gradimir V. Milovanovi\'{c}]{Gradimir V. Milovanovi\'{c}$^\dag$}

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

\thanks{$^\dag$Serbian Academy of Sciences and Arts, Beograd, Serbia  \& State University of Novi Pazar, Serbia, e-mail: {\tt gvm@mi.sanu.ac.rs}.}
            %                End variant 1



\date{<to be inserted by the editor>}

\begin{abstract}
Weighted quadrature formulas on the half line
$(a,+\infty)$,  $a>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<1$, $m\in \NN_0$)
are considered and several numerical examples are  included.
\end{abstract}

\subjclass{65D30, 65D32}

\keywords{Gaussian quadrature rules, nodes, Christoffel numbers, non-exponentially decreasing integrands.}

\dedicatory{Dedicated to prof. Ion P\u{a}v\u{a}loiu on the occasion of his 75th anniversary}

\maketitle

\section{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){\D}x=\sum_{k=1}^n A_k f(x_k)+R_n(f),
\end{equation}
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 {\E}^{-x}$, $\alpha>-1$,  is the well known {\it generalized Gauss-Laguerre} quadrature rule (cf. \cite[p.~325]{gm-gvm-2008}), which is exact for all algebraic polynomials of degree at most $2n-1$, i.e., when $f\in \PP_{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
 \cite{Usp}. Otherwise, the corresponding problems for quadrature
 rules on finite intervals was studied much earlier by \cite{Posse}, Markov \cite{Markov1885}, Stieltjes \cite{Stieltjes1884}, etc. On some new results in this directions see books \cite{GauBook04} and \cite{gm-gvm-2008}, including the so-called truncated quadrature rules obtained by ignoring the last part of its nodes (see Mastroianni and Monegato \cite{gm-gm-2003}).

 Very recently Gautschi \cite{Gautschi2} has constructed a  special logarithmically weighted quadrature formula on $(0,+\infty)$, when $x\mapsto (x-1-\log x)\E^{-x}$. Also, Xu and Milovanovi\'c
 \cite{Xu_GVM} have developed  generalized Gaussian quadrature rules of the form \eqref{GQCBS}, with $x\mapsto w(x)=\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\PP_{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 \cite{ FILOMAT15} and \cite{JCAM15}.



On the other side, a large number of integrals of the form $\int_a^{+\infty}F(x){\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 \cite{Evans2005}). As a starting simple example, Evans \cite{Evans2005} 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 \eqref{GQCBS} of Gaussian type on $(a,+\infty)$, $a>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<1$, $m\in \NN_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.



\section{Generalized Weighted Gaussian Rules}\label{GWGR}
\setcounter{theorem}{0} \setcounter{equation}{0}

Suppose $a>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}\tfrac{w(x)}{x^2}{\D}x<+\infty.	
\end{equation}
Following Evans \cite{Evans2005}, we develop a general approach
for constructing generalized Gaussian quadrature formulas of the
form \eqref{GQCBS}. In the cases of integrals on $(\alpha, +\infty)$, when $\alpha< a$, we simply take
\[\int_\alpha^{+\infty}w(x)f(x){\D}x=\int_\alpha^{a}w(x)f(x){\D}x+
\int_a^{+\infty}w(x)f(x){\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 \eqref{GQCBS}, which is exact for all functions of the form
\[x\mapsto \tfrac1{x^2}P_m\Bigl(\tfrac1{x}\Bigr),\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) \tfrac1{x^2}P_m\Bigl(\tfrac1{x}\Bigr){\D}x=\sum_{k=1}^n  \tfrac{A_k}{x_k^2}P_m\Bigl(\tfrac1{x_k}\Bigr),
\quad m=0,1,\ldots,2n-1.
\end{equation}

\begin{unremark} 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$.\hfil\qed
\end{unremark}	


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\Bigl(\tfrac1{t}\Bigr)g(t){\D}t=
\sum_{k=1}^n B_k g(\tau_k)+R_n^G(g),	
\end{equation}
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 \eqref{uslovi}, such quadrature formulas exist  uniquely, because the
all moments $\mu_k=\int_0^{1/a} w(1/t)t^k{\D}t$, $k\ge 0$, exist
and $\mu_0>0$.

It is known  that the nodes $\tau_k$ in \eqref{IntegFin} are eigenvalues of the following symmetric tridiagonal Jacobi matrix (cf. \cite[pp.~325--328]{gm-gvm-2008})
\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}
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\NN_0}$ orthogonal with respect to the inner product
\begin{equation}\label{inner}
(p,q)=\int_0^{1/a}w\Bigl(\tfrac1{t}\Bigr)p(t)q(t){\D}t.	
\end{equation}
In fact, $\pi_n(t)=(t-\tau_1)\cdots(t-\tau_n)$.


The weight coefficients $B_k$ in (\ref{IntegFin}) 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){\D}t$.

The most popular method for solving this eigenvalue problem is the Golub-Welsch procedure, obtained by a simplification of QR algorithm \cite{GoWe}. This procedure is implemented in several packages
    including the most known {\tt ORTPOL} given by Gautschi \cite{gautschi3}.

 As we can see from \eqref{Jac-matr}, for constructing Gauss--Christoffel quadratures
 \eqref{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 \eqref{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 {\it strongly non-classical polynomials}, these recursion coefficients must be
constructed numerically (cf. \cite{GAU82},  \cite{GauBook04}, \cite[pp.~159--166]{gm-gvm-2008}).
However, recent progress  in  symbolic computation   and  variable-precision arithmetic  today makes it possible to generate  the recursive coefficients in \eqref{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 {\sc Mathematica} package  {\tt OrthogonalPolynomials} (see \cite{Aca} and \cite{MilCvet-MB12}), is downloadable from the web site {\tt  http://www.mi.sanu.ac.rs/\~{}gvm/}. Also, there is Gautschi's software in  {\sc Matlab} (packages {\tt OPQ} and {\tt SOPQ}).
\smallskip

Now, we can give our main result:

\begin{theorem} Let $x\mapsto w(x)$ be a weight function on $(a,+\infty)$, $a>0$,  such that the condition \eqref{uslovi} holds. Assume also that $\tau_k$ and $B_k$, $k=1,\ldots,n$,  are nodes and Christoffel numbers of the Gaussian quadrature formula \eqref{IntegFin},  respectively. Then there exists the generalized Gaussian quadrature formula
\begin{equation}\label{novaQ}
\int_a^{+\infty}w(x) f(x){\D}x=\sum_{k=1}^n A_k f(x_k)+R_n(f),	
\end{equation}
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}
which is exact for all functions of the form
$f(x)=x^{-2}P(x^{-1})$, where $P\in\PP_{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})$.
\end{theorem}

\begin{proof} We start with the system of $2n$ nonlinear equations \eqref{SystemEq}, whose solution determines the parameters of
the quadrature formula \eqref{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 \eqref{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) \tfrac1{x^2}\pi_m\Bigl(\tfrac1{x}\Bigr){\D}x&=&
\int_0^{1/a}w\Bigl(\tfrac1{t}\Bigr)\pi_m(t){\D}t\\[1mm]
&=&(\pi_0,\pi_m)\\[1mm]
&=&\mu_0\delta_{0,m}, 	
\end{eqnarray*}
where the inner product is defined by \eqref{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 \tfrac{1}{x_k^2}\pi_m\Bigl(\tfrac1{x_k}\Bigr)=
\mu_0\delta_{0,m},\quad m=0,1,\ldots,2n-1,	
\end{equation}
but, by an application of the Gaussian rule \eqref{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}
because $R_n^G(g)=0$ for each $g\in\PP_{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
\eqref{IntegFin}.

Since the systems of equations \eqref{sys1} and \eqref{sys2} are equivalent, the statement of this theorem follows directly.	
\end{proof}
	


\section{ Special cases and numerical examples}\label{SpecialC}
\setcounter{theorem}{0} \setcounter{equation}{0}
\setcounter{table}{0} \setcounter{figure}{0}


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<1$ and $m\in \NN_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<1$.

The condition \eqref{uslovi} is satisfied, because
\[\int_a^{+\infty} \tfrac{w(x)}{x^2}{\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){\D}x=a\int_1^{+\infty}f(ax){\D}x,\]
we see that for this important case the following statement holds.

\begin{proposition}\label{propKOEF} Let
\begin{equation}\label{const}
\int_a^{+\infty}f(x){\D}x=\sum_{k=1}^n A_k(a)f(x_k(a))+R_n(f;a),\quad a>0,	
\end{equation}
be a generalized Gaussian quadrature \eqref{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.\]
\end{proposition}

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

Recursive coefficients in \eqref{TTRR}, in this case for translated monic Legendre polynomials, are
\[\alpha_k=\tfrac12,\ \ k\ge0,\quad  \beta_0=1,\ \ \beta_k=\tfrac{k^2}{4(4k^2-1)},\ \ k\ge1.\]
Otherwise, it can be obtained using our {\sc Mathematica} Package  {\tt Orthogonal\-Polynomials} in symbolic form (see \cite{Aca} and \cite{MilCvet-MB12}). 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. \cite[160--162]{gm-gvm-2008}:
{\tt
 \begin{verbatim}
  << orthogonalPolynomials`
  mom=Table[1/(k+1), {k,0,79}];
  {al,be} = aChebyshevAlgorithm[mom, Algorithm -> Symbolic]
\end{verbatim}}
These recursive coefficients enable us to construct quadrature formulas
\eqref{IntegFin} for any number of nodes up to $40$.

However, in this Legendre case (translated to $(0,1)$) we can directly use
{\tt aGaussianNodesWeights} routine to construct nodes and weights in the Gaussian quadrature formula \eqref{IntegFin}, as well as ones in
 the quadrature formula  \eqref{novaQ}:
{\tt
 \begin{verbatim}
  << 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]];]
\end{verbatim}}
Thus, in this way for $a=1$, we obtain quadrature parameters
$x_k$ and $A_k$ for each $n=2(2)40$.

\begin{example} In order to show the efficiency of our quadrature rule \eqref{novaQ} we apply it to the integral
\[J(a;c)=\int_a^{+\infty}\tfrac{1}{(x-2)^2+c^2}{\D}x
=\tfrac{1}{2 c}\left[\pi -2 \arctan\left(\tfrac{a-2}{c}\right)\right],\]
for different values of $a>0$ and $c>0$. In Figure \ref{figures12} we present graphics of the function $x\mapsto f(x;c)=1/((x-2)^2+c^2)$ for
$c=1/8,1/4,1/2$, 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 \eqref{const}, we apply it to $J(a;1)$ for $a=1/2$, $1$, $2$, $3$, $4$, and $8$, when $n=2(2)40$.
\begin{figure}[h]
\begin{center}
$$ \includegraphics[width=0.50\textwidth]{GraffunBOJA.pdf}%
\includegraphics[width=0.50\textwidth]{TVRintBOJA.pdf}$$
\end{center}
\caption{The function $x\mapsto f(x;c)$ (left)  and  the integral $a\mapsto J(a;c)$ (right) for $c=1/8$ (blue line), $c=1/4$ (black line), $c=1/2$ (brown line), and $c=1$ (red line).}\label{figures12}%
\end{figure}

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
\[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~\ref{figure3} in a log-scale.
\begin{figure}[h]
\begin{center}
$$ \includegraphics[width=0.80\textwidth]{err_graf1.pdf}$$
\end{center}
\caption{Relative errors $err_n(f(\cdot;1);a)$ in quadrature sums $Q_n(f(\cdot;1);a)$.}\label{figure3}%
\end{figure}
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 {\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~\ref{TableEx1}, where  quadrature approximations and corresponding relative errors are presented for
$(a,c)=(1,1/4)$, $(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 \eqref{const} is very fast.
\begin{table}[h]
 \centering %\vskip -4pt {%\footnotesize
{\small \begin{tabular}
 [c]{|r|r|l||l|l||c|}\hline
 $n$ & \multicolumn{2}{|c||}{$(a,c)=(1,1/4)$} & \multicolumn{2}{|c||}{$(a,c)=(21/10,10^{-6})$} &    $(a,c)=(4,10^{-6})$\\
 \hline
 $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)$  \\
\hline
 \end{tabular}
 }%footnotesize
 \caption{Quadrature sums $Q_n(f(\cdot;c);a)$  and their relative errors $err_n(f(\cdot;c);a)$ for integrals $J(a;c)$.}\label{TableEx1}
 \end{table}
\end{example}

In the sequel we consider quadrature rules with respect to the weight function
$x\mapsto w(x)=x^\beta\log x$, $0\le \beta<1$.
 Here we suppose that $a\ge 1$. The condition
\eqref{uslovi} is satisfied, because
\begin{equation}\label{U1}
0<\int_a^{+\infty}\tfrac{w(x)}{x^2}{\D}x= \tfrac{a^{\beta-1}}{(1-\beta)^2}\left[1+(1-\beta)\log a\right].	
\end{equation}
In this case, the moments
\[\mu_k=\int_0^{1/a}w(1/t)t^k{\D}t=\int_0^{1/a}t^{k-\beta}\log\tfrac1t{\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\ge0.
\end{equation}

Taking the first one hundred  moments ({\tt mom}) (e.g. for
$a=1$ and $\beta=1/4$) and using  {\sc Mathematica} Package  {\tt OrthogonalPolynomials}, we can get the first fifty recurrence coefficients $\alpha_k$ and $\beta_k$ (denoted  by {\tt \{al1,be1\}}) in the three-term recurrence relation (\ref{TTRR})  in a symbolic form
{\tt
 \begin{verbatim}
  << 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]
\end{verbatim}}
For example, first four coefficients are
{\small\[\alpha_0=\frac{9}{49},\,\alpha_1=\frac{209897}{452025},\,  \alpha_2=\frac{6582284926939}{13538179995075},\,  \alpha_3=\frac{7618613698603068100869609}{15464687102113919816429449}
 \]}
and
{\small
\[\beta_0=\frac{16}{9},\  \,\beta_1=\frac{11808}{290521},\  \,\beta_2=\frac{213147564896}{3717280400625},\  \,\beta_3=\frac{421267942813254097088}{6997413354065613077481}.\]}

\begin{example} As a test example we consider the function $x\mapsto f(x)=(x+1)^{-2}$ and integral (see \cite{Evans2005})
\begin{align*}
I(f;a)=&\int_a^{+\infty}\tfrac{x^{1/4}\log x}{(x+1)^2}{\D}x\\
=&\tfrac{1}{36 a^\frac74 (a+1)}\left\{-9 (a+1) \Phi\Bigl(-\tfrac{1}{a},2,\tfrac{7}{4}\Bigr)\right.\qquad\quad\\[2mm]
&\qquad \left.+4 a \left[3 (a+1)(\log a+4) \, _2F_1\Bigl(\tfrac{3}{4},1;\tfrac{7}{4};-\tfrac{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} \tfrac{z^k}{(k+a)^s}\quad\mbox{and}\quad
_2F_1(a,b;c;z)=\sum_{k=0}^{+\infty}\tfrac{(a)_k(b)_k}{(c)_k}\tfrac{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={\E}$, i.e.,
\[ I(f;1)=1.35974328097600895\ldots\ \ \mbox{and}\ \
  I(f;\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 \ref{Tab:one}. As we can see only two two decimal digits are true in quadrature sum with $2048$ nodes!

\begin{table}[h]
  \centering %\vskip -4pt {%\footnotesize
 \begin{tabular}
 [c]{|c|c|c|c|c|c|}\hline
 $n=2$ & $n=8$ & $n=32$ & $n=128$ & $n=512$ & $n=2048$\\
 \hline
 $6.72(-1)$ & $3.60(-1)$ & $1.64(-1)$ & $7.00(-2)$ & $2.90(-2)$ & $1.18(-2)$ \\
 \hline
 \end{tabular}
% }%footnotesize
 \caption{Relative errors in Gauss-Laguerre quadrature sums with
 $n=2,8,32,128,512$ and $2048$ nodes.}\label{Tab:one}
 \end{table}

Now, we apply our quadrature formula \eqref{novaQ}, with parameters
given by \eqref{PARAM}, to $I(f;1)$ and $I(f;\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)$,
\[err_n(f;a)=\Bigl|\frac{Q_{n}(f;a)-I(f;a)}{I(f;a)}\Bigr|, \]
are presented in Table~\ref{Tab:two}.
\begin{table}[h]
 \centering %\vskip -4pt {%\footnotesize
 \begin{tabular}
 [c]{|c|c|c|c|c|c|c|}\hline
 $a$ & $n=2$ & $n=4$ & $n=6$ & $n=8$ & $n=10$ & $n=12$\\
 \hline
 $1$ & $2.94(-3)$ & $4.24(-6)\ \,$ & $5.15(-9)\ $ & $5.72(-12)$ & $4.74(-13)$ & $7.07(-13)$ \\
  ${\E}$ & $2.40(-4)$ & $1.64(-8)\ \,$ & $\ 8.91(-13)$ & $8.83(-14)$ & $5.31(-14)$ & $3.80(-14)$ \\
   ${\E}^2$ & $7.18(-6)$ & $1.28(-11)$ & $\ 3.10(-14)$ &  &  &  \\
 \hline
 \end{tabular}
% }%footnotesize
 \caption{Relative errors $err_n(f;a)$ in quadrature sums $Q_n(f;a)$ for
 different number of nodes $n$ and three values of $a$ $(=1,{\E},\ \mbox{and}\ {\E}^2)$.} \label{Tab:two}
 \end{table}

 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={\E}^2$ and $n=2,4$, and $6$.
\end{example}

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<1,\ m=2,3,\ldots\ ,\]
on the interval $(a,+\infty)$, with $a\ge 1$.

The condition
\eqref{uslovi} for $U_m=\int_a^{+\infty}x^{-2}w_m(x){\D}x$
is also satisfied, because
\[U_m=\tfrac1{1-\beta}\left[mU_{m-1}+a^{\beta-1}\log^m a\right],\quad m=2,3,\ldots\,,\]
where $U_1$ is given in \eqref{U1}. The corresponding moments
\[\mu_k^{[m]}=\int_0^{1/a}t^{k-\beta}\log^m\tfrac1{t}{\D}t,\quad k=0,1,\ldots\,,\]
can be expressed recursively in terms of the moments $\mu_k^{[m-1]}$,
\[\mu_k^{[m]}=\tfrac{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 \eqref{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\ge0.\]


The corresponding recursive coefficients, for example for $\beta=0$
and $a=1$, are
{\small\[\alpha_0=\frac{1}{8},\,\alpha_1=\frac{115}{296},\,  \alpha_2=\frac{28200187}{62721512},\,  \alpha_3=\frac{28003451041760695}{59414538084233528}, \ \ldots
 \]}
and
{\small
\[\beta_0=2,\,\beta_1=\frac{37}{1728},\ \,\beta_2=\frac{211897}{4620375},\  \,\beta_3=\frac{945381680572419}{17600932734728000},\ \ldots\
.\]}

\begin{example} We consider the function $x\mapsto 1/(1+x^2)$ and the corresponding weighted integral over $(a,+\infty)$
\[I(f;a)=\int_a^{+\infty}\tfrac{\log^2x}{1+x^2}{\D}x,\]
for two values of $a$ $(=1$ and $={\E})$, for which
\[I(f;1)=1.93789229251873876096726969169\ldots\]
and
\[\ \ I(f;{\E})=\numprint{1.80988687939786942602016447246}\ldots\ .\]

Applying  our quadrature formula \eqref{novaQ}, with parameters
given by \eqref{PARAM}, to $I(f;1)$ and $I(f;\E)$, with  $n=2(2)12$ nodes, we get the corresponding quadrature approximations
$Q_n(f;a)$ with the relative errors $err_n(f;a)$ presented in Table~\ref{Tab:three}.
\begin{table}[h]
 \centering %\vskip -4pt {%\footnotesize
 \begin{tabular}
 [c]{|c|c|c|c|c|c|c|}\hline
 $a$ & $n=2$ & $n=4$ & $n=6$ & $n=8$ & $n=10$ & $n=12$\\
 \hline
 $1$ & $1.66(-4)$ & $1.31(-6)\ \,$ & $1.98(-10)\ $ & $5.73(-12)$ & $2.08(-15)$ & $2.56(-17)$ \\
  ${\E}$ & $5.33(-5)$ & $5.04(-10)\ \,$ & $\ 1.86(-13)$ & $2.05(-17)$ & $1.22(-21)$ & $3.30(-26)$ \\
 \hline
 \end{tabular}
% }%footnotesize
 \caption{Relative errors $err_n(f;a)$ in quadrature sums $Q_n(f;a)$ for  different number of nodes $n$ and two
 values of $a$ $(=1$ and $={\E})$.}\label{Tab:three}
 \end{table}

Here also we can note a faster convergence when $a$  has a  larger value.
\end{example}

\begin{thebibliography}{9}

{\footnotesize

\bibitem{Aca}
{\sc A.S. Cvetkovi\'c} and {\sc G.V. Milovanovi\'c}, {\it  The Mathematica Package ``OrthogonalPolynomials''}, Facta Univ. Ser. Math. Inform.,  {\bf19}, pp. 17--36,  2004.

\bibitem{Evans2005}
{\sc G.A. Evans}, {\it  Some new thoughts on Gauss-Laguerre quadrature}, Int. J. Comput. Math. {\bf 82}, pp. 721--730,   2005.

\bibitem{GAU82}
{\sc W. Gautschi}, {\it On generating orthogonal polynomials}, SIAM J. Sci. Statist. Comput., {\bf3}, pp. 289--317, 1982.

\bibitem{gautschi3}
{\sc W. Gautschi}, {\it Algorithm 726: ORTHPOL -- A package of routines for generating orthogonal polynomials and Gauss-type quadrature rules},  ACM Trans. Math. Software, {\bf 20}, pp. 21--62,  1994.


\bibitem{GauBook04}
{\sc W. Gautschi},{\it Orthogonal Polynomials: Computation and
Approximation},  Clarendon Press, Oxford, 2004.

\bibitem{Gautschi2}
{\sc W. Gautschi},{\it Gauss quadrature routines for two classes of logarithmic weight functions},  Numer. Algor. {\bf 55}, pp. 265--277, 2010.

\bibitem{GoWe}
{\sc G. Golub} and {\sc J.H. Welsch},
 {\it Calculation of Gauss quadrature rules}, Math. Comp.
{\bf 23}, pp. 221--230, 1969.

\bibitem{Markov1885}
{\sc A. Markov}, {\it Sur la m\'ethode de Gauss pour le calcul approch\'e  des int\'egrales},  Math. Ann., {\bf25}, pp.~427--432, 1885.

\bibitem{gm-gm-2003} {\sc G. Mastroianni} and {\sc G. Monegato},
{\it Truncated quadrature rules over $(0, \infty)$ and
Nystr\"om-type methods}, SIAM J. Numer. Anal. {\bf41}, pp. 1870--1892, 2003.

\bibitem{gm-gvm-2008} {\sc G. Mastroianni} and {\sc G.V. Milovanovi\'c}, {\it  Interpolation Processes -- Basic Theory and Applications},  Springer-Verlag, Berlin -- Heidelberg -- New York, 2008.

\bibitem{Stud15}
{\sc G.V. Milovanovi\'c}, {\it  Construction and applications of Gaussian quadratures with nonclassical and exotic weight function}, Stud. Univ. Babe\c s-Bolyai Math., {\bf  60}, pp. 211--233, 2015.

\bibitem{FILOMAT15}
{\sc G.V. Milovanovi\'c}, {\it  Generalized Gaussian quadratures for integrals with logarithmic singularity}, FILOMAT  (to appear).


 \bibitem{MilCvet-MB12}  {\sc G.V. Milovanovi\'c} and  {\sc A.S. Cvetkovi\'c}, {\it  Special classes of orthogonal polynomials and corresponding quadratures of Gaussian type}, Math. Balkanica, {\bf 26}, pp. 169--184, 2012.


\bibitem{JCAM15}
{\sc G.V. Milovanovi\'c, T.S. Igi\'c} and {\sc D. Turni\'c},  {\it Generalized quadrature rules of Gaussian type for numerical evaluation of singular integrals}, J. Comput. Appl. Math. {\bf 278}, pp. 306--25, 2015.

\bibitem{Posse}
{\sc C. Posse}, {\it Sur les quadratures}, Nouv. Ann. Math. (2) {\bf14}, pp. 49--62, 1875.

\bibitem{Stieltjes1884}
{\sc T.J. Stieltjes}, {\it Quelques recherches sur la th\'eorie des quadratures dites m\'ecaniques},  Ann. Sci. \'Ec. Norm. Paris, S\'er. 2, {\bf1}, pp. 409--426, 1884.

 \bibitem{Usp}  {\sc J.V. Uspensky}, {\it On the convergence of quadrature formulas related to an infinite interval}, Trans. Amer. Math. Soc., {\bf 30}, pp. 542--559, 1928.

\bibitem{Xu_GVM}
{\sc Zhenhua Xu} and {\sc G.V. Milovanovi\'c}, {\it
 Efficient method for the computation of oscillatory Bessel  transform and Bessel Hilbert transform} (submitted).

}%\footnotesize

\end{thebibliography}
\end{document}