# Efficient generation of associated Legendre Polynomials

For the calculation of a high order spherical harmonics expansion model, I would need a matrix of associated Legendre Polynomials ${P}_{l}^{m}(x)$ with a constant $x$ and $\{l,m\}$ running from $\{\{2...lmax\}, \{0...l\}\}$ with $lmax$ in the order of 2000. Thus just using the straightforward approach to use Table

Table[LegendreP[l,m,Cos[37. \[Degree]]],{l,2,2000},{m,0,l}]


takes a few minutes (here I used $x=\cos(37°)$ as an example to get numerical results). However since there are recurrence relations for associated Legendre Polynomials like $$(l-m+1){P}_{l+1}^{m}(x)=(2l+1)x {P}_{l}^{m}(x)-(l+m){P}_{l-1}^{m}(x)$$ and $$\sqrt{1-x^2}{P}_{l}^{m+1}(x)=(l-m)x{P}_{l}^{m}(x)-(l+m){P}_{l-1}^{m}(x)$$ which would make the calculation of a large matrix for a defined $x$ very efficient, I was trying to use the built in function RecurrenceTable in Mathematica

x=Cos[37. \[Degree]];
RecurrenceTable[{p[l, m + 1] ==1/Sqrt[1 - x^2] ((l - m)
x p[l, m] - (l + m) p[l - 1, m]), p[l + 1, m] ==
1/(l - m + 1) ((2 l + 1) x p[l, m] - (l + m) p[l - 1, m]),
p[1, 0] == LegendreP[1, 0, x],
p[2, 0] == LegendreP[2, 0, x]}, p, {l, 2, 2000}, {m, 0, l}]


but this function throws the error that the system is overdetermined, which I believe it is not.

I did not try out NestList so far since the two-dimensional parameter space $\{l,m\}$ makes it hard to implement.

Any ideas?

You can make use of functions, which remember their values. Something like

p[l_, m_] := (p[l, m] =
1/Sqrt[1 -
x^2] ((l - (m - 1)) x p[l, (m - 1)] - (l + (m - 1)) p[
l - 1, (m - 1)])) /; (l > 0) && m > 0;

p[l_, m_] := (p[l, m] =
1/((l - 1) - m +
1) ((2 ( l - 1) + 1) x p[l - 1, m] - ((l - 1) + m) p[l - 2,
m])) /; l > 2;

p[1, 0] := LegendreP[1, 0, x];
p[1, 1] := LegendreP[1, 1, x];
p[2, 0] := LegendreP[2, 0, x];

Table[p[l, m], {l, 2, 2000}, {m, 0, l}]; // AbsoluteTiming


{59.8004, Null}

Note, however, that numerically I would not expect this to give the same results, because of precision issues. Please also check if I correctly subtracted all indices.

• Thanks for the hint, that did the job. By the way I updated the recursion I mentioned in my question with a much faster and better one taken from Mohlenkamp (ohio.edu/people/mohlenka/research/uguide.pdf) – Rainer Dec 24 '16 at 11:53
• @Rainer Could you like to share your code here? It is very useful to calculate Legendre functions numerical in an efficient way by Mathematica. – tanghe2014 May 22 '17 at 10:11

(This is a comment that got too long for the comment box, but I might edit this later with some sundry code.)

First things first: please at least try to read Gautschi's classic paper if you'll be in the business of implementing three-term recurrences for special functions. For the specific case of the associated Legendre functions/Ferrers functions, these two papers by Gil and Segura are of interest.

Briefly put: the recurrences for the regular solutions $P$ are stable for increasing $\ell$ and decreasing $m$ (the reverse is true for the irregular solutions $Q$). (In the parlance of recurrences: $P$ is dominant and $Q$ is minimal for the recurrence on the degree $\ell$, while $Q$ is dominant and $P$ is minimal for the recurrence on the order $m$.) So, generating $P$ by starting from $P_{0}^{m}(x)$ and $P_{1}^{m}(x)$ and going up is fine, but one cannot get accurate results by recursing upward from $P_{\ell}^{0}(x)$ and $P_{\ell}^{1}(x)$.

You might also consider incorporating any desired normalization factors into your recurrence at the outset.