sphunif (version 1.4.3)

Gauss_Legen: 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)

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).

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.

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))