Gegenbauer {sphunif}

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

• Gegen_polyn: a matrix of size c(length(theta), length(k)) containing the evaluation of the length(k) Gegenbauer polynomials at theta.

• Gegen_coefs: a vector of size length(k) containing the coefficients b_{k, p}.

• Gegen_series: the evaluation of the truncated series expansion, a vector of size length(theta).

• Gegen_norm: the Gegenbauer norm of the truncated series, a scalar if cumulative = FALSE, otherwise a vector of size length(k).

• Gegen_polyn_2d: a 4-dimensional array of size c(length(theta_1), length(theta_2), length(k), length(m)) containing the evaluation of the length(k) * length(m) 2-dimensional Gegenbauer polynomials at the bivariate grid spanned by theta_1 and theta_2.

• Gegen_coefs_2d: a matrix of size c(length(k), length(m)) containing the coefficients b_{k, m, p}.

• Gegen_series_2d: the evaluation of the truncated series expansion, a matrix of size c(length(theta_1), length(theta_2)).

• Gegen_norm_2d: the 2-dimensional Gegenbauer norm of the truncated series, a scalar.

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)