Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2007
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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.




  • Prev by Date: Re: Changing Textbook.nb to label equations
  • Next by Date: Re: Automated documentation (like javadoc) for packages?
  • Previous by thread: Re: Re: LegendreP error (bug?) in Mathematica
  • Next by thread: Re: LegendreP error (bug?) in Mathematica