cheb2leg gives negatives for odd polynomials

[erny123]

Posted on the Julia discord: https://discourse.julialang.org/t/error-in-fasttransforms-jl-for-cheb2leg/108995

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", phikleg[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?

[MikaelSlevinsky]

It's because you're using your own second-kind Chebyshev points as opposed to chebyshevpoints(Float64, nn, Val(2)), which are decreasing. Once I change that in your second code block, the output is correct.

[MikaelSlevinsky]

See https://juliaapproximation.github.io/FastTransforms.jl/dev/generated/chebyshev/

[erny123]

@MikaelSlevinsky you're right! Thanks for the help!