# Error in FastTransforms.jl for cheb2leg

**URL:** https://discourse.julialang.org/t/error-in-fasttransforms-jl-for-cheb2leg/108995
**Category:** General Usage
**Tags:** fftw, math, fft
**Created:** [January 19, 2024, 2:28am UTC](https://discourse.julialang.org/t/error-in-fasttransforms-jl-for-cheb2leg/108995 "2024-01-19T02:28:56Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![erny123](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erny123/32/52255_2.png) [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [January 19, 2024, 2:28am UTC](https://discourse.julialang.org/t/error-in-fasttransforms-jl-for-cheb2leg/108995/1 "2024-01-19T02:28:56Z")

</div>

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](https://en.wikipedia.org/wiki/Finite_Legendre_transform) 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.

```julia
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

```

```julia
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

```

```julia
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

```

```julia
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

```

```julia
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:

```julia
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?

---

<div class="post-metadata">

### Author: ![erny123](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erny123/32/52255_2.png) [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [January 22, 2024, 10:20pm UTC](https://discourse.julialang.org/t/error-in-fasttransforms-jl-for-cheb2leg/108995/2 "2024-01-22T22:20:31Z")

</div>

I posted on the github page:

> <https://github.com/JuliaApproximation/FastTransforms.jl/issues/241>
>
> Posted on the Julia discord: \[https://discourse.julialang.org/t/error-in-fasttra…nsforms-jl-for-cheb2leg/108995\](https://discourse.julialang.org/t/error-in-fasttransforms-jl-for-cheb2leg/108995)
> 
> I’m working with Legendre transforms and according to \[FastTransforms.jl\](https://juliahub.com/ui/Packages/FastTransforms), 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\](https://en.wikipedia.org/wiki/Finite\_Legendre\_transform) 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?

And **[MikaelSlevinsky](https://github.com/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.