Gauss_Legen function

Gauss--Legendre quadrature

Convenience for computing the nodes $x_k$ and weights $w_k$ of the Gauss--Legendre quadrature formula in $(a, b)$: [REMOVE_ME]$\int_a^b f(x) w(x)\,\mathrm{d}x\approx\sum_{k=1}^N w_k f(x_k).\int_a^b f(x) dx\approx\sum_{k=1}^N w_k f(x_k) [REMOVE_ME_2]$.

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.

Returns

A matrix of size c(N, 1) with the nodes $x_k$

(Gauss_Legen_nodes) or the corresponding weights $w_k$

(Gauss_Legen_weights).

Description

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

.

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 c("\n", "list("statmod")") package (Smyth, 1998), and have $10^{-15}$ precision.

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

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.