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.