Re: LegendreP error (bug?) in Mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg81188] Re: LegendreP error (bug?) in Mathematica
- From: sashap <pavlyk at gmail.com>
- Date: Sat, 15 Sep 2007 04:07:39 -0400 (EDT)
- References: <fcb3qh$fk1$1@smc.vnet.net><fcdfqu$pub$1@smc.vnet.net>
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.