Re: LegendreP error (bug?) in Mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg81231] Re: LegendreP error (bug?) in Mathematica
- From: Roman <rschmied at gmail.com>
- Date: Mon, 17 Sep 2007 03:30:16 -0400 (EDT)
- References: <fcb3qh$fk1$1@smc.vnet.net><fcg3uu$rdh$1@smc.vnet.net>
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) of http://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. 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.