       Re: LegendreP error (bug?) in Mathematica

• To: mathgroup at smc.vnet.net
• Subject: [mg81511] Re: LegendreP error (bug?) in Mathematica
• From: Roman <rschmied at gmail.com>
• Date: Wed, 26 Sep 2007 06:43:00 -0400 (EDT)
• References: <fcb3qh\$fk1\$1@smc.vnet.net><fd5893\$28n\$1@smc.vnet.net>

```Here is the same algorithm giving not only the final Y(L,M,x)
computed, but the entire series of intermediate results Y(L1,M,x) with
L1=|M|..L. And the same for the associated Legendre polynomials.

Mmax = 1000;

YCh = Table[(-1)^L(2L-1)!!Sqrt[(2L+1)/(4Pi(2L)!)],{L,0,Mmax}]//N;

YSC = Compile[{{L,_Integer},{M,_Integer},{x,_Real}},
Block[{A},
Which[
M>L,{},
M==L,{YCh[[L+1]](1-x^2)^(L/2)},
M==L-1,YCh[[L]](1-x^2)^((L-1)/2){1,Sqrt[2L+1]x},
True,A=YCh[[M+1]](1-x^2)^(M/2);
Prepend[Transpose[NestList[{#[]+1,#[],
Sqrt[(2#[]+1)/((#[]-M)(#[]+M))]
(Sqrt[2#[]-1]x#[]
-Sqrt[((#[]-1+M)(#[]-1-M))/
(2#[]-3)]#[])}&,
{M+2,A,Sqrt[2M+3]x A},
L-M-1]][],A]]],
{{YCh,_Real,1}}];

YS[L_Integer?NonNegative,M_Integer?NonNegative,
x_Real/;MachineNumberQ[x]&&-1<=x<=1]/;
Developer`MachineIntegerQ[L]&&
Developer`MachineIntegerQ[M]&&
M<=Mmax := YSC[L,M,x]
YS[L_Integer?NonNegative,M_Integer?Negative,
x_Real/;MachineNumberQ[x]&&-1<=x<=1]/;
Developer`MachineIntegerQ[L]&&
Developer`MachineIntegerQ[M]&&
-M<=Mmax := (-1)^M YSC[L,-M,x]
YS[L_,M_,x_] :=
Table[SphericalHarmonicY[L1,M,ArcCos[x],0],{L1,M,L}]

PS[L_,M_,x_] := Table[Sqrt[4Pi/(2L1+1)(L1+M)!/(L1-M)!],
{L1,Abs[M],L}]*YS[L,M,x]

Roman.

On Sep 23, 10:31 am, Roman <rschm... at gmail.com> wrote:
> 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,#[],
>             Sqrt[(2#[]+1)/((#[]-M)(#[]+M))]
>             (Sqrt[2#[]-1]x#[]
>             -Sqrt[((#[]-1+M)(#[]-1-M))/
>             (2#[]-3)]#[])}&,
>             {M+2,A,Sqrt[2M+3]x A},
>             L-M-1][]]],
>    {{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,#[],
> >      Sqrt[(2#[]+1)/((#[]+M)(#[]-M))](x
> > Sqrt[2#[]-1]#[]-
> >      Sqrt[(#[]+M-1)(#[]-M-1)/(2#[]-3)]#[])}&,
> >      {M+2,Y[M,M,x],Y[M+1,M,x]},L-M-1][]
>
> > 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, #[],
> > > >      (x (2 #[] + 1) #[] - (#[] + M) #[])/(#[] + 1 - M)}
> > > > &,
> > > >      {M + 1, P[M, M, x], P[M + 1, M, x]}, L - M - 1][]
>
> > > > 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[], (
> > >          z (2 #1[] + 1) #1[] - (#1[] + m) #1[])/(#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][]]]];
>
> > > 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 := LegendreP[200, 43, 4/5] // N
> > > > > > Out = 2.9256424676613492`*^97
>
> > > > > > In := LegendreP[200, 43, 0.8]
> > > > > > Out = 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:= N[10^17 + 23] - N[10^17]
>
> > > > > Out= 16.
>
> > > > > Mathematica's arbitrary precision arithmetic tells you
> > > > > that the result is not to be trusted:
>
> > > > > In:= N[10^17 + 23, 16] - N[10^17, 16]
>
> > > > > Out= 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: Any Mathematica 6 book yet?
• Next by Date: Re: how to get random numbers from a distribution
• Previous by thread: Re: LegendreP error (bug?) in Mathematica
• Next by thread: Re: PlotLegend in Mathematica 6.0.1 doesn´t work