# Error in FastTransforms.jl for cheb2leg

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?

1 Like

I posted on the github page:

And MikaelSlevinsky pointed out that the guass-lobatto points for the FastTransforms.jl are set from positive max to decreasing. So, using my code, I would have to reverse the chebyshev points array:

`xcheb = xcheb[end:-1:1]`

or use the FastTransforms function:

`xcheb = chebyshevpoints(Float64, nn+1, Val(2))`

This gives the correct legendre transform.