MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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.






  • Prev by Date: Transpose[ InterpolatingFunction, { n1, n2, ...} ]
  • Next by Date: Re: Re: Problem with inverse laplace transform (FIX)
  • Previous by thread: Re: Re: LegendreP error (bug?) in Mathematica
  • Next by thread: Re: LegendreP error (bug?) in Mathematica