Gegenbauer {sphunif}R Documentation

Gegenbauer polynomials and coefficients

Description

The Gegenbauer polynomials \{C_k^{(\lambda)}(x)\}_{k = 0}^\infty form a family of orthogonal polynomials on the interval [-1, 1] with respect to the weight function (1 - x^2)^{\lambda - 1/2}, for \lambda > -1/2, \lambda \neq 0. They usually appear when dealing with functions defined on S^{p-1} := \{{\bf x} \in R^p : ||{\bf x}|| = 1\} with index \lambda = p / 2 - 1.

The Gegenbauer polynomials are somehow simpler to evaluate for x = \cos(\theta), with \theta \in [0, \pi]. This simplifies also the connection with the Chebyshev polynomials \{T_k(x)\}_{k = 0}^\infty, which admit the explicit expression T_k(\cos(\theta)) = \cos(k\theta). The Chebyshev polynomials appear as the limit of the Gegenbauer polynomials (divided by \lambda) when \lambda goes to 0, so they can be regarded as the extension by continuity of \{C_k^{(p/2 - 1)}(x)\}_{k = 0}^\infty to the case p = 2.

For a reasonably smooth function \psi defined on [0, \pi],

\psi(\theta) = \sum_{k = 0}^\infty b_{k, p} C_k^{(p/2 - 1)}(\cos(\theta)),

provided that the coefficients

b_{k, p} := \frac{1}{c_{k, p}} \int_0^\pi \psi(\theta) C_k^{(p/2 - 1)}(\cos(\theta)) (\sin(\theta))^{p - 2}\,\mathrm{d}\theta

are finite, where the normalizing constants are

c_{k, p} := \int_0^\pi (C_k^{(p/2 - 1)}(\cos(\theta)))^2 (\sin(\theta))^{p - 2} \,\mathrm{d}\theta.

The (squared) "Gegenbauer norm" of \psi is

\|\psi\|_{G, p}^2 := \int_0^\pi \psi(\theta)^2 C_k^{(p/2 - 1)}(\cos(\theta)) (\sin(\theta))^{p - 2}\,\mathrm{d}\theta.

The previous expansion can be generalized for a 2-dimensional function \psi defined on [0, \pi] \times [0, \pi]:

\psi(\theta_1, \theta_2) = \sum_{k = 0}^\infty \sum_{m = 0}^\infty b_{k, m, p} C_k^{(p/2 - 1)}(\cos(\theta_1)) C_k^{(p/2 - 1)}(\cos(\theta_2)),

with coefficients

b_{k, m, p} := \frac{1}{c_{k, p} c_{m, p}} \int_0^\pi\int_0^\pi \psi(\theta_1, \theta_2) C_k^{(p/2 - 1)}(\cos(\theta_1)) C_k^{(p/2 - 1)}(\cos(\theta_2)) (\sin(\theta_1))^{p - 2} (\sin(\theta_2))^{p - 2}\,\mathrm{d}\theta_1\,\mathrm{d}\theta_2.

The (squared) "Gegenbauer norm" of \psi is

\|\psi\|_{G, p}^2 := \int_0^\pi\int_0^\pi \psi(\theta_1, \theta_2)^2 C_k^{(p/2 - 1)}(\cos(\theta_1)) C_k^{(p/2 - 1)}(\cos(\theta_2)) (\sin(\theta_1))^{p - 2} (\sin(\theta_2))^{p - 2} \,\mathrm{d}\theta_1\,\mathrm{d}\theta_2.

Usage

Gegen_polyn(theta, k, p)

Gegen_coefs(k, p, psi, Gauss = TRUE, N = 320, normalize = TRUE,
  only_const = FALSE, tol = 1e-06, ...)

Gegen_series(theta, coefs, k, p, normalize = TRUE)

Gegen_norm(coefs, k, p, normalize = TRUE, cumulative = FALSE)

Gegen_polyn_2d(theta_1, theta_2, k, m, p)

Gegen_coefs_2d(k, m, p, psi, Gauss = TRUE, N = 320, normalize = TRUE,
  only_const = FALSE, tol = 1e-06, ...)

Gegen_series_2d(theta_1, theta_2, coefs, k, m, p, normalize = TRUE)

Gegen_norm_2d(coefs, k, m, p, normalize = TRUE)

Arguments

theta, theta_1, theta_2

vectors with values in [0, \pi].

k, m

vectors with the orders of the Gegenbauer polynomials. Must be integers larger or equal than 0.

p

integer giving the dimension of the ambient space R^p that contains S^{p-1}.

psi

function defined in [0, \pi] and whose Gegenbauer coefficients are to be computed. Must be vectorized. For Gegen_coefs_2d, it must return a matrix of size c(length(theta_1), length(theta_2)).

Gauss

use a Gauss–Legendre quadrature rule of N nodes in the computation of the Gegenbauer coefficients? Otherwise, call integrate. Defaults to TRUE.

N

number of points used in the Gauss–Legendre quadrature for computing the Gegenbauer coefficients. Defaults to 320.

normalize

consider normalized coefficients (divided by c_{k, p})? Defaults to TRUE.

only_const

return only the normalizing constants c_{k, p}? Defaults to FALSE.

tol

tolerance passed to integrate's rel.tol and abs.tol if Gauss = FALSE. Defaults to 1e-6.

...

further arguments to be passed to psi.

coefs

for Gegen_series and Gegen_norm, a vector of coefficients b_{k, p} with length length(k). For Gegen_series_2d and Gegen_norm_2d, a matrix of coefficients b_{k, m, p} with size c(length(k), length(m)). The order of the coefficients is given by k and m.

cumulative

return the cumulative norm for increasing truncation of the series? Defaults to FALSE.

Details

The Gegen_polyn function is a wrapper to the functions gegenpoly_n and gegenpoly_array in the gsl-package, which they interface the functions defined in the header file gsl_sf_gegenbauer.h (documented here) of the GNU Scientific Library.

Note that the function Gegen_polyn computes the regular unnormalized Gegenbauer polynomials.

For the case p = 2, the Chebyshev polynomials are considered.

Value

References

Galassi, M., Davies, J., Theiler, J., Gough, B., Jungman, G., Alken, P., Booth, M., and Rossi, F. (2009) GNU Scientific Library Reference Manual. Network Theory Ltd. http://www.gnu.org/software/gsl/

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/

Examples

## Representation of Gegenbauer polynomials (Chebyshev polynomials for p = 2)

th <- seq(0, pi, l = 500)
k <- 0:3
old_par <- par(mfrow = c(2, 2))
for (p in 2:5) {
  matplot(th, t(Gegen_polyn(theta = th, k = k, p = p)), lty = 1,
          type = "l", main = substitute(p == d, list(d = p)),
          axes = FALSE, xlab = expression(theta), ylab = "")
  axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
       labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
  axis(2); box()
  mtext(text = expression({C[k]^{p/2 - 1}}(cos(theta))), side = 2,
        line = 2, cex = 0.75)
  legend("bottomleft", legend = paste("k =", k), lwd = 2, col = seq_along(k))
}
par(old_par)

## Coefficients and series in p = 2

# Function in [0, pi] to be projected in Chebyshev polynomials
psi <- function(th) -sin(th / 2)

# Coefficients
p <- 2
k <- 0:4
(coefs <- Gegen_coefs(k = k, p = p, psi = psi))

# Series
plot(th, psi(th), type = "l", axes = FALSE, xlab = expression(theta),
      ylab = "", ylim = c(-1.25, 0))
axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
     labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
axis(2); box()
col <- viridisLite::viridis(length(coefs))
for (i in seq_along(coefs)) {
  lines(th, Gegen_series(theta = th, coefs = coefs[1:(i + 1)], k = 0:i,
                         p = p), col = col[i])
}
lines(th, psi(th), lwd = 2)

## Coefficients and series in p = 3

# Function in [0, pi] to be projected in Gegenbauer polynomials
psi <- function(th) tan(th / 3)

# Coefficients
p <- 3
k <- 0:10
(coefs <- Gegen_coefs(k = k, p = p, psi = psi))

# Series
plot(th, psi(th), type = "l", axes = FALSE, xlab = expression(theta),
      ylab = "", ylim = c(0, 2))
axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
     labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
axis(2); box()
col <- viridisLite::viridis(length(coefs))
for (i in seq_along(coefs)) {
  lines(th, Gegen_series(theta = th, coefs = coefs[1:(i + 1)], k = 0:i,
                         p = p), col = col[i])
}
lines(th, psi(th), lwd = 2)

## Surface representation

# Surface in [0, pi]^2 to be projected in Gegenbauer polynomials
p <- 3
psi <- function(th_1, th_2) A_theta_x(theta = th_1, x = cos(th_2),
                                      p = p, as_matrix = TRUE)

# Coefficients
k <- 0:20
m <- 0:10
coefs <- Gegen_coefs_2d(k = k, m = m, p = p, psi = psi)

# Series
th <- seq(0, pi, l = 100)
col <- viridisLite::viridis(20)
old_par <- par(mfrow = c(2, 2))
image(th, th, A_theta_x(theta = th, x = cos(th), p = p), axes = FALSE,
      col = col, zlim = c(0, 1), xlab = expression(theta[1]),
      ylab = expression(theta[2]), main = "Original")
axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
     labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
axis(2, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
     labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
box()
for(K in c(5, 10, 20)) {
  A <- Gegen_series_2d(theta_1 = th, theta_2 = th,
                       coefs = coefs[1:(K + 1), ], k = 0:K, m = m, p = p)
  image(th, th, A, axes = FALSE, col = col, zlim = c(0, 1),
        xlab = expression(theta[1]), ylab = expression(theta[2]),
        main = paste(K, "x", m[length(m)], "coefficients"))
  axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
       labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
  axis(2, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
       labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
  box()
}
par(old_par)
[Package sphunif version 1.4.0 Index]