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: [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.




  • Prev by Date: Re: Replace in Modules
  • Next by Date: Re: LegendreP error (bug?) in Mathematica
  • Previous by thread: Re: LegendreP error (bug?) in Mathematica
  • Next by thread: Re: LegendreP error (bug?) in Mathematica