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?


2 Answers 2


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.

  • $\begingroup$ 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) $\endgroup$
    – Rainer
    Dec 24, 2016 at 11:53
  • $\begingroup$ @Rainer Could you like to share your code here? It is very useful to calculate Legendre functions numerical in an efficient way by Mathematica. $\endgroup$
    – tanghe2014
    May 22, 2017 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.