Re: LegendreP error (bug?) in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg81271] Re: LegendreP error (bug?) in Mathematica*From*: sashap <pavlyk at gmail.com>*Date*: Tue, 18 Sep 2007 00:39:21 -0400 (EDT)*References*: <fcb3qh$fk1$1@smc.vnet.net><fclahp$f37$1@smc.vnet.net>

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.