I've been triying several packages to compute the orthonormal Gegenbauer polynomials basis which seems to be orthonormal at [-1,1]. Using the package midasmlI get however a different result:
library(midasml)
x <- seq(-1, 1, length.out = 100)
matplot(B <- gb(degree = 3, alpha = 1, a = -1, b = 1, jmax = NULL, X = x), type="l")
round(t(B) %*% B, 3) #not orthonormal
any reasonement why this is happening? any suggestions for an alternative approach? Thanks!
Orthonormality doesn't mean that this matrix is orthogonal. It means that the integral from -1 to 1 of P_i(x)P_j(x)w(x) is 0 if i != j and 1 if i = j, where w is a weight function. The weight function for the Gegenbauer polynomials with parameter alpha is (1 - x^2)^(alpha - 1/2).
You can get these polynomials with the orthopolynom package.
library(orthopolynom)
alpha <- 1
gp.list <- gegenbauer.polynomials(3, alpha, normalized = FALSE)
print(gp.list)
# [[1]]
# 1
#
# [[2]]
# 2*x
#
# [[3]]
# -1 + 4*x^2
#
# [[4]]
# -4*x + 8*x^3
w <- function(x) { # the weight function
(1 - x^2)^(alpha - 1/2)
}
GP2 <- function(x) {
-1 + 4*x^2
}
GP3 <- function(x) {
-4*x + 8*x^3
}
f <- function(x) {
GP2(x) * GP3(x) * w(x)
}
integrate(f, lower = -1, upper = 1) # should be zero
# 0 with absolute error < 1.4e-14
f <- function(x) {
GP3(x) * GP3(x) * w(x)
}
integrate(f, lower = -1, upper = 1) # should be pi/2
# 1.570796 with absolute error < 6e-07
Here the integral of P_i(x)²w(x) is not 1 because I took normalized = FALSE.