Gauss Quadrature for Integrals and Sums

1. Introduction

Approximating integrals, especially intractable ones, by finite sums is as old of a problem as calculus itself. Most approximation procedures follow the structure of the original classic discretized approximation of the integral of a real function f(x) defined in $x_{-}\le x\le x_{+}$ as

$${\displaystyle {\int }_{x_{-}}^{x_{+}}f(x)\hspace{0.17em}dx}\cong {\displaystyle \sum _{n=0}^{N-1}{\Delta }_{n}f(x_n)}$$

(0)

for some large enough integer N. The accuracy of the approximation and its convergence have always been at the core of the development of such approximation schemes. One of the most prominent techniques used in these approximations in the celebrated Gauss quadrature integral approximation [1,2,3]. It exploits the orthogonality and recursion relation of orthogonal polynomials. The quadrature has also been extended to the approximation of infinite sums (see, for example, [4] and [5] and references therein). In this work, we combine the two. That is, we propose an extended Gauss quadrature associated with integral measures that consist of discrete as well as continuous components. In Section 2, we summarize the conventional Gauss quadrature (GQ) integral approximation where the integration measure is purely continuous. In Section 3, we summarize GQ for approximating infinite sums by finite ones. Finally, in Section 4, we present our main findings where we combine the two.

2. The Conventional GQ: Purely Continuous Measure

Let $f(x)$ be an integrable real valued function with respect to the positive measure $d\mu (x)=$ $\rho (x)\hspace{0.17em}dx$ in a real vector space, which is spanned by the complete set of orthonormal basis ${\left\{p_n(x)\right\}}_{n=0}^{\infty }$, where $x\in [x_{-},x_{+}]\subset ℝ$. Orthogonality is defined as

$${\displaystyle {\int }_{x_{-}}^{x_{+}}\rho (x)p_n(x)p_m(x)\hspace{0.17em}dx}={\delta }_{n,m}$$

(1)

This suggests that we could take $p_n(x)$ to be an orthogonal polynomial of degree n in x with pure continuous spectrum and where $\rho (x)$ is the associated positive weight function [6,7,8]. The spectral theorem (a.k.a Favard’s theorem) dictates that such polynomials satisfy the following symmetric three-term recursion relation

$$x\hspace{0.17em}{p}_n(x)=a_n{p}_n(x)+b_{n-1}{p}_{n-1}(x)+b_n{p}_{n+1}(x)$$

(2)

for $n=1,2,3,\mathrm{...}$ and where the “recursion coefficients” $\left\{a_n,b_n\right\}$ are real constants such that $b_n\ne 0$ for all $n\ge 0$. This recursion gives all the polynomials of any degree starting with the two seed values $p_0(x)=1$ and $p_1(x)=(x-a_0)/b_0$ (i.e., $b_{-1}\equiv 0$). Associated with this space is an infinite dimensional real tridiagonal symmetric matrix (Jacobi matrix) whose elements are

$$J_{n,m}=a_n{\delta }_{n,m}+b_{n-1}{\delta }_{n,m+1}+b_n{\delta }_{n,m-1}={\displaystyle {\int }_{x_{-}}^{x_{+}}\rho (x)\hspace{0.17em}{p}_n(x)\hspace{0.17em}x\hspace{0.17em}{p}_m(x)\hspace{0.17em}dx},$$

(3a)

$$J=\left(\begin{array}{ccccc}{a}_0& b_0& & & \\ b_0& a_1& b_1& & \\ & b_1& a_2& b_2& \\ & & \times & \times & \times \\ & & & \times & \times \end{array}\right)$$

(3b)

For numerical computations, however, the space is truncated to a finite N-dimensional subspace spanned by ${\left\{p_n(x)\right\}}_{n=0}^{N-1}$. The tridiagonal matrix (3) becomes a finite N×N matrix J. The real N distinct eigenvalues of J, which we designate as the set ${\left\{{\epsilon }_n\right\}}_{n=0}^{N-1}$, are the zeros of the polynomial $p_N(x)$ [i.e. $p_N({\epsilon }_n)=0$]. Let ${\left\{{\Lambda }_{m,n}\right\}}_{m=0}^{N-1}$ be the normalized eigenvector of J associated with the eigenvalue ${\epsilon }_n$. In this setting, Gauss quadrature integral approximation states that [1,2,3]

$${\displaystyle {\int }_{x_{-}}^{x_{+}}\rho (x)f(x)\hspace{0.17em}dx}\cong {\displaystyle \sum _{n=0}^{N-1}{\omega }_{n}f({\epsilon }_n)},$$

(4)

where the “numerical weights” could be evaluated as ${\omega }_n={\Lambda }_{0,n}^2$. Due to the lower numerical cost in computing matrix eigenvalues instead of eigenvectors, we can also write these numerical weights in terms of ${\left\{{\epsilon }_n\right\}}_{n=0}^{N-1}$ and another set of eigenvalues ${\left\{{\widehat{\epsilon }}_m\right\}}_{m=0}^{N-2}$ as

$${\omega }_n=\frac{{\displaystyle \prod _{m=0}^{N-2}\left({\epsilon }_n-{\widehat{\epsilon }}_m\right)}}{{\displaystyle \prod _{\begin{array}{l}k=0\\ k\ne n\end{array}}^{N-1}\left({\epsilon }_n-{\epsilon }_k\right)}},$$

(5)

where ${\left\{{\widehat{\epsilon }}_m\right\}}_{m=0}^{N-2}$ is the set of eigenvalues of a submatrix of J obtained by deleting the first (zeroth) row and first column. If sorted, these eigenvalues interlace as ${\epsilon }_0<{\widehat{\epsilon }}_0<{\epsilon }_1<{\widehat{\epsilon }}_1<{\epsilon }_2\mathrm{....}{\widehat{\epsilon }}_{N-2}<{\epsilon }_{N-1}$. The integral approximation (4) becomes exact if $f(x)$ is a polynomial in x of a degree less than or equal to $2N-1$. If, instead of the integral (4), we have ${\displaystyle {\int }_{x_{-}}^{x_{+}}f(x)dx}$ then this could be approximated as follows

$${\displaystyle {\int }_{x_{-}}^{x_{+}}f(x)\hspace{0.17em}dx}={\displaystyle {\int }_{x_{-}}^{x_{+}}\rho (x)\frac{f(x)}{\rho (x)}dx}\cong {\displaystyle \sum _{n=0}^{N-1}{\omega }_n\frac{f({\epsilon }_n)}{\rho ({\epsilon }_n)}}:={\displaystyle \sum _{n=0}^{N-1}{\tilde{\omega }}_{n}f({\epsilon }_n)},$$

(6)

where $\left\{{\tilde{\omega }}_n\right\}$ are named the “derivative weights” and ${\tilde{\omega }}_n={\omega }_n/\rho ({\epsilon }_n)$.

Very often, one encounters integrals that represent matrix elements of functions in the basis set $\left\{{\varphi }_n(x)\right\}=\left\{\sqrt{\rho (x)}\hspace{0.17em}{p}_n(x)\right\}$ (for example, the matrix elements of a potential function in quantum mechanics) of the form

$$f_{n,m}={\displaystyle {\int }_{x_{-}}^{x_{+}}{\varphi }_n(x)f(x){\varphi }_m(x)\hspace{0.17em}dx}={\displaystyle {\int }_{x_{-}}^{x_{+}}\rho (x)p_n(x)f(x)p_m(x)\hspace{0.17em}dx}.$$

(7)

Using Eq. (4), this integral is approximated as follows

$$f_{n,m}\cong {\displaystyle \sum _{k=0}^{N-1}{\omega }_k{p}_n({\epsilon }_k)f({\epsilon }_k)p_m({\epsilon }_k)},$$

(8)

Now, it is well established that $p_n({\epsilon }_k)={\Lambda }_{n,k}/{\Lambda }_{0,k}$ for $n,k=0,1,2,\mathrm{..},N-1$. Therefore, with ${\omega }_k={\Lambda }_{0,k}^2$ this could be rewritten in matrix form as

$$f_{n,m}\cong {\displaystyle \sum _{k=0}^{N-1}{\Lambda }_{n,k}f({\epsilon }_k){\Lambda }_{m,k}}={\left(\Lambda F{\Lambda }^{\mathrm{T}}\right)}_{n,m},$$

(9)

where F is a diagonal matrix with elements $F_{k,j}=f({\epsilon }_k)\hspace{0.17em}{\delta }_{k,j}$. Therefore, to obtain an approximate evaluation of integrals using Gauss quadrature associated with orthogonal polynomials satisfying (1) and (2), one needs only the tridiagonal symmetric matrix J, which is constructed using the recursion coefficients $\left\{a_n,b_n\right\}$, and possibly the weight function $\rho (x)$ for integrals of the type (6).

3. GQ for a Purely Discrete Measure

If instead of the integral ${\displaystyle {\int }_{x_{-}}^{x_{+}}\rho (x)f(x)dx}$, we were tasked with approximating the infinite weighted sum ${\displaystyle {\sum }_{k=0}^{\infty }{\xi }_{k}f(x_k)}$, for a given discrete set $\left\{{\xi }_k,x_k\right\}$, by a finite sum similar to that of Eq. (4). Then, we could approach the solution to the problem as follows. We can convert the integral ${\displaystyle {\int }_{x_{-}}^{x_{+}}\rho (x)f(x)dx}$ into the sum ${\displaystyle {\sum }_{k=0}^{\infty }{\xi }_{k}f(x_k)}$ by discretizing the weight function and writing $\rho (x)=\chi (x){\displaystyle {\sum }_{k=0}^{\infty }\delta (x-x_k)}$ with ${\xi }_k=\chi (x_k)$. Performing this same weight function discretization in the orthogonality (1) turns it into

$${\displaystyle \sum _{k=0}^{\infty }{\xi }_k{p}_n(x_k)p_m(x_k)}={\delta }_{n,m}.$$

(10)

There are many discrete polynomials that satisfy this kind of orthogonality. The Appendix lists relevant properties of few of these discrete hypergeometric orthogonal polynomials that belong to the Askey scheme [9]. Some of these polynomials constitute an infinite discrete sequence whereas others form finite sequences. Examples of the former are the Charlier and Meixner polynomials. Whereas examples of the latter include the Racah, dual Hahn and Krawtchouk polynomials [9]. All these polynomials do satisfy symmetric three-term recursion relations similar to (2) that read

$$x_k{p}_n(x_k)=a_n{p}_n(x_k)+b_{n-1}{p}_{n-1}(x_k)+b_n{p}_{n+1}(x_k).$$

(11)

Consequently, we can write the following approximation

$${\displaystyle \sum _{k=0}^M{\xi }_{k}f(x_k)}\cong {\displaystyle \sum _{n=0}^{N-1}{\omega }_{n}f({\epsilon }_n)},$$

(12)

where M can be a finite or an infinite integer. However, typically $M>>N$ otherwise the approximation is not called for. Moreover, the set ${\left\{{\epsilon }_n,{\omega }_k\right\}}_{n=0}^{N-1}$ in (12) are calculated using the tridiagonal symmetric matrix J obtained from (11) exactly as in Section 2 and shown by Eq. (3b). Additionally, the approximation (12) becomes exact if $f(x)$ is a polynomial in x of a degree less than or equal to $2N-1$. Here too, if instead of the sum (12) we have ${\displaystyle {\sum }_{k=0}^{M}f(x_k)}$, then we can write

$${\displaystyle \sum _{k=0}^{M}f(x_k)}\cong {\displaystyle \sum _{n=0}^{N-1}{\tilde{\omega }}_{n}f({\epsilon }_n)},$$

(13)

where the derivative weights are calculated as ${\tilde{\omega }}_n={\omega }_n/\chi ({\epsilon }_n)$. As an example, we consider the function $f(x)=r^x/\Gamma (x+1)$ with $(x,r)>0$ and obtain the approximation of its infinite sum (13) using GQ associated with the Charlier and Meixner polynomials. Since $x_k=k=0,1,2,\mathrm{...}$ for both polynomials, then the infinite sum has the exact value $e^r$. Table 1 demonstrates the validity of formula (13) and illustrates the relative accuracy and convergence associated with the two polynomials for this particular function. In another example, we consider the finite sum (13) for the function $f(x)=(x+1)\hspace{0.17em}{r}^{x+1}/\Gamma (x+r+2)$ with $M=100$ and employ GQ associated with the finite Krawtchouk polynomial. Since $x_k=k$, then this finite sum has the following exact closed form (see formula 3, Section 1.3.9 in Ref. [6])

$${\displaystyle \sum _{k=0}^M\frac{(k+1)\hspace{0.17em}{r}^{k+1}}{\Gamma (k+r+2)}}=\frac{1}{\Gamma (r)}-\frac{r^{M+2}}{\Gamma (M+r+2)}.$$

(14)

Table 2 demonstrates the validity of formula (13) and illustrates its accuracy as well as its convergence as we vary N and the polynomial parameter γ.

4. GQ for a Measure with a Mix of Continuous and Discrete Spectra

If the integral measure $d\mu (x)=\rho (x)\hspace{0.17em}dx$ contains a mix of discrete as well as continuous spectra, then the weight function could be written as follows

$$\rho (x)=\sigma (x)+\chi (x){\displaystyle \sum _{k=0}^M\delta (x-x_k)},$$

(15)

where $\sigma (x)$ is the continuous component of the weight function and $\chi (x_k):={\xi }_k$ is the discrete component. The discrete spectrum (a.k.a mass points) is the set ${\left\{x_k\right\}}_{k=0}^M$ of size $M+1$, which could be finite or infinite. The orthogonality (1) becomes [6,7,8]

$${\displaystyle {\int }_{x_{-}}^{x_{+}}\sigma (x)p_n(x)p_m(x)\hspace{0.17em}dx}+{\displaystyle \sum _{k=0}^M{\xi }_k{p}_n(x_k)p_m(x_k)}={\delta }_{n,m}.$$

(16)

These polynomials do also satisfy a symmetric three-term recursion relation similar to (2) where the associated tridiagonal symmetric matrix J is endowed with both continuous as well as discrete eigenvalues. Appendix A lists the relevant properties of two of these hypergeometric orthogonal polynomials that belong to the Askey scheme [9].

If the function $f(x)$ is defined over $x\in [x_{-},x_{+}]\cup {\left\{x_k\right\}}_{k=0}^M$ then Gauss quadrature integral approximation (4) could be extended to read

$${\displaystyle {\int }_{x_{-}}^{x_{+}}\sigma (x)f(x)\hspace{0.17em}dx}+{\displaystyle \sum _{k=0}^M{\xi }_{k}f(x_k)}\cong {\displaystyle \sum _{n=0}^{N-1}{\omega }_{n}f({\epsilon }_n)},$$

(17)

where the set ${\left\{{\epsilon }_n,{\omega }_k\right\}}_{n=0}^{N-1}$ in (17) are calculated using the tridiagonal symmetric matrix J exactly as in Section 2. Here too, the integral approximation (17) becomes exact if $f(x)$ is a polynomial in x of a degree less than or equal to $2N-1$. If the sum in (17) is finite such that M is of the same order as N (i.e., $M\underset{\sim }{<}N$), then one could think of (17) as an integral approximation by two finite sums as follows

$${\displaystyle {\int }_{x_{-}}^{x_{+}}\sigma (x)f(x)\hspace{0.17em}dx}\cong {\displaystyle \sum _{n=0}^{N-1}{\omega }_{n}f({\epsilon }_n)}-{\displaystyle \sum _{k=0}^M{\xi }_{k}f(x_k)}$$

(18)

As an example, we consider the function $f(x)=x^r{e}^{-x/2}$ with $r>0$ and obtain the approximation of its integral using GQ associated with the continuous dual Hahn polynomial $S_n^{\mu }(x^2;\alpha ,\alpha )$ with a mixed spectrum (i.e., $\mu <0$). Since the polynomial argument is not x but $x^2$and $x_k^2=-{(k+\mu )}^2$, then the GQ integral formula (18) must be revised to read

$${\displaystyle {\int }_0^{\infty }\sigma (x)f(x^2)\hspace{0.17em}dx}\cong {\displaystyle \sum _{n=0}^{N-1}{\omega }_{n}f({\epsilon }_n)}-{\displaystyle \sum _{k=0}^M{\xi }_{k}f(x_k^2)}$$

(19)

Table 3 demonstrates the validity this formula and illustrates its convergence for a given μ and several values of α. We compare GQ results to the numerical integration of ${\displaystyle {\int }_0^{\infty }\sigma (x)f(x^2)\hspace{0.17em}dx}$ obtained by the built-in Mathcad integration scheme with an accuracy set by the precision variable (TOL) set at ${10}^{-14}$ and tested by showing the relative error in the calculation of ${\displaystyle {\int }_0^{\infty }\sigma (x)\hspace{0.17em}dx}+{\displaystyle \sum _{k=0}^M{\xi }_k}$ whose exact value is one. It is worth observing that the accuracy of the approximation diminishes if the quadrature degree N is increased beyond a certain critical maximum value that depends on α for a fixed μ. For comparison purposes, Table 4 gives the integral approximation of the same function using GQ associated with the continuous dual Hahn polynomial $S_n^{\mu }(x^2;\alpha ,\alpha )$ but with pure continuous spectrum (i.e., $\mu \mapsto \left|\mu \right|>0$ and $\alpha -\mu \mapsto \alpha$). That is, we use the GQ formula (4). Evidently, the results in Table 4 are more accurate than those in Table 3 for this particular function. However, one should be aware that the comparison is imprecise because of the parameter change $\alpha -\mu \mapsto \alpha$.

We plan to follow this work by another in which we attempt at approximating the following integral-sum combination that replaces (17)

$${\displaystyle {\int }_{x_{-}}^{x_{+}}f(x)\hspace{0.17em}dx}+{\displaystyle \sum _{k=0}^{M}f(x_k)}\cong {\displaystyle \sum _{n=0}^{N-1}{\tilde{\omega }}_{n}f({\epsilon }_n)},$$

(20)

where we determine the derivative weights ${\tilde{\omega }}_n$ in terms of ${\omega }_n$, $\sigma ({\epsilon }_n)$ and $\chi ({\epsilon }_n)$.