Re: LegendreP error (bug?) in Mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg81362] Re: LegendreP error (bug?) in Mathematica
- From: Roman <rschmied at gmail.com>
- Date: Thu, 20 Sep 2007 03:55:32 -0400 (EDT)
- References: <fcb3qh$fk1$1@smc.vnet.net><fcnle7$sac$1@smc.vnet.net>
Bobby,
The only time I have used them was in quantum mechanics, concerning
systems of spherical symmetry. There you only need the spherical
harmonics though, not the associated Legendre polynomials, and so
there are no normalization issues. Very large values of the arguments
come in when you try to go to the classical limit of quantum
mechanics, where huge sums over huge basis sets are needed to
represent quantities accurately. Think of a Rydberg atom or molecule,
where an electron is on a quasi-classical trajectory of very large
angular momentum (in units of PlanckConstantReduced) around an ion. In
order to describe this system quantum mechanically, you'll need lots
of spherical harmonics with parameters on the order of this
dimensionless angular momentum. It is then very instructive to see how
Kepler's laws follow from quantum mechanics in this limit of large
angular momentum.
On the other hand, semi-classical techniques have been developed back
in pre-computer days when people wanted to do the above calculation
without manually evaluating spherical harmonics with large parameters.
In that sense, I don't know where you would *really* need them.
Roman.
On Sep 19, 11:49 am, DrMajorBob <drmajor... at bigfoot.com> wrote:
> Just curious, mind you...
>
> What are these huge values of LegendreP likely to be used for?
>
> Bobby
>
>
>
> On Tue, 18 Sep 2007 04:54:01 -0500, 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]]+1,#[[3]],
> > Sqrt[(2#[[1]]+1)/((#[[1]]+M)(#[[1]]-M))](x
> > Sqrt[2#[[1]]-1]#[[3]]-
> > Sqrt[(#[[1]]+M-1)(#[[1]]-M-1)/(2#[[1]]-3)]#[[2]])}&,
> > {M+2,Y[M,M,x],Y[M+1,M,x]},L-M-1][[3]]
>
> > 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]] + 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.
>
> --
>
> DrMajor... at bigfoot.com