Re: LegendreP error (bug?) in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg81423] Re: LegendreP error (bug?) in Mathematica*From*: Roman <rschmied at gmail.com>*Date*: Sun, 23 Sep 2007 04:27:58 -0400 (EDT)*References*: <fcb3qh$fk1$1@smc.vnet.net><fco79o$3p9$1@smc.vnet.net>

Here's a stable and compiled version for spherical harmonics and associated Legendre polynomials. This works for |x|<=1 and |M|<=L, with L and M integers. Thanks to Oleksandr for input on the user interface. As before, Y[L, M, x] numerically computes SphericalHarmonicY[L, M, ArcCos[x], 0], and P[L, M, x] computes LegendreP[L, M, x]. If you need SphericalHarmonicY[L, M, theta, phi], then just use Y[L, M, ArcCos[theta]] * Exp[I*M*phi] instead, for better speed and accuracy. Y[L,M,x] is computed iteratively as Y[M,M,x], Y[M+1,M,x], Y[M +2,M,x], ..., Y[L-1,M,x], Y[L,M,x]. Therefore, if you need all these spherical harmonics for given M and x, then you'd be much better off rewriting this algorithm to spit out all intermediate results, instead of re-doing the whole iteration for each value of L! (unless timing is of no importance). Some coefficients are pre-computed once and for all, and thus we need to specify an upper limit on |M| of what will ever be needed. I'm sorry about this; if you have a more elegant way of dealing with this, please let us know! Mmax = 1000; precompute coefficients: this cannot be done with machine numbers, because the factorials become too large for large indices. YCh = Table[(-1)^L(2L-1)!!Sqrt[(2L+1)/(4Pi(2L)!)],{L,0,Mmax}]//N; compiled code for numerically computing spherical harmonics from recursion relation: YC = Compile[{{L,_Integer},{M,_Integer},{x,_Real}}, Block[{A}, Which[ M>L,0., M==L,YCh[[L+1]](1-x^2)^(L/2), M==L-1,Sqrt[2L+1]YCh[[L]]x(1-x^2)^((L-1)/2), True,A=YCh[[M+1]](1-x^2)^(M/2); Nest[{#[[1]]+1,#[[3]], Sqrt[(2#[[1]]+1)/((#[[1]]-M)(#[[1]]+M))] (Sqrt[2#[[1]]-1]x#[[3]] -Sqrt[((#[[1]]-1+M)(#[[1]]-1-M))/ (2#[[1]]-3)]#[[2]])}&, {M+2,A,Sqrt[2M+3]x A}, L-M-1][[3]]]], {{YCh,_Real,1}}]; user interface, which does some range checking: (defaults to built- in algorithm if conditions are not fulfilled) Y[L_Integer?NonNegative,M_Integer?NonNegative, x_Real/;MachineNumberQ[x]&&-1<=x<=1]/; Developer`MachineIntegerQ[L]&& Developer`MachineIntegerQ[M]&& M<=Mmax := YC[L,M,x] Y[L_Integer?NonNegative,M_Integer?Negative, x_Real/;MachineNumberQ[x]&&-1<=x<=1]/; Developer`MachineIntegerQ[L]&& Developer`MachineIntegerQ[M]&& -M<=Mmax := (-1)^M YC[L,-M,x] Y[L_,M_,x_] := SphericalHarmonicY[L,M,ArcCos[x],0] compute associated Legendre polynomials from the above algorithm for spherical harmonics: P[L_,M_,x_] := Sqrt[4Pi/(2L+1)(L+M)!/(L-M)!]Y[L,M,x] Cheers! Roman. On Sep 18, 11:55 am, Roman <rschm... at gmail.com> wrote: > Oleksandr, et al., > > Your compiled code works only for small enough L and M, since the > associated Legendre polynomials quickly exceed the range of machine > numbers. It then defaults to slow evaluation. To speed things up you > can start from a stable recursion relation for the spherical harmonics > (like the GSL people do), something like this for Y[L,M,x] = > SphericalHarmonicY[L,M,ArcCos[x],0]: > > Y[L_, M_, x_] := 0/;M>L > Y[L_, M_, x_] := (-1)^M Y[L,-M,x]/;M<0 > Y[L_, M_, x_] := (-1)^L(2L-1)!!Sqrt[(2L+1)/(4Pi(2L)!)](1-x^2)^(L/ > 2)/;M==L > Y[L_, M_, x_] := x Sqrt[2L+1]Y[L-1,L-1,x]/;M==L-1 > Y[L_, M_, x_] := Nest[{#[[1]]+1,#[[3]], > Sqrt[(2#[[1]]+1)/((#[[1]]+M)(#[[1]]-M))](x > Sqrt[2#[[1]]-1]#[[3]]- > Sqrt[(#[[1]]+M-1)(#[[1]]-M-1)/(2#[[1]]-3)]#[[2]])}&, > {M+2,Y[M,M,x],Y[M+1,M,x]},L-M-1][[3]] > > If you compile this like you did, you'll get a good algorithm for the > spherical harmonics which never overflows the machine numbers. From > those you'd use arbitrary precision math to compute > P[L_, M_, x_] := Sqrt[4Pi(L+M)!/((2L+1)(L-M)!)]Y[L,M,x] > which may easily overflow the machine numbers. > > Roman. > > On Sep 18, 6:50 am, sashap <pav... at gmail.com> wrote: > > > On Sep 17, 2:32 am, Roman <rschm... at gmail.com> wrote: > > > > John, Oleksandr, Bob, Andrzej, et al., > > > > If the problem here is arithmetic cancellations, as Oleksandr claims, > > > then Mathematica is definitely not using the best algorithm out there! > > > It is true that when you write out an associated Legendre polynomial > > > explicitely and then plug in numbers, the evaluation of the polynomial > > > becomes terribly unstable because of extreme cancellations. But that's > > > NOT how you're supposed to compute Legendre polynomials; instead, you > > > can use recursion relations, which are numerically stable. > > > > Andrzej: Thanks for pointing out the issue with MachinePrecision > > > numbers lacking precision tracking. But again, this just explains > > > Mathematica's behavior but not why such a bad algorithm was chosen in > > > the first place. > > > > I agree that by upping the WorkingPrecision you can make the > > > Mathematica built-in algorithm work (otherwise something would be > > > seriously amiss), but this does not explain why the built-in algorithm > > > is so bad in the first place that it requires extra precision. Here's > > > a more stable algorithm that I've assembled from equations (66,69,70) > > > of MathWorld's Legendre page (http://mathworld.wolfram.com/ > > > LegendrePolynomial.html): (no guarantees on speed) > > > > P[L_, M_, x_] := 0 /; M > L > > > P[L_, M_, x_] := (-1)^L (2 L - 1)!! (1 - x^2)^(L/2) /; M == L > > > P[L_, M_, x_] := x (2 L - 1) P[L - 1, L - 1, x] /; M == L - 1 > > > P[L_, M_, x_] := Nest[{#[[1]] + 1, #[[3]], > > > (x (2 #[[1]] + 1) #[[3]] - (#[[1]] + M) #[[2]])/(#[[1]] + 1 - M)} > > > &, > > > {M + 1, P[M, M, x], P[M + 1, M, x]}, L - M - 1][[3]] > > > > It works only for integer L and M, and does not do any checks on the > > > parameters. You can see that this is more stable by plotting something > > > like > > > > Plot[P[200, 43, x], {x, -1, 1}] > > > Plot[LegendreP[200, 43, x], {x, -1, 1}] > > > > >From this you can calculate the spherical harmonics using equation (6) > > > > ofhttp://mathworld.wolfram.com/SphericalHarmonic.html > > > > Y[L_, M_, th_, ph_] := Sqrt[((2 L + 1)/(4 Pi))*((L - M)!/(L + M)!)] * > > > P[L, M, Cos[th]] E^(I*M*ph) > > > > Roman. > > > Roman, > > > Thank you for a very useful observation. Version 6 is using > > hypergeometric function to compute LegendreP, which is > > suboptimal for machine numbers. > > > Using Compile makes your evaluation scheme quite fast: > > > lp = Compile[{{n, _Integer}, {m, _Integer}, {z, _Real}}, > > Block[{P, L, M, x}, > > Which[m > n, 0., m == n, (-1)^m (2 m - 1.)!! (1 - z^2)^(m/2), > > m == n - 1, z (2 n - 1) (-1)^m (2 m - 1.)!! (1 - z^2)^(m/2), > > True, Nest[{#1[[1]] + 1, #1[[3]], ( > > z (2 #1[[1]] + 1) #1[[3]] - (#1[[1]] + m) #1[[2]])/(#1[[ > > 1]] + 1 - m)} &, ({m + 1, #1, z (2 m + 1) #1} &)[(-1)^ > > m (2 m - 1.)!! (1 - z^2)^(m/2)], n - m - 1][[3]]]]]; > > > Clear[P] > > P[n_Integer?NonNegative, m_Integer?NonNegative, > > z_Real /; MachineNumberQ[z] && -1 <= z <= 1] /; > > Developer`MachineIntegerQ[n] && Developer`MachineIntegerQ[m] := > > lp[n, m, z] > > P[n_, m_, z_] := LegendreP[n, m, z] > > > z is restricted to -1<=z<=1, because LegendreP grows too fast > > outside of it, so Compile will throw an exception. Try > > > lp[200, 43, -1.2] > > > With the definition above Plot[P[200,43,x],{x,-1,1}] takes > > half a second and produces reliable result. > > > Future version of Mathematica will use this method for > > numerical evaluation of LegendreP polynomials. > > > Regards, > > Oleksandr Pavlyk > > Special Functions Developer > > Wolfram Research > > > > On Sep 15, 10:09 am, sashap <pav... at gmail.com> wrote: > > > > > On Sep 14, 3:13 am, Roman <rschm... at gmail.com> wrote: > > > > > > I confirm the problem. Just as an example, > > > > > > In[1] := LegendreP[200, 43, 4/5] // N > > > > > Out[1] = 2.9256424676613492`*^97 > > > > > > In[2] := LegendreP[200, 43, 0.8] > > > > > Out[2] = 6.151579920980095`*^118 > > > > > > give strikingly different results! (The former result is accurate.) > > > > > The cause for such a discrepancy is arithmetic cancellations. > > > > You can see it by subtracting two close and sufficiently large > > > > numbers: > > > > > In[9]:= N[10^17 + 23] - N[10^17] > > > > > Out[9]= 16. > > > > > Mathematica's arbitrary precision arithmetic tells you > > > > that the result is not to be trusted: > > > > > In[10]:= N[10^17 + 23, 16] - N[10^17, 16] > > > > > Out[10]= 0.*10^1 > > > > > Something similar happens inside LegendreP algorithm. > > > > The cure is simple, try using arbitrary precision arithmetic. > > > > > Compare > > > > > Plot[LegendreP[200, 43, x], {x, 0, 1}] > > > > > vs. > > > > > Plot[LegendreP[200, 43, x], {x, 0, 1}, WorkingPrecision -> 75] > > > > > You can even do > > > > > Plot[LegendreP[200, 43, x], {x, 0, 1}, > > > > WorkingPrecision -> Infinity] > > > > > then all arguments x are rationalized (with SetPrecision[x, Infinity]) and > > > > LegendreP is evaluated exactly, and then turned into a number from > > > > plotting. > > > > > Oleksandr Pavlyk > > > > Special Functions developer > > > > Wolfram Research > > > > > > It seems that this problem occurs only for the associated Legendre > > > > > polynomials with large m; for m=0 the numerical result is accurate. > > > > > MathWorld (http://mathworld.wolfram.com/LegendrePolynomial.html) gives > > > > > a recursion relation (Eq. 66) for the associated Legendre polynomials, > > > > > and I was under the impression that this gave stable results. John, > > > > > maybe you can use this recursion relation to get better results, or > > > > > you can call the GNU Scientific Library through MathLink (http://www.gnu.org/software/gsl/). Bhuvanesh, I am very curious how you > > > > > explain this behavior. > > > > > > Roman.