I’m working with Legendre transforms and according to FastTransforms.jl, one can do a Chebyshev Transform and use the cheb2leg() function to transform from Chebyshev coefficients to Legendre coefficients.
I’ve compared this method to the pure naive transform given by wikipedia and it seems like every odd coefficient comes out to be the negative of what an actual Legendre transform should give.
I am using the Gauss-Lobatto points. The functions below calculate the Gauss-Lobatto points and weights for the Legendre and Chebyshev functions.
using FFTW
using LinearAlgebra
using Printf
using FastTransforms
function LegendrePolynomial(k::Int64, x::Float64)
if k==0
return 1.0
elseif k==1
return x
end
Lk = 0.0
Lkm2 = 1 # L_{k-2}
Lkm1 = x #L_{k-1}
for j in 2:k
Lk = (2.0*j - 1)/j *x *Lkm1 - (j-1)/j * Lkm2
Lkm2 = Lkm1
Lkm1 = Lk
end
return Lk
end
function qAndLEvaluation(N::Int64, x::Float64)
k=2
Lnm2 = 1.0
Lnm1 = x
Ln =0.0
dLn = 0.0
dLnm2 = 0.0
dLnm1 = 1.0
for j in 2:N
Ln = (2.0*j - 1)/j *x *Lnm1 - (j-1)/j * Lnm2
dLn = dLnm2 + (2.0*j-1.0)*Lnm1
Lnm2 = Lnm1
Lnm1 = Ln
dLnm2 = dLnm1
dLnm1 = dLn
#Lnsum = Lksum + Lk
end
q = dLn
#second derivative from https://en.wikipedia.org/wiki/Legendre_polynomials#Rodrigues'_formula_and_other_explicit_formulas
dq = (2*x*dLn - N*(N+1)*Ln)/(1-x^2) #ddLn
return q, dq, Ln
end
function LegendreGaussLobattoNodesAndWeights(N::Int64, nit::Int64, TOL::Float64)
x = zeros(Float64,N+1)
w = zeros(Float64,N+1)
xold = 0.0
if N ==1
x[1] = -1.0
w[1] = 1.0
x[2] = 1.0
w[2] = w[1]
else
x[1] = -1.0
w[1] = 2.0/(N*(N+1))
x[end] = 1.0
w[end] = w[1]
#xold = 0.0
for j in 1:(floor(Int64,((N+1)/2))-1)
x[j+1] = - cos((j+0.25)/N*pi - 3.0/(8.0*N*pi) /(j + 0.25))
#println(j)
#x[j+1] = -cos(pi*j/N)
for k in 0:nit
q, dq, Ln = qAndLEvaluation(N, x[j+1])
delta = -q/dq
xold = x[j+1]
x[j+1] = x[j+1] + delta
if abs(x[j+1] - xold) < TOL #abs(delta) < TOL*abs(x[j])
break
end
if k==nit-1
println("Reached Iterator Max")
end
end
q, dq, Ln = qAndLEvaluation(N, x[j+1])
x[end-j] =-x[j+1]
w[j+1] = 2.0/(N*(N+1)) / Ln^2
w[end-j] = w[j+1]
end
end
if mod(N,2) == 0
q, dq, Ln = qAndLEvaluation(N, 0.0)
x[ceil(Int64,N/2)+1] = 0.0
w[ceil(Int64,N/2)+1] = 2.0/(N*(N+1)) / Ln^2
end
return x, w
end
function ChebyshevGaussLobattoNodesAndWeights(N::Int64)
x = zeros(Float64,N+1)
w = zeros(Float64,N+1)
for j in 0:N
x[j+1] = -cos(pi*j/N)
w[j+1] = pi/(N)
end
w[1] = w[1]/2.0
w[end] = w[end]/2.0
return x, w
end
function LegendreTransform(Phi::Array{Float64}, x::Array{Float64}, w::Array{Float64})
N = length(Phi)
Phik = zeros(Float64, N)
Lk = zeros(Float64,N)
Lkk = zeros(Float64,N)
ll = 0.0
den = 0.0
for k in 1:N
@. Lk = 0.0
@. Lkk = 0.0
#calculate legendre functions at x[k]
den = 0.0
for j in 1:N
#Lk[j] = LegendrePolynomial(k-1, x[j])
ll = LegendrePolynomial(k-1, x[j])
Phik[k] = Phik[k] + w[j]*Phi[j]*ll
den = den + ll^2 *w[j]
end
Phik[k] = Phik[k]/den
end
return Phik
end
Below, here I try to do a legendre transform of a function with the second, third, fourth, and fifth Legendre polynomials:
nn = 10
xcheb, wcheb = ChebyshevGaussLobattoNodesAndWeights(nn);
xleg, wleg = LegendreGaussLobattoNodesAndWeights(nn, 10000, 1.0e-6);
phi0cheb = @. 0.5*(3*xcheb^2 - 1) + 0.5*(5*xcheb^3 - 3*xcheb)+ 1/8*(35*xcheb^4 - 30*xcheb^2 + 3) + 1/8.0*(63*xcheb^5 -70*xcheb^3 + 15*xcheb)
#phi0[end] = 0.0
phi0leg = @. 0.5*(3*xleg^2 -1) + 0.5*(5*xleg^3 - 3*xleg) + 1/8*(35*xleg^4 - 30*xleg^2 + 3) + 1/8.0*(63*xleg^5 -70*xleg^3 + 15*xleg);
#set up plans
phicopy = deepcopy(phi0cheb);
pcheb = FastTransforms.plan_chebyshevtransform(phicopy,Val(2),flags=FFTW.MEASURE);
pctol = FastTransforms.plan_cheb2leg(phicopy,normcheb=false, normleg=false);
#naive direct legendre transform
phikleg = LegendreTransform(phi0leg, xleg, wleg);
#chebyshev transform
phikcheb = pcheb*phi0cheb;
phikchebtoleg = pctol*phikcheb;
for i in 1:length(phi0kcl)
@printf("%.1e, \t\t %.1e\n", phiklegl[i], phikchebtoleg[i])
end
The results are the following:
6.9e-17, -9.5e-17 2.9e-16, 3.8e-16 1.0e+00, 1.0e+00 1.0e+00, -1.0e+00 1.0e+00, 1.0e+00 1.0e+00, -1.0e+00 9.9e-16, 7.4e-16 -2.1e-16, -3.8e-16 3.5e-16, -4.8e-17 -6.6e-16, 3.6e-16 -2.8e-16, -3.8e-16
You can see that the odd polynomial coefficients (n=odd) are the negative of the actual legendre coefficient.
Anyone have an idea if this is actually what’s happening?