Re: LegendreP error (bug?) in Mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg81203] Re: LegendreP error (bug?) in Mathematica
- From: Roman <rschmied at gmail.com>
- Date: Sat, 15 Sep 2007 04:15:29 -0400 (EDT)
- References: <fcdf41$pjg$1@smc.vnet.net>
Bob,
Of course it is a precision problem. The question is, why does
Mathematica's numerical algorithm for associated Legendre polynomials
give wrong results for large m? The inaccuracies are far beyond what
you'd expect from the limited precision of the input arguments. To be
more specific, look at my previous example and series-expand around
x=0.8:
In[1] := Series[LegendreP[200, 43, x], {x, 4/5, 1}] // Normal // N
Out[1] = 2.92564*10^97 + 1.37446*10^100 (-0.8 + x)
>From the ratio of the coefficients you see that, roughly speaking, if
I insert x=0.8 with a relative error of 10^(-16), then I expect a
result with a relative error of 4*10^(-14). In general, I expect the
relative error to be around 400 times greater in the result than in
the argument. But not so for Mathematica's numerical algorithm. If you
increase the relative precision, as you say,
In[2] := LegendreP[200, 43, 0.8`50]
Out[2] =
2.925642467661265646732164377813044273198`24.332438293726835*^97
you get a result which is ok but more than 25 orders of magnitude less
precise than the argument, far worse than the factor of 400 expected
from the slope of the Legendre polynomial! So if you start with a
machine number with 16 significant figures (MachinePrecision), and
lose 25 in the calculation, then naturally you end up with garbage.
But, instead of telling you that you are getting garbage, the
algorithm returns a number with exaggerated precision:
In[3] := LegendreP[200, 43, 0.8] // Precision
Out[3] = MachinePrecision
Thus I agree with John's complaint that Mathematica's numerical
algorithm is not good. At least you should be getting a warning of
some sort if the results are so bad.
Unfortunately, it's no better if you use the GNU Scientific Library
through MathLink, since the GSL manual says: "The following functions
compute the associated Legendre Polynomials P_l^m(x). Note that this
function grows combinatorially with l and can overflow for l larger
than about 150. There is no trouble for small m, but overflow occurs
when m and l are both large. Rather than allow overflows, these
functions refuse to calculate P_l^m(x) and return GSL_EOVRFLW when
they can sense that l and m are too big."
Roman.
On Sep 14, 10:01 am, Bob Hanlon <hanl... at cox.net> wrote:
> In version 6 I do not see a problem with the first two examples. In the third example, l is undefined. Defining l and changing 3.7 to 37/10 (use rational numbers to maintain high precision) works fine. Alternatively, specify higher precision (e.g., 3.7`25). Also, in version 6, Plot has an option to change (increase) the WorkingPrecision.
>
> Bob Hanlon
>
> ---- John Ralston <rals... at ku.edu> wrote:
>
> > LegendreP[ l, m] and SphericalHarmonicY[ t, p, l, m] go
> > wrong for large index l .
>
> > For l> 40 or so neither can be used reliably everywhere.
>
> > To see the breakdown plot the functions. Not all index cases fail.
> > Here's some examples:
>
> > Plot[LegendreP[ 44, 13, x] , {x, -1, 1} ]
>
> > Plot[LegendreP[ 66, 9, x] , {x, -1, 1} ]
>
> > ListPlot[ Table[ {j,
> > Sqrt[ 4Pi/ (
> > 2l + 1) ] Abs[ SphericalHarmonicY[ j, 0, 0,
> > 3.7 ]]}, {j, 1, 55}], PlotJoined -> True]
> ;
> > Has anyone fixed this? Does anyone care?
> > I need l ranges above 200.
>
> > thanks
> > John Ralston
>
> > > I find serious bugs in Mathematica 5.1 LegendreP
> > > and SphericalHarmonicY. It is not a matter of
> > > definitions or syntax, but a catastrophic failure
> > > easy to establish. Math archives show a history
> > > of discussion, but invariably centered on variations
> > > of definition.
>
> > > Wolfram Research shows little interest. I'm
> > > wondering if the failure
> > > to perfom is well known and has been repaired by
> > > someone
> > > re-writing the commands. I don't often follow this
> > > forum, but I've joined to either get access to code
> > > that works, or to inform people so they can go after
> > > the problem. Does anyone care?
>
> > > John Ralston
- Follow-Ups:
- Re: Re: LegendreP error (bug?) in Mathematica
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Re: Re: LegendreP error (bug?) in Mathematica