Legendre polynomials in MatLab

So, due to a HW assignment I work on, I needed to have a fast code that computes the nth degree Legendre polynomial. Looking around, I found out that Matlab has a built-in function, y = legendre(n,x), that calculates the value of  associated Legendre polynomials up to n for x. If x is a scalar, y is a vector. If x is a vector, y is a matrix and so on. I only needed the first elements each time (because those corresponded to the Legendre polynomial).

Now, that was a huge waste of both time and memory resources for me, because I wanted to calculate a lot of integrals and I didn’t have any use for the extra values returned. I then looked for something else at Matlab Central, but all the algorithms there for computing the coefficients were unstable. Take n large enough and we diverge to infinity as we get to 1 (instead of having the value of 1).

Looking at one of my favorites textbooks, Numerical Recipes in C, I found a short code that does the same thing as the Legendre function in Matlab. I translated it to Matlab and it works really well, so here it is.

function pn = legPoly(x,i)
% LEGPOLY calculates the legendre polynomial of i-th order in x
% x can be a scalar, a vector or a matrix

% pn : nth legendre polynomial
% pn_1 : n-1 legendre polynomial
% pn_p1 : n+1 legendre polynomial

% Konstantinos G.
% Modified version of the C code in "Numerical recipes in C"

pn_p1 = 1;
pn = ones(size(x,1), size(x,2));

if(i>0)
 pn_1 = x.*pn;
 if(i==1)
 pn = pn_1;
 else
 for i = 2 : i
 pn_p1 = (x .* (2 * i - 1).*pn_1 - (i - 1)*pn)/i;
 pn = pn_1;
 pn_1 = pn_p1;
 end
 pn = pn_p1;
 end
end
end

I don’t claim that this code is either written in the best way or that it is bug-free. From what little testing I did, it works and can help if you just want to compute the values of Legendre polynomials fast.

By the way, if want to write matlab code in your blog as I did, with all the  fancy markup, you can learn the secret at Andres Marrugo’s blog.

Advertisements

8 Comments

  1. Hi, I want to generate the first 4 legendre polynomials and then plot them and I dont know how to to them. Can you help me pls?

    1. Hello Naomi,

      I can help you but you need to answer a couple of questions first.

      1) Do you want to do it Matlab or do you have to do it in another language?

      2) Are you stuck in generating the polynomials, plotting them or both?

  2. %initialisation
    L{1} = [1 0];
    L{2} = [3/2 0 -1/2];

    %construction des poplynomes de legendre par récurrence
    for i = 2 : 1 : n
    L{i+1} = ([L{i} 0]*(2*i+1) – [0 0 L{i-1}]*i)/(i+1);
    end

    here is my matlab code to generate the legendres polynomes. (L0 is missing).
    PLotting them should not be a problem (help polyval if trouble).

    When i’m trying to calculate the roots, I got complex values for n >= 44 (for n< 44 I have the good values) …

    1. Hi Quentin,

      You are using Bonnet’s relation, the same relation that I am using in the code above. The difference is that I am using the relation to calculate the value at a particular x and you are using it to compute the coefficients of the nth polynomial. This is, unfortunately, an unstable way to do it and it breaks down for large n.

      Here’s an easy, first test that you can try to see if you are on the right track. Remember that Legendre polynomials are normalized, so that P_{n}(1)=1 for any n. If your polynomials don’t become 1 at 1, then something’s wrong.

      I am guessing that if you change your code so as to return the value of the nth polynomial instead of returning the coefficients, that you will get the correct values.

  3. Hello sir,
    I am doing project on Legendre Moment. I am unable to applying it on any gray scale image of let size N*N. I am doing it in MATLAB. Please help me .

    1. Hi i can help you but you need to answer a couple of questions.

      1) Why you use Legendre Polynomials?
      2) What is your purpose?

      1. Sir as i mentioned i m doing my project in which i have to work on Legendre Moment which is calculated by using Legendre function. But Iam new to MATLAB and this Topic also and not get it how to do this. Please sir help me.

Comments are closed.