Gauss_Legen {sphunif}

Gauss–Legendre quadrature

Description

Convenience for computing the nodes x_k and weights w_k of the Gauss–Legendre quadrature formula in (a, b):

\int_a^b f(x) w(x)\,\mathrm{d}x\approx\sum_{k=1}^N w_k f(x_k).

.

Usage

Gauss_Legen_nodes(a = -1, b = 1, N = 40L)

Gauss_Legen_weights(a = -1, b = 1, N = 40L)

Arguments

a, b

scalars giving the interval (a, b). Defaults to (-1, 1).

N

number of points used in the Gauss–Legendre quadrature. The following choices are supported: 5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, and 5120. Defaults to 40.

Details

For C^\infty functions, Gauss–Legendre quadrature can be very efficient. It is exact for polynomials up to degree 2N - 1.

The nodes and weights up to N = 80 were retrieved from NIST and have 10^{-21} precision. For N = 160 onwards, the nodes and weights were computed with the gauss.quad function from the statmod package (Smyth, 1998), and have 10^{-15} precision.

Value

A matrix of size c(N, 1) with the nodes x_k (Gauss_Legen_nodes) or the corresponding weights w_k (Gauss_Legen_weights).

References

NIST Digital Library of Mathematical Functions. Release 1.0.20 of 2018-09-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, and B. V. Saunders, eds. https://dlmf.nist.gov/

Smyth, G. K. (1998). Numerical integration. In: Encyclopedia of Biostatistics, P. Armitage and T. Colton (eds.), Wiley, London, pp. 3088-3095.

Examples

## Integration of a smooth function in (1, 10)

# Weights and nodes for integrating
x_k <- Gauss_Legen_nodes(a = 1, b = 10, N = 40)
w_k <- Gauss_Legen_weights(a = 1, b = 10, N = 40)

# Check quadrature
f <- function(x) sin(x) * x^2 - log(x + 1)
integrate(f, lower = 1, upper = 10, rel.tol = 1e-12)
sum(w_k * f(x_k))

# Exact for polynomials up to degree 2 * N - 1
f <- function(x) (((x + 0.5) / 1e3)^5 - ((x - 0.5)/ 5)^4 +
  ((x - 0.25) / 10)^2 + 1)^20
sum(w_k * f(x_k))
integrate(f, lower = -1, upper = 1, rel.tol = 1e-12)

## Integration on (0, pi)

# Weights and nodes for integrating
th_k <- Gauss_Legen_nodes(a = 0, b = pi, N = 40)
w_k <- Gauss_Legen_weights(a = 0, b = pi, N = 40)

# Check quadrature
p <- 4
psi <- function(th) -sin(th / 2)
w <- function(th) sin(th)^(p - 2)
integrate(function(th) psi(th) * w(th), lower = 0, upper = pi,
          rel.tol = 1e-12)
sum(w_k * psi(th_k) * w(th_k))

# Integral with Gegenbauer polynomial
k <- 3
C_k <- function(th) drop(Gegen_polyn(theta = th, k = k, p = p))
integrate(function(th) psi(th) * C_k(th) * w(th), lower = 0, upper = pi,
          rel.tol = 1e-12)
th_k <- drop(Gauss_Legen_nodes(a = 0, b = pi, N = 80))
w_k <- drop(Gauss_Legen_weights(a = 0, b = pi, N = 80))
sum(w_k * psi(th_k) * C_k(th_k) * w(th_k))