https://github.com/cran/pracma
Tip revision: 7e2814c76ec9bc4503af661775ea73f11ec16537 authored by Hans W. Borchers on 15 November 2014, 00:00:00 UTC
version 1.7.9
version 1.7.9
Tip revision: 7e2814c
chebyshev.R
##
## c h e b y s h e v . R
##
chebPoly <- function(n, x = NULL) {
stopifnot(is.numeric(n), length(n) == 1,
floor(n) == ceiling(n), n >= 0)
N <- max(2, n+1)
T <- matrix(0, N, N)
# Preset degree 0 and 1
T[1, 1] <- 1
T[2, 2] <- 1
# Use recursion formula
if (n >= 2)
for (i in 3:N)
T[i, ] <- 2 * c(0, T[i-1, 1:(N-1)]) - T[i-2, ]
# Reverse with highest coefficient first
T <- T[ ,ncol(T):1]
if (is.null(x)) {
return(T)
} else if (is.numeric(x)) {
return(polyval(T[n+1, ], x))
} else
stop("Argument 'x' must be a numeric vector.")
}
chebCoeff <- function(fun, a, b, n) {
N <- n+1
# Map interval [a, b] to [-1, 1]
c1 <- 0.5*(b-a)
c2 <- 0.5*(b+a)
# Evaluate function at Chebyshev points (in [a, b])
k <- c(0:(N-1))
y <- cos(pi*(k+0.5)/N)
fy <- fun(c1*y + c2)
# Now compute the Chebyshev coefficients
c0 <- 2.0 / N
K <- matrix(k+0.5, N, 1)
cheb <- c0 * fy %*% cos(pi/N * K %*% k)
# Remove too small coefficients and return
eps <- 2 * .Machine$double.eps # machine precision
cheb <- ifelse(abs(cheb) < eps, 0, cheb)
return(drop(cheb))
}
chebApprox <- function(x, fun, a, b, n) {
# Compute the Chebyshev polynomials of fun in [a, b]
cP <- chebPoly(n)
cC <- chebCoeff(fun, a, b, n)
p <- drop(cC %*% cP)
c0 <- cC[1]
# Map x into [-1, 1] and evaluate the Chebyshev polynomial
xx <- (2*x - (b+a))/(b-a)
yy <- polyval(p, xx) - c0/2
return(yy)
}