Search for the roots of the Legendre polynomial in python

Question

I am writing a program that solves the Legendre-Gauss integral. The nth order quadrature algorithm requires, at some point, finding the roots of the nth order Legendre polynomial, Pn (x), assigning them to the Absc array (for the “abscissa”). Pn is an nth-order polynomial with n independent real roots on the interval [-1,1]. I would like to be able to calculate the roots, and not just import them from any library. I can create an array that gives polynomial coefficients, which I call PCoeff. To find the roots I tried

Absc = numpy.roots(PCoeff) 

This works until about n = 40, but beyond that it starts to fail, giving complex roots when it really shouldn't be. I also tried to define a polynomial using

 P = numpy.poly1d(PCoeff) Absc = Pr 

but this gives the same problems, apparently because it uses the same numpy root search algorithm.

Another method that seemed promising was to use scipy.optimize.fsolve (Pn, x0), where x0 is an array of n-elements of my guessing in the roots. The problem is that depending on my x0 options, this method can produce one specific root more than once instead of other roots. I tried filling x0 as equidistant points on [-1,1]

 x0 = numpy.zeros(n) step = 2./(n-1) for i in xrange(n): x0[i] = -1. + i*step 

but as soon as I get to n = 5, fsolve gives some repeating roots and neglects others. I also tried using the results of numpy.roots as x0. However, in the problem area where np.roots gives complex values, this causes an error in fsolve

 TypeError: array cannot be safely cast to required type 

I saw online that there is a scipy.optimize.roots () routine that can work, but it is not in the scipy library on my computer. Updating is a problem because I do not have permission to download things on this computer.

I would like to be able to run the quadrature in 64th order for high accuracy, but this root finds the cause of the failure. Any ideas?

Answer 1

Since np.roots rely on "searching for eigenvalues ​​of the Companion matrix", as indicated in the documentation, you may run into the problem of error propagation, which leads to a nonzero imaginary part of the roots. Perhaps you can just drop the imaginary part using the np.real function.

You can try to calculate the roots differently using the Taylor approximation of the roots:

https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial

Answer 2

You can find a simple example of implementing your SymPy in its documentation :

 >>> for n in range(5): ... nprint(polyroots(taylor(lambda x: legendre(n, x), 0, n)[::-1])) ... [] [0.0] [-0.57735, 0.57735] [-0.774597, 0.0, 0.774597] [-0.861136, -0.339981, 0.339981, 0.861136] 

An example, as the previous answer suggests, uses the taylor polynomial decomposition.