Possible improvement for the quadrature approximation routines.

[Flintro]

In the code ported from the old MATLAB libraries referenced here the Gauss-Legendre Quadrature approximation might have a few issues. When making this approximation we generate n legendre polynomials. There are a few issues.

From what the comments in the code suggest, for i = 1,...,n the function qnwlege(n::Int, a::Real, b::Real) takes its nodes to be the ith roots of the ith legendre polynomial, whereas the optimal approximation would take the n roots of the nth legendre polynomial. This happens in the lines:

        # p1 is now a vector of Legendre polynomials of degree 1..n
        # pp will be the deriative of each p1 at the nth zero of each
        pp = n*(z.*p1-p2)./(z.*z-1)
        z1 = z
        z = z1 - p1./pp  # newton's method

I think the actual code is fine and we just need to change the comment.

Another issue is that on my machine the fix() function used is not recognised, and I'm not sure what it is supposed to do. I replaced it with the ceil() function and it seems to work.

The final and perhaps most difficult problem is that for many choices of n we hit max iteration limits due to the implementation of newton's method. I have attempted to address this with a different approach below, but I sometimes get wrong answers due to the primitive nature of what I've done.

My approach was something like:

function GLQ_Nodes(n, steps)

    # generate a matrix whose rows will be points of the ith legendre polynomial
    x = collect(linspace(-1, 1, steps))
    P = Matrix(n+1,steps)
    P[1,] = ones(steps)
    P[2,] = x

   # create the ith legendre polynomials and put them in the matrix
    for i in 3:n+1
        P[i,] =  ((2*(i-1)+1)*x.*P[i-1,]-(i-1)*P[i-2,])./(i)
    end
   
   # take the last polynomial and find its roots
    nth_legendre = P[n+1,]
    y = []
    abs_legendre = abs.(nth_legendre)
    for i in 1:n
        y = append!(y, indmin(abs_legendre))
        abs_legendre = deleteat!(abs_legendre, indmin(abs_legendre))
    end

    return x[y]

end

GLQ_Nodes(10,1000)

If someone knows how to fix the max iteration problem then that would be great. Otherwise I'll try to fix my approach.