MathGroup Archive 2007

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

Search the Archive

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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg81229] Re: [mg81203] Re: LegendreP error (bug?) in Mathematica
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Mon, 17 Sep 2007 03:29:12 -0400 (EDT)
  • References: <fcdf41$pjg$1@smc.vnet.net> <200709150815.EAA28245@smc.vnet.net>

If an expression contains any MachinePrecision number  Mathematica  
does not track precision at all and there is no guarantee that the  
outcome of your computations is not total nonsense. In such cases any  
precision tracking has to be done "by hand".
There is completely no point asking for th Precision of such a  
computation because such Mathematica's notion of precision in such  
situations is purely formal - they always have MachinePrecision.
If you want Mathematica to really try to track precision you need to  
use arbitrary precision input. This is the basic Mathematica design  
and there is no point complaining about it: it certainly is not going  
to be changed.

By the way: I am not aware of any system other than Mathematica that  
offers any automatic precision tracking; certianly Mathematica's best  
known rival does not do it.

Andrzej Kozlowski


On 15 Sep 2007, at 17:15, Roman 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
>
>
>



  • Prev by Date: Button Help Example
  • Next by Date: Re: Help: Font problem
  • Previous by thread: Re: LegendreP error (bug?) in Mathematica
  • Next by thread: Re: LegendreP error (bug?) in Mathematica