MathGroup Archive 2007

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

Search the Archive

Re: LegendreP error (bug?) in Mathematica

  • To: mathgroup at smc.vnet.net
  • Subject: [mg81287] Re: LegendreP error (bug?) in Mathematica
  • From: Roman <rschmied at gmail.com>
  • Date: Tue, 18 Sep 2007 00:47:39 -0400 (EDT)
  • References: <fcdf41$pjg$1@smc.vnet.net><fclajn$f4a$1@smc.vnet.net>

John,
I've assembled a MathLink program to compute spherical harmonics with
high L and M by calling the GNU Scientific Library. You can use this
to compute the associated Legendre polynomials by using Eq. (6) of
http://mathworld.wolfram.com/SphericalHarmonic.html

Just write a file with the following contents:

----- start of sph.tm
--------------------------------------------------------
#include <gsl/gsl_sf_legendre.h>

:Begin:
:Function:      gsl_sf_legendre_sphPlm
:Pattern:       YLM[L_?IntegerQ,M_?IntegerQ,x_?NumericQ]/;((L>=0)&&(-
L<=M<=L)&&(-1<=x<=1))
:Arguments:     {L,M,x}
:ArgumentTypes: {Integer,Integer,Real}
:ReturnType:    Real
:End:

:Evaluate:      YLM::usage = "YLM[L,M,x] computes
SphericalHarmonicY[L,M,ArcCos[x],0] by calling the GNU Scientific
Library. Multiply by E^(I*M*phi) to get
SphericalHarmonicY[L,M,ArcCos[x],phi]."

int main(int argc, char *argv[]) {
  return MLMain(argc,argv);
}
----- end of sph.tm
--------------------------------------------------------

You obviously need to have the GSL installed to compile this. Here's
how to compile:

- on a Macintosh, run these two commands in this order:
/Applications/Mathematica.app/SystemFiles/Links/MathLink/DeveloperKit/
CompilerAdditions/mprep sph.tm -o sphtm.c
gcc -Wall -O2 -I/Applications/Mathematica.app/SystemFiles/Links/
MathLink/DeveloperKit/CompilerAdditions -lgsl -L/Applications/
Mathematica.app/SystemFiles/Links/MathLink/DeveloperKit/
CompilerAdditions -lMLi3 sphtm.c -o sph

- on 64-bit Linux, run these:
/usr/local/Wolfram/Mathematica/6.0/SystemFiles/Links/MathLink/
DeveloperKit/Linux-x86-64/CompilerAdditions/mprep sph.tm -o sphtm.c
gcc -Wall -O2 -I/usr/local/Wolfram/Mathematica/6.0/SystemFiles/Links/
MathLink/DeveloperKit/Linux-x86-64/CompilerAdditions -lgsl -lgslcblas -
L/usr/local/Wolfram/Mathematica/6.0/SystemFiles/Links/MathLink/
DeveloperKit/Linux-x86-64/CompilerAdditions -lML64i3 -lm -lpthread -
lrt sphtm.c -o sph

On other types of machines, you'll have to figure out the exact
library locations and necessary libraries somehow. There are examples
in the directory /Applications/Mathematica.app/SystemFiles/Links/
MathLink/DeveloperKit/MathLinkExamples (or similar on a different kind
of machine), and the Makefile contains some hints on how to compile.

Once that's done, you have an executable named "sph" which you can
call from Mathematica:

Install["~/the/pertinent/directory/sph"]
?YLM
SphericalHarmonicY[200, 43, ArcCos[0.8], 0]
SphericalHarmonicY[200, 43, ArcCos[4/5], 0]//N
YLM[200, 43, 0.8]

Cheers!
Roman.

On Sep 17, 9: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




  • Prev by Date: Re: Defining a Function with an Indeterminate Number of Arguments
  • Next by Date: Re: coupled map lattice problem
  • Previous by thread: Re: LegendreP error (bug?) in Mathematica
  • Next by thread: Re: LegendreP error (bug?) in Mathematica