Gauss Quadrature for Integrals and Sums

Preprint

Article

Gauss Quadrature for Integrals and Sums

Altmetrics

Downloads

304

Views

Comments

A peer-reviewed article of this preprint also exists.

Abdulaziz D. Alhaidari^*

This version is not peer-reviewed

Submitted:

16 February 2023

Posted:

17 February 2023

You are already at the latest version

Alerts

Abstract

Gauss quadrature integral approximation is extended to include integrals with a measure consisting of a continuous as well as a discrete component. That is, we give an approximation for the integral of a function plus its sum over a discrete weighted set.

Keywords:

Subject: Computer Science and Mathematics - Computational Mathematics

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_{-} \leq x \leq x_{+}

\int_{x_{-}}^{x_{+}} f (x) d x ≅ \sum_{n = 0}^{N - 1} Δ_{n} f (x_{n}),

(1)

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 μ (x) =

ρ (x) d x

in a real vector space, which is spanned by the complete set of orthonormal basis

{\{p_{n} (x)\}}_{n = 0}^{\infty}

, where

x \in [x_{-}, x_{+}] \subset ℝ

. Orthogonality is defined as

\int_{x_{-}}^{x_{+}} ρ (x) p_{n} (x) p_{m} (x) d x = δ_{n, m} .

(2)

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

ρ (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 p_{n} (x) = a_{n} p_{n} (x) + b_{n - 1} p_{n - 1} (x) + b_{n} p_{n + 1} (x),

(3)

for

n = 1, 2, 3, ...

and where the “recursion coefficients”

\{a_{n}, b_{n}\}

are real constants such that

b_{n} \neq 0

for all

n \geq 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 (a Jacobi matrix) whose elements are

J_{n, m} = a_{n} δ_{n, m} + b_{n - 1} δ_{n, m + 1} + b_{n} δ_{n, m - 1} = \int_{x_{-}}^{x_{+}} ρ (x) p_{n} (x) x p_{m} (x) d x,

(4a)

J = (\begin{matrix} a_{0} & b_{0} \\ b_{0} & a_{1} & b_{1} \\ b_{1} & a_{2} & b_{2} \\ \times & \times & \times \\ \times & \times \end{matrix}) .

(4b)

One can easily show that by iterating (4a) for x,

x^{2}

x^{3}

,... and using the three-term recursion relation (3) repeatedly, we obtain

\int_{x_{-}}^{x_{+}} ρ (x) p_{n} (x) x^{k} p_{m} (x) d x = {(J^{k})}_{n, m} .

(4c)

For numerical computations, however, the space is truncated to a finite N-dimensional subspace spanned by

{\{p_{n} (x)\}}_{n = 0}^{N - 1}

. The tridiagonal matrix (4) becomes a finite N×N matrix J. The real N distinct eigenvalues of J, which we designate as the set

{\{ε_{n}\}}_{n = 0}^{N - 1}

, are the zeros of the polynomial

p_{N} (x)

[i.e.

p_{N} (ε_{n}) = 0

]. Let

{\{Λ_{m, n}\}}_{m = 0}^{N - 1}

be the normalized eigenvector of J associated with the eigenvalue

ε_{n}

. In this setting, Gauss quadrature integral approximation states that [1,2,3]

\int_{x_{-}}^{x_{+}} ρ (x) f (x) d x ≅ \sum_{n = 0}^{N - 1} ω_{n} f (ε_{n}),

(5)

where the “numerical weights” could be evaluated as

ω_{n} = Λ_{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

{\{ε_{n}\}}_{n = 0}^{N - 1}

and another set of eigenvalues

{\{{\hat{ε}}_{m}\}}_{m = 0}^{N - 2}

ω_{n} = \frac{\prod_{m = 0}^{N - 2} (ε_{n} - {\hat{ε}}_{m})}{\prod_{\begin{array}{l} k = 0 \\ k \neq n \end{array}}^{N - 1} (ε_{n} - ε_{k})},

(6)

where

{\{{\hat{ε}}_{m}\}}_{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

ε_{0} < {\hat{ε}}_{0} < ε_{1} < {\hat{ε}}_{1} < ε_{2} .... {\hat{ε}}_{N - 2} < ε_{N - 1}

. The integral approximation (5) becomes exact if

f (x)

is a polynomial in x of a degree less than or equal to

2 N - 1

. If, instead of the integral (5), we have

\int_{x_{-}}^{x_{+}} f (x) d x

then this could be approximated as follows

\int_{x_{-}}^{x_{+}} f (x) d x = \int_{x_{-}}^{x_{+}} ρ (x) \frac{f (x)}{ρ (x)} d x ≅ \sum_{n = 0}^{N - 1} ω_{n} \frac{f (ε_{n})}{ρ (ε_{n})} : = \sum_{n = 0}^{N - 1} {\tilde{ω}}_{n} f (ε_{n}),

(7)

where

\{{\tilde{ω}}_{n}\}

are named the “derivative weights” and

{\tilde{ω}}_{n} = ω_{n} / ρ (ε_{n})

Very often, one encounters integrals that represent matrix elements of functions in the basis set

\{ϕ_{n} (x)\} = \{\sqrt{ρ (x)} p_{n} (x)\}

(for example, the matrix elements of a potential function in quantum mechanics) of the form

f_{n, m} = \int_{x_{-}}^{x_{+}} ϕ_{n} (x) f (x) ϕ_{m} (x) d x = \int_{x_{-}}^{x_{+}} ρ (x) p_{n} (x) f (x) p_{m} (x) d x .

(8)

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

f_{n, m} ≅ \sum_{k = 0}^{N - 1} ω_{k} p_{n} (ε_{k}) f (ε_{k}) p_{m} (ε_{k}) .

(9)

Now, it is well established that

p_{n} (ε_{k}) = Λ_{n, k} / Λ_{0, k}

for

n, k = 0, 1, 2, .., N - 1

. Therefore, with

ω_{k} = Λ_{0, k}^{2}

this could be rewritten in matrix form as

f_{n, m} ≅ \sum_{k = 0}^{N - 1} Λ_{n, k} f (ε_{k}) Λ_{m, k} = {(Λ F Λ^{T})}_{n, m},

(10)

where F is a diagonal matrix with elements

F_{k, j} = f (ε_{k}) δ_{k, j}

. Therefore, to obtain an approximate evaluation of integrals using Gauss quadrature associated with orthogonal polynomials satisfying (2) and (3), one needs only the tridiagonal symmetric matrix J, which is constructed using the recursion coefficients

\{a_{n}, b_{n}\}

, and possibly the weight function

ρ (x)

for integrals of the type (7).

3. GQ for a purely discrete measure

If instead of the integral

\int_{x_{-}}^{x_{+}} ρ (x) f (x) d x

, we were tasked with approximating the infinite weighted sum

\sum_{k = 0}^{\infty} ξ_{k} f (x_{k})

, for a given discrete set

\{ξ_{k}, x_{k}\}

, by a finite sum similar to that of Eq. (5). Then, we could approach the solution to the problem as follows. We can convert the integral

\int_{x_{-}}^{x_{+}} ρ (x) f (x) d x

into the sum

\sum_{k = 0}^{\infty} ξ_{k} f (x_{k})

by discretizing the weight function and writing

ρ (x) = χ (x) \sum_{k = 0}^{\infty} δ (x - x_{k})

with

ξ_{k} = χ (x_{k})

. Performing this same weight function discretization in the orthogonality (2) turns it into

\sum_{k = 0}^{\infty} ξ_{k} p_{n} (x_{k}) p_{m} (x_{k}) = δ_{n, m} .

(11)

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 (3) 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}) .

(12)

Consequently, we can write the following approximation

\sum_{k = 0}^{M} ξ_{k} f (x_{k}) ≅ \sum_{n = 0}^{N - 1} ω_{n} f (ε_{n}),

(13)

where M can be a finite or an infinite integer. However, typically

M > > N

otherwise the approximation is not called for. Moreover, the set

{\{ε_{n}, ω_{k}\}}_{n = 0}^{N - 1}

in (13) are calculated using the tridiagonal symmetric matrix J obtained from (12) exactly as in Section 2 and shown by Eq. (4b). Additionally, the approximation (13) becomes exact if

f (x)

is a polynomial in x of a degree less than or equal to

2 N - 1

. Here too, if instead of the sum (13) we have

\sum_{k = 0}^{M} f (x_{k})

, then we can write

\sum_{k = 0}^{M} f (x_{k}) ≅ \sum_{n = 0}^{N - 1} {\tilde{ω}}_{n} f (ε_{n}),

(14)

where the derivative weights are calculated as

{\tilde{ω}}_{n} = ω_{n} / χ (ε_{n})

. As an example, we consider the function

f (x) = r^{x} / Γ (x + 1)

with

(x, r) > 0

and obtain the approximation of its infinite sum (14) using GQ associated with the Charlier and Meixner polynomials. Since

x_{k} = k = 0, 1, 2, ...

for both polynomials, then the infinite sum has the exact value

e^{r}

. Table 1 demonstrates the validity of formula (14) and illustrates the relative accuracy and convergence associated with the two polynomials for this particular function. In another example, we consider the finite sum (14) for the function

f (x) = (x + 1) r^{x + 1} / Γ (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. [10])

\sum_{k = 0}^{M} \frac{(k + 1) r^{k + 1}}{Γ (k + r + 2)} = \frac{1}{Γ (r)} - \frac{r^{M + 2}}{Γ (M + r + 2)} .

(15)

Table 2 demonstrates the validity of formula (14) 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 μ (x) = ρ (x) d x

contains a mix of discrete as well as continuous spectra, then the weight function could be written as follows

ρ (x) = σ (x) + χ (x) \sum_{k = 0}^{M} δ (x - x_{k}),

(16)

where

σ (x)

is the continuous component of the weight function and

χ (x_{k}) : = ξ_{k}

is the discrete component. The discrete spectrum (a.k.a mass points or bound states) is the set

{\{x_{k}\}}_{k = 0}^{M}

of size

M + 1

, which could be finite or infinite. The orthogonality (2) becomes [6,7,8]

\int_{x_{-}}^{x_{+}} σ (x) p_{n} (x) p_{m} (x) d x + \sum_{k = 0}^{M} ξ_{k} p_{n} (x_{k}) p_{m} (x_{k}) = δ_{n, m} .

(17)

These polynomials do also satisfy a symmetric three-term recursion relation similar to (3) 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 {\{x_{k}\}}_{k = 0}^{M}

then Gauss quadrature integral approximation (5) could be extended for this scenario to read

\int_{x_{-}}^{x_{+}} σ (x) f (x) d x + \sum_{k = 0}^{M} ξ_{k} f (x_{k}) ≅ \sum_{n = 0}^{N - 1} ω_{n} f (ε_{n}),

(18)

where the set

{\{ε_{n}, ω_{k}\}}_{n = 0}^{N - 1}

in (18) are calculated using the tridiagonal symmetric matrix J exactly as in Section 2. Here too, the integral approximation (18) becomes exact if

f (x)

is a polynomial in x of a degree less than or equal to

2 N - 1

. If the sum in (18) is finite such that M is of the same order as N (i.e.,

M \underset{\sim}{<} N

), then one could think of (18) as an integral approximation by the difference between two finite sums as follows

\int_{x_{-}}^{x_{+}} σ (x) f (x) d x ≅ \sum_{n = 0}^{N - 1} ω_{n} f (ε_{n}) - \sum_{k = 0}^{M} ξ_{k} f (x_{k}) .

(19)

As an example, we consider the function

f (x) = x^{m} e^{- x / 2}

where m is a positive integer and obtain an approximation of its integral-sum (18) using GQ associated with the continuous dual Hahn polynomial

S_{n}^{μ} (x^{2}; α, α)

with a mixed spectrum (i.e.,

μ < 0

). Since the polynomial argument is not x but

x^{2}

and

x_{k}^{2} = - {(k + μ)}^{2}

, then the GQ integral formula (18) must be revised to read

\int_{0}^{\infty} σ (x) f (x^{2}) d x + \sum_{k = 0}^{M} ξ_{k} f (x_{k}^{2}) ≅ \sum_{n = 0}^{N - 1} ω_{n} f (ε_{n}) .

(20)

Table 3 demonstrates the validity this formula and illustrates its convergence for a given μ and several values of α. The GQ result is compared to the exact result

{[f (J)]}_{00}

[11].

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

\int_{x_{-}}^{x_{+}} f (x) d x + \sum_{k = 0}^{M} f (x_{k}) ≅ \sum_{n = 0}^{N - 1} {\tilde{ω}}_{n} f (ε_{n}),

(21)

where we determine the derivative weights

{\tilde{ω}}_{n}

in terms of

ω_{n}

σ (ε_{n})

and

χ (ε_{n})

Appendix: Relevant orthogonal polynomials

In this Appendix, we give the most essential properties of orthogonal polynomials that are relevant to our work. All these polynomials belong to the hypergeometric class in the Askey scheme [9]. We consider the orthonormal version of these polynomials. That is, where the corresponding three-term recursion relation is symmetric as given by (3) and the orthogonality is a pure

δ_{n, m}

similar to those given by Eqs. (2), (11) and (17). Therefore, such polynomials are well-defined once we give the recursion coefficients

\{a_{n}, b_{n}\}

and the weight function

ρ (x)

\{ξ_{k}\}

, or the pair

\{σ (x), ξ_{k}\}

corresponding to continuous, discrete, or mixed spectra, respectively. Moreover, we always take the initial polynomials as

p_{0} (x) = 1

and

p_{1} (x) =

(x - a_{0}) / b_{0}

A.1 Discrete weight function

In this section of the Appendix, we give the essential properties of three of the discrete polynomials. The first two form infinite sequences whereas the third is a finite sequence. However, all have

x_{k} = k

The Charlier polynomial:

C_{n}^{μ} (k) = \sqrt{\frac{μ^{n}}{n!}} {}_{2}F_{0} (\begin{matrix} - n, - k \\ ¯ \end{matrix}| - 1 / μ),

(A1)

where

μ > 0

and

k, n = 0, 1, 2, ...

. The recursion coefficients are

a_{n} = n + μ, b_{n} = - \sqrt{μ (n + 1)} .

(A2)

The discrete weight function reads as follows

ξ_{k} = μ^{k} / e^{μ} k! .

(A3)

The Meixner polynomial:

M_{n}^{μ} (k; β) = \sqrt{\frac{{(2 μ)}_{n}}{n!}} β^{n / 2} {}_{2}F_{1} (\begin{matrix} - n, - k \\ 2 μ \end{matrix}| 1 - β^{- 1}),

(A4)

where

μ > 0

1 > β > 0

and

k, n = 0, 1, 2, ...

. Moreover,

{(a)}_{n} = a (a + 1) (a + 2) ... (a + n - 1) =

\frac{Γ (n + a)}{Γ (a)}

is the Pochhammer symbol (the shifted factorial). The recursion coefficients are

a_{n} = \frac{n (1 + β) + 2 μ β}{1 - β}, b_{n} = - \frac{\sqrt{β}}{1 - β} \sqrt{(n + 1) (n + 2 μ)} .

(A5)

The discrete weight function reads as follows

ξ_{k} = {(1 - β)}^{2 μ} {(2 μ)}_{k} \frac{β^{k}}{k!} .

(A6)

The Krawtchouk polynomial:

K_{n}^{M} (k; γ) = \sqrt{\frac{M!}{n! (M - n)!}} {(\frac{γ}{1 - γ})}^{n / 2} {}_{2}F_{1} (\begin{matrix} - n, - k \\ - M \end{matrix}| γ^{- 1}),

(A7)

where

1 > γ > 0

and

k, n = 0, 1, 2, .., M

. The recursion coefficients are

a_{n} = M γ + n (1 - 2 γ), b_{n} = - \sqrt{(n + 1) (M - n) γ (1 - γ)} .

(A8)

The discrete weight function reads as follows

ξ_{k} = {(1 - γ)}^{N - k} \frac{Γ (M + 1) γ^{k}}{Γ (M - k + 1) k!} .

(A9)

A.2 Discrete and continuous weight function mix

In this section of the Appendix, we give the essential properties of two orthogonal polynomials that have either a pure continuous measure (if

μ > 0

) or a mix of continuous and discrete measures (if

μ < 0

). The discrete spectrum points for both are

x_{k}^{2} = - {(k + μ)}^{2}

, which is finite and of size equal to

⌊- μ⌋

where

⌊z⌋

is the largest integer less than z.

The continuous dual Hahn polynomial:

S_{n}^{μ} (x^{2}; α, β) = \sqrt{\frac{{(μ + α)}_{n} {(μ + β)}_{n}}{n! {(α + β)}_{n}}} {}_{3}F_{2} (\begin{matrix} - n, μ + i x, μ - i x \\ μ + α, μ + β \end{matrix}| 1) .

(A10)

It is a polynomial in

x^{2}

of degree n with

x_{-} = 0

and

x_{+} = + \infty

. If

μ > 0

then α and β must be positive or a pair of complex conjugates with positive real parts. However, if

μ < 0

then

α + μ

and

β + μ

should be positive or a pair of complex conjugates with positive real parts. The recursion coefficients are

a_{n} = (n + μ + α) (n + μ + β) + n (n + α + β - 1) - μ^{2},

(A11a)

b_{n} = - \sqrt{(n + 1) (n + α + β) (n + μ + α) (n + μ + β)} .

(A11b)

The continuous component of the weight function reads as follows

σ (x) = \frac{1}{2 π} \frac{{|Γ (μ + i x) Γ (α + i x) Γ (β + i x) / Γ (2 i x)|}^{2}}{Γ (μ + α) Γ (μ + β) Γ (α + β)} .

(A12)

The discrete component of the weight function is

ξ_{k} = 2 \frac{Γ (α - μ) Γ (β - μ)}{Γ (α + β) Γ (1 - 2 μ)} \frac{(- μ - k)}{{(- 1)}^{k} k!} \frac{{(μ + α)}_{k} {(μ + β)}_{k} {(2 μ)}_{k}}{{(μ - α + 1)}_{k} {(μ - β + 1)}_{k}},

(A13)

where

k = 0, 1, 2, .., ⌊- μ⌋

The Wilson polynomial:

\begin{matrix} W_{n}^{μ} (x^{2}; ν; α, β) = \\ = \sqrt{(\frac{2 n + μ + ν + α + β - 1}{n + μ + ν + α + β - 1}) \frac{{(μ + α)}_{n} {(μ + β)}_{n} {(μ + ν)}_{n} {(μ + ν + α + β)}_{n}}{{(ν + α)}_{n} {(ν + β)}_{n} {(α + β)}_{n} n!}} {}_{4}F_{3} (\begin{matrix} - n, n + μ + ν + α + β - 1, μ + i x, μ - i x \\ μ + ν, μ + α, μ + β \end{matrix}| 1) \end{matrix}

(A14)

It is a polynomial in

x^{2}

of degree n with

x_{-} = 0

and

x_{+} = + \infty

. If

μ > 0

then

\{ν, α, β\}

should be positive or complex conjugates with positive real parts. However, if

μ < 0

then

\{ν, α, β\} + μ

must be positive or complex conjugates with positive real parts. The recursion coefficients are

a_{n} = \frac{(n + μ + ν) (n + μ + α) (n + μ + β) (n + μ + ν + α + β - 1)}{(2 n + μ + ν + α + β) (2 n + μ + ν + α + β - 1)} + \frac{n (n + ν + α - 1) (n + ν + b - 1) (n + α + b - 1)}{(2 n + μ + ν + α + β - 1) (2 n + μ + ν + α + β - 2)} - μ^{2},

(A15a)

b_{n} = - \frac{1}{2 n + μ + ν + α + β} \sqrt{\frac{(n + 1) (n + μ + ν) (n + α + β) (n + μ + α) (n + μ + β) (n + ν + α) (n + ν + β) (n + μ + ν + α + β - 1)}{(2 n + μ + ν + α + β - 1) (2 n + μ + ν + α + β + 1)}} .

(A15b)

The continuous component of the weight function reads as follows

σ (x) = \frac{1}{2 π} \frac{Γ (μ + ν + α + β) {|Γ (μ + i x) Γ (ν + i x) Γ (α + i x) Γ (β + i x) / Γ (2 i x)|}^{2}}{Γ (μ + ν) Γ (α + β) Γ (μ + α) Γ (μ + β) Γ (ν + α) Γ (ν + β)} .

(A16)

The discrete component of the weight function is

ξ_{k} = 2 \frac{Γ (μ + ν + α + β) Γ (ν - μ) Γ (α - μ) Γ (β - μ)}{Γ (- 2 μ + 1) Γ (α + β) Γ (α + ν) Γ (β + ν)} \frac{(- μ - k) {(2 μ)}_{k} {(μ + ν)}_{k} {(μ + α)}_{k} {(μ + β)}_{k}}{{(μ - ν + 1)}_{k} {(μ - α + 1)}_{k} {(μ - β + 1)}_{k} k!},

(A17)

where

k = 0, 1, 2, .., ⌊- μ⌋

References

H. Stroud and D. Secrest, Gaussian Quadrature Formulas (Prentice-Hall, 1966).
N. Kovvali, Theory and Applications of Gaussian Quadrature Methods (Springer, 2011).
H. Brass and K. Petras, Quadrature Theory: The Theory of Numerical Integration on a Compact Interval (AMS, 2011).
H., Moniem. Gaussian Quadrature for sums: A rapidly convergent summation scheme. Math. Comput. 2010, 79, 857. [Google Scholar]
S. Engblom, Gaussian Quadratures with Respect to Discrete Measures (Uppsala University, Department of information technology, Technical report 2006-007, 2006).
G. Szegő, Orthogonal Polynomials (American Mathematical Society, 1939).
T. S. Chihara, An Introduction to Orthogonal Polynomials (Dover, 2011).
M. E. H. Ismail, Classical and Quantum orthogonal polynomials in one variable (Cambridge University press, 2009).
R. Koekoek, P. A. Lesky and R. F. Swarttouw, Hypergeometric Orthogonal Polynomials and Their q-Analogues (Springer, 2010).
D. Zwillinger (ed.), Standard Mathematical Tables and Formulas, 31^st edition (CRC press, 2003).
Using the orthogonality (17) and in analogy to (4c), we can write $\int_{x_{-}}^{x_{+}} σ (x) x^{l} p_{n} (x) p_{m} (x) d x + \sum_{k = 0}^{M} ξ_{k} x_{k}^{l} p_{n} (x_{k}) p_{m} (x_{k}) = {(J^{l})}_{n, m}$ Therefore, if we approximate the function $f (x^{2})$ by a Taylor series up to ${(x^{2})}^{l}$ then $f (x^{2})$ becomes a polynomial in $x^{2}$ of degree l, which we call $F_{l} (x^{2})$ . Thus, we can write $\begin{matrix} \int ρ (x) F_{l} (x^{2}) p_{n} (x^{2}) p_{m} (x^{2}) d x = \int_{x_{-}}^{x_{+}} σ (x) F_{l} (x^{2}) p_{n} (x^{2}) p_{m} (x^{2}) d x \\ + \sum_{k = 0}^{M} ξ_{k} F_{l} (x_{k}^{2}) p_{n} (x_{k}^{2}) p_{m} (x_{k}^{2}) = {[F_{l} (J)]}_{n, m} \end{matrix}$ Hence, the left side of Eq. (20) has the exact value: $\lim_{l \to \infty} {[F_{l} (J)]}_{0, 0} = {[f (J)]}_{0, 0}$ . To calculate the matrix $F_{l} (J)$ , we proceed as follows. Let ${\{Λ_{m, n}\}}_{m = 0}^{K - 1}$ be the normalized eigenvector of the $K \times K$ finite submatrix of J corresponding to the eigenvalue $ε_{n}$ . Then, we can write $F_{l} (J) = Λ W Λ^{T}$ where W is a diagonal matrix whose elements are: $W_{n, m} =$ $δ_{n, m} F_{l} (ε_{n})$ . For improved accuracy we take the matrix size K as large as numerically possible.

Table 1. The relative error in the evaluation of the infinite sum (14) for the function

f (x) = r^{x} / Γ (x + 1)

with

r = 3.0

using Gauss quadrature associated with the Charlier and Meixner polynomials for several values of N. The polynomial parameter μ is fixed as

μ = 2.0

for both polynomials. However, we took several values for the Meixner polynomial parameter β. The relative error is calculated as

|\frac{e^{r} - Σ}{e^{r} + Σ}|

, where Σ is the Gauss quadrature sum. For numerical computations, we used Mathcad^® software version 14.0 with single precision.

Table 1. The relative error in the evaluation of the infinite sum (14) for the function

f (x) = r^{x} / Γ (x + 1)

with

r = 3.0

using Gauss quadrature associated with the Charlier and Meixner polynomials for several values of N. The polynomial parameter μ is fixed as

μ = 2.0

for both polynomials. However, we took several values for the Meixner polynomial parameter β. The relative error is calculated as

|\frac{e^{r} - Σ}{e^{r} + Σ}|

, where Σ is the Gauss quadrature sum. For numerical computations, we used Mathcad^® software version 14.0 with single precision.

GQ	$N = 2$	$N = 4$	$N = 7$	$N = 10$	$N = 15$
Charlier	5.694 $\times 10^{- 3}$	6.525 $\times 10^{- 6}$	4.165 $\times 10^{- 11}$	2.653 $\times 10^{- 16}$	8.844 $\times 10^{- 17}$
Meixner ( $β = 0.2$ )	6.943 $\times 10^{- 3}$	1.231 $\times 10^{- 4}$	1.964 $\times 10^{- 7}$	1.522 $\times 10^{- 10}$	1.946 $\times 10^{- 15}$
Meixner ( $β = 0.4$ )	3.900 $\times 10^{- 2}$	2.272 $\times 10^{- 3}$	3.192 $\times 10^{- 5}$	8.121 $\times 10^{- 7}$	1.1969 $\times 10^{- 9}$
Meixner ( $β = 0.6$ )	9.541 $\times 10^{- 2}$	5.266 $\times 10^{- 3}$	1.131 $\times 10^{- 3}$	2.588 $\times 10^{- 5}$	8.008 $\times 10^{- 6}$

Table 2. The relative error in the calculation of the finite sum (14) for the function

f (x) = (x + 1) r^{x + 1} / Γ (x + r + 2)

with

r = 3.0

using Gauss quadrature associated with the Krawtchouk polynomial for

M = 100

and for several values of N. We took several values for the polynomial parameter γ. The relative error is calculated as

|\frac{Exc - Σ}{Exc + Σ}|

, where Exc stands for the exact value obtained by formula (15) and Σ is the Gauss quadrature sum.

Table 2. The relative error in the calculation of the finite sum (14) for the function

f (x) = (x + 1) r^{x + 1} / Γ (x + r + 2)

with

r = 3.0

using Gauss quadrature associated with the Krawtchouk polynomial for

M = 100

and for several values of N. We took several values for the polynomial parameter γ. The relative error is calculated as

|\frac{Exc - Σ}{Exc + Σ}|

, where Exc stands for the exact value obtained by formula (15) and Σ is the Gauss quadrature sum.

γ	$N = 10$	$N = 20$	$N = 30$	$N = 40$	$N = 50$
0.01	4.002 $\times 10^{- 11}$	7.725 $\times 10^{- 13}$	9.770 $\times 10^{- 15}$	2.220 $\times 10^{- 16}$	5.329 $\times 10^{- 15}$
0.10	3.600 $\times 10^{- 2}$	8.826 $\times 10^{- 6}$	2.469 $\times 10^{- 11}$	6.222 $\times 10^{- 12}$	5.390 $\times 10^{- 13}$
0.20	8.514 $\times 10^{- 1}$	4.065 $\times 10^{- 2}$	1.075 $\times 10^{- 4}$	9.438 $\times 10^{- 9}$	1.799 $\times 10^{- 14}$
0.30	9.999 $\times 10^{- 1}$	6.666 $\times 10^{- 1}$	4.314 $\times 10^{- 2}$	2.807 $\times 10^{- 4}$	8.968 $\times 10^{- 8}$

Table 3. The relative error in the calculation of the integral-sum combination on the left side of Eq. (20) for the function

f (x) = x^{m} e^{- x / 2}

with

m = 3

. We used QG associated with the continuous dual Hahn polynomial

S_{n}^{μ} (x^{2}; α, α)

for

μ = - 3.5

and for several values of the parameter α. The relative error is calculated as

|\frac{Exc - Σ}{Exc + Σ}|

where Σ is the Gauss quadrature and Exc is

{[f (J)]}_{00}

with a matrix size of 200 for J.

Table 3. The relative error in the calculation of the integral-sum combination on the left side of Eq. (20) for the function

f (x) = x^{m} e^{- x / 2}

with

m = 3

. We used QG associated with the continuous dual Hahn polynomial

S_{n}^{μ} (x^{2}; α, α)

for

μ = - 3.5

and for several values of the parameter α. The relative error is calculated as

|\frac{Exc - Σ}{Exc + Σ}|

where Σ is the Gauss quadrature and Exc is

{[f (J)]}_{00}

with a matrix size of 200 for J.

$α + μ$	$N = 10$	$N = 20$	$N = 30$	$N = 50$	$N = 100$
1.0	6.752 $\times 10^{- 5}$	4.338 $\times 10^{- 7}$	1.169 $\times 10^{- 8}$	6.258 $\times 10^{- 11}$	3.594 $\times 10^{- 12}$
2.0	2.012 $\times 10^{- 3}$	2.577 $\times 10^{- 5}$	9.999 $\times 10^{- 7}$	7.667 $\times 10^{- 9}$	2.048 $\times 10^{- 12}$
3.0	1.713 $\times 10^{- 2}$	4.119 $\times 10^{- 4}$	2.255 $\times 10^{- 5}$	2.584 $\times 10^{- 7}$	1.289 $\times 10^{- 10}$
4.0	7.529 $\times 10^{- 2}$	3.168 $\times 10^{- 3}$	2.403 $\times 10^{- 4}$	4.043 $\times 10^{- 6}$	3.385 $\times 10^{- 9}$
5.0	2.197 $\times 10^{- 1}$	1.494 $\times 10^{- 2}$	1.539 $\times 10^{- 3}$	3.743 $\times 10^{- 5}$	4.938 $\times 10^{- 8}$

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Gauss Quadrature for Integrals and Sums

Abstract

1. Introduction

2. The conventional GQ: Purely continuous measure

3. GQ for a purely discrete measure

4. GQ for a measure with a mix of continuous and discrete spectra

Appendix: Relevant orthogonal polynomials

A.1 Discrete weight function

A.2 Discrete and continuous weight function mix

References

MDPI Initiatives

Important Links

Subscribe