MathGroup Archive 2007

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

Search the Archive

Re: LegendreP error (bug?) in Mathematica

On Sep 14, 3:13 am, Roman <rschm... at> 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

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.


Plot[LegendreP[200, 43, x], {x, 0, 1}]


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

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 ( 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 ( 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