       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 := LegendreP[200, 43, 4/5] // N
> Out = 2.9256424676613492`*^97
>
> In := LegendreP[200, 43, 0.8]
> Out = 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:= N[10^17 + 23] - N[10^17]

Out= 16.

Mathematica's arbitrary precision arithmetic tells you
that the result is not to be trusted:

In:= N[10^17 + 23, 16] - N[10^17, 16]

Out= 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