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.