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

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.

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.

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

Run this code

## 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)
# \donttest{
## 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)
# }

Run the code above in your browser using DataLab