Re: Re: LegendreP error (bug?) in Mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg81333] Re: [mg81290] Re: LegendreP error (bug?) in Mathematica
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Wed, 19 Sep 2007 05:36:15 -0400 (EDT)
- References: <fcb3qh$fk1$1@smc.vnet.net> <fcnle7$sac$1@smc.vnet.net> <21737589.1190175665164.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
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 <rschmied 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.
>
>
>
>
--
DrMajorBob at bigfoot.com