Re: LegendreP error (bug?) in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg81203] Re: LegendreP error (bug?) in Mathematica*From*: Roman <rschmied at gmail.com>*Date*: Sat, 15 Sep 2007 04:15:29 -0400 (EDT)*References*: <fcdf41$pjg$1@smc.vnet.net>

Bob, Of course it is a precision problem. The question is, why does Mathematica's numerical algorithm for associated Legendre polynomials give wrong results for large m? The inaccuracies are far beyond what you'd expect from the limited precision of the input arguments. To be more specific, look at my previous example and series-expand around x=0.8: In[1] := Series[LegendreP[200, 43, x], {x, 4/5, 1}] // Normal // N Out[1] = 2.92564*10^97 + 1.37446*10^100 (-0.8 + x) >From the ratio of the coefficients you see that, roughly speaking, if I insert x=0.8 with a relative error of 10^(-16), then I expect a result with a relative error of 4*10^(-14). In general, I expect the relative error to be around 400 times greater in the result than in the argument. But not so for Mathematica's numerical algorithm. If you increase the relative precision, as you say, In[2] := LegendreP[200, 43, 0.8`50] Out[2] = 2.925642467661265646732164377813044273198`24.332438293726835*^97 you get a result which is ok but more than 25 orders of magnitude less precise than the argument, far worse than the factor of 400 expected from the slope of the Legendre polynomial! So if you start with a machine number with 16 significant figures (MachinePrecision), and lose 25 in the calculation, then naturally you end up with garbage. But, instead of telling you that you are getting garbage, the algorithm returns a number with exaggerated precision: In[3] := LegendreP[200, 43, 0.8] // Precision Out[3] = MachinePrecision Thus I agree with John's complaint that Mathematica's numerical algorithm is not good. At least you should be getting a warning of some sort if the results are so bad. Unfortunately, it's no better if you use the GNU Scientific Library through MathLink, since the GSL manual says: "The following functions compute the associated Legendre Polynomials P_l^m(x). Note that this function grows combinatorially with l and can overflow for l larger than about 150. There is no trouble for small m, but overflow occurs when m and l are both large. Rather than allow overflows, these functions refuse to calculate P_l^m(x) and return GSL_EOVRFLW when they can sense that l and m are too big." Roman. On Sep 14, 10:01 am, Bob Hanlon <hanl... at cox.net> wrote: > In version 6 I do not see a problem with the first two examples. In the third example, l is undefined. Defining l and changing 3.7 to 37/10 (use rational numbers to maintain high precision) works fine. Alternatively, specify higher precision (e.g., 3.7`25). Also, in version 6, Plot has an option to change (increase) the WorkingPrecision. > > Bob Hanlon > > ---- John Ralston <rals... at ku.edu> wrote: > > > LegendreP[ l, m] and SphericalHarmonicY[ t, p, l, m] go > > wrong for large index l . > > > For l> 40 or so neither can be used reliably everywhere. > > > To see the breakdown plot the functions. Not all index cases fail. > > Here's some examples: > > > Plot[LegendreP[ 44, 13, x] , {x, -1, 1} ] > > > Plot[LegendreP[ 66, 9, x] , {x, -1, 1} ] > > > ListPlot[ Table[ {j, > > Sqrt[ 4Pi/ ( > > 2l + 1) ] Abs[ SphericalHarmonicY[ j, 0, 0, > > 3.7 ]]}, {j, 1, 55}], PlotJoined -> True] > ; > > Has anyone fixed this? Does anyone care? > > I need l ranges above 200. > > > thanks > > John Ralston > > > > I find serious bugs in Mathematica 5.1 LegendreP > > > and SphericalHarmonicY. It is not a matter of > > > definitions or syntax, but a catastrophic failure > > > easy to establish. Math archives show a history > > > of discussion, but invariably centered on variations > > > of definition. > > > > Wolfram Research shows little interest. I'm > > > wondering if the failure > > > to perfom is well known and has been repaired by > > > someone > > > re-writing the commands. I don't often follow this > > > forum, but I've joined to either get access to code > > > that works, or to inform people so they can go after > > > the problem. Does anyone care? > > > > John Ralston

**Follow-Ups**:**Re: Re: LegendreP error (bug?) in Mathematica***From:*Andrzej Kozlowski <akoz@mimuw.edu.pl>