Gaussian Quadrature ( Legendre Polynomials )

I have tried to create a function that will find the coefficients of the nth order Legendre polynomial without using syms x, but I have got stuck on how to actually get the coefficients with 2 unknowns in my equation... Can anyone help a bit??

Here's my code :

function [y] = LegendrePoly(n)
if n == 0
    y = 1;
elseif n == 1
    y = x;
else 
    x = zeros(n+1,1);
    x(n+1) = 1;
    x_1 = zeros(n+1,1);
    x_1(n) = 1;
    for i=2:n;
        y(i)=((2*x_1(n))/x_1(n))*x*y(i-1)-(x(n-1)/x_1(n))*y(i-2);
    end
end

[Karan Gill]

There's an example on this topic: https://www.mathworks.com/help/symbolic/examples/gauss-laguerre-quadrature-evaluation-points-and-weights.html. Is it helpful?

Answer 1

[John D'Errico]

NO NO NO. You say that you don't want to use syms here. So why in the name of god and little green apples are you doing things like this?

y = x;

All you need to store are the COEFFICIENTS!!!!!!! It is a polynomial. You have not passed x into the function to evaluate. That you can do trivially with polyval anyway.

Set up LegendrePoly to return a vector of length n+1. The coefficients of that polynomial.

I'm feeling too lazy to crack open Abramowitz & Stegun, so this page gives the recurrence relation and the first two polynomials.

https://en.wikipedia.org/wiki/Legendre_polynomials

Thus, P_0 = 1. P_1 = x, and:

(n+1)*P_(n+1) = (2*n+1)*P_n*x - n*P_(n-1)

Next, DON'T EVER write recurrence relations like this recursively, unless you use a memoized form. The problem is that the number of recursive calls grows exponentially. It is far simpler to just use a loop anyway.

function coefs = LegendrePoly(n)
if n == 0
  coefs = 1;
elseif n == 1
  coefs = [1 0];
else
  P_nm1 = 1;
  P_n = [1 0];

    for i=1:(n-1);
      P_np1 = ((2*i+1)*[P_n,0] - i*[0,0,P_nm1])/(i+1); % recurrence
      [P_nm1,P_n] = deal(P_n,P_np1); % shift
    end
    coefs = P_np1;
  end

Think about what [P_n,0] does, in terms of polynomial coefficients.

Thus we know that P_2=(3*x^2-1)/2. For n==2, we get

coefs =
  1.5   0  -0.5

Likewise, for n==5, we get

coefs
coefs =
   7.875     0    -8.75    0    1.875   0

From the table of polynomials, I got what I expected.

[63/8 0 -70/8 0 15/8 0]
ans =
   7.875     0    -8.75    0    1.875    0

You can evaluate that 5th degree polynomial simply as:

polyval(coefs,0.5)
ans =
   0.08984375

Finally, in order to use them as polynomials for Gaussian quadrature, you will need the derivative polynomials too.

polyder(coefs)
ans =
   39.375     0    -26.25      0     1.875

Its been a while since I had to derive the Gaussian quadrature but you need some roots too. Again, trivial.

roots(coefs)
ans =
                       0
      -0.906179845938664
      -0.538469310105683
       0.906179845938664
       0.538469310105683

[Ho Nam Ernest Yim]

THANK you so much