Re: LegendreP error (bug?) in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg81275] Re: LegendreP error (bug?) in Mathematica*From*: Yaroslav Bulatov <yaroslavvb at gmail.com>*Date*: Tue, 18 Sep 2007 00:41:25 -0400 (EDT)*References*: <fcdf41$pjg$1@smc.vnet.net><fclajn$f4a$1@smc.vnet.net>

I agree with Roman that Mathematica's implementation of LegendreP is suboptimal. Consider the precision of native LegendreP[n,x] implementation vs 3-term recursion formula as a function of n legendre3term[0, x] = 1; legendre3term[1, x] = x; legendre3term[n_, x_] := Module[{p1 = 1, p2 = x}, For[i = 2, i < n + 1, i++, {p1, p2} = {p2, ((2*i - 1)*x*p2 - (i - 1)*p1)/i}; ]; p2 ]; ints = Block[{$MinPrecision = 50, $MaxPrecision = 50}, x = N[Interval[1/10], $MinPrecision]; {legendre3term[#, x], LegendreP[#, x]} & /@ Range[1, 100] ] /. {Interval[{a_, b_}] :> b - a}; ListLinePlot[Transpose[-Log[10, ints]], PlotLabel -> "Precision of LegendreP and 3-term recursion for large n", AxesLabel -> {"n", "Precision"}] in this example, native LegendreP implementation loses .55 digits of precision each time n increases by 1, while 3-term recursion loses 0.041 legendre3term also runs faster than LegendreP for a few values I tried Yaroslav On Sep 17, 12:33 am, Roman <rschm... at gmail.com> wrote: > John, > If you only need the spherical harmonics, but not the associated > Legendre polynomials, then it turns out that the GNU Scietific Library > does very well: > http://www.gnu.org/software/gsl/manual/html_node/Associated-Legendre-... > I've tried this with very large numbers L and M, and the results were > very accurate (unlike the Mathematica results). > You can write a small C wrapper program that you then call using > MathLink. If you're not familiar with this, just ask. > Roman. > > On Sep 15, 10:24 am, Roman <rschm... at gmail.com> wrote: > > > 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

**Re: LegendreP error (bug?) in Mathematica**

**Re: Single-step evaluation in Mathematica**

**Re: LegendreP error (bug?) in Mathematica**

**Re: LegendreP error (bug?) in Mathematica**