From: Mandeep Singh on
Hello everyone,

I have been writing a code to do Legendre-Gauss Quadrature. In order to get the weights and the roots, I need to solve Legendre polynomials. The problem I am having is that after n = 6, the system produces some answers with a very small imaginary component. For n >= 7, the system does not produce the answers at all.

function [Pnx Roots weights] = symle(n)
format long
syms x
j = n;
n2 = n - 1;
j2 = n2;
Pn = ((2^j*factorial(j)))^(-1)*(subs((diff('(x^2-1)^j','x',j)),'j',n)); %Rodrigues' formula
Pnp = factor(diff(Pn,'x'));
Pnx = factor(simplify((Pn)));
Roots1 = solve((Pnx));
Roots = factor(simplify(real(Roots1)));
Pn_1 = subs(((1/(2^(j2)*factorial(j2)))*diff('(x^2-1)^(j2)','x',j2)),'j2',n2);
weights = factor(real((2)./(n*subs(Pn_1,Roots).*subs(Pnp,'x',Roots))));

Also, the way results are displayed is really bad as well. I want decimal answers or reduced fractions, but what I get is a mess with cosines and sines.

Example:

[a b c] = (symle(3))
a =
(x*(5*x^2 - 3))/2
b =
0
-15^(1/2)/5
15^(1/2)/5
c =
2^3*3^(-2)
5*3^(-2)
5*3^(-2)

This is the result at n = 3. This is exactly what I want. If you go ahead and use n = 7 or 8, the difference is huge.

Anyone know what I can do?

Thanks!