MathGroup Archive 2007

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

Search the Archive

Re: LegendreP Evaluation Mystery

  • To: mathgroup at smc.vnet.net
  • Subject: [mg75134] Re: LegendreP Evaluation Mystery
  • From: Roman <rschmied at gmail.com>
  • Date: Wed, 18 Apr 2007 05:01:21 -0400 (EDT)
  • References: <evvaqs$95r$1@smc.vnet.net>

Antonio,

The evaluation of Legende polynomials is a tricky business. For high
indices n and m, the explicit writing as a polynomial, as you do in
P[y], is in fact numerically unstable as you have noticed. If you look
at the resulting polynomial P[y], it contains very large terms, but we
know that the result must be between -1 and +1, so there is a very
delicate cancellation of huge terms going on. Numerically this is just
asking for trouble.

The function LegendreP[...] does not simply evaluate a polynomial but
rather uses a numerically stable recursion relation, which then turns
out to be slower but much better behaved. See, e.g.,
   http://en.wikipedia.org/wiki/Legendre_polynomials
at the section "Additional properties of Legendre polynomials"

The same is probably true for Bessel functions. Maybe you can
implement your own Legendre and Bessel procedures starting from the
"Numerical Recipes" book, chapter "Special Functions". Or you could
call the GNU Scientific Library through MathLink
   http://www.gnu.org/software/gsl/
But I doubt you'll do better than the Mathematica built-in numerical
routines, except if you need a whole series of LegendreP[n,m,x] for
fixed x but varying n,m (the recursion relations tend to produce not
only the LegendreP[n,m,x] you are asking for, but all those with
smaller indices as well).

Or you can just define P[y] etc. through SetDelayed (:=), as Norbert
Marxer points out. Or do some steps of your calculation analytically
instead of numerically, if possible.

Cheers!
Roman


On Apr 16, 10:06 am, "Antonio" <ane... at gmail.com> wrote:
> Dear all,
>
> I have found that Mathematica v5.2 evaluates BesselJ and LegendreP
> differently. Why is this? I have written this test below to illustrate
> the timing.
>
> In[11]:=
> k=10000;
> n=51;
> m=Random[Integer,{0,n}];
> x=Random[Real,{0,1}];
> P[y_]=LegendreP[n,m,y];
> J[y_]=BesselJ[n,y];
> Timing[Do[LegendreP[n,m,x],{k}]]
> Timing[Do[P[x],{k}]]
> Timing[Do[BesselJ[n,x],{k}]]
> Timing[Do[J[x],{k}]]
>
> Out[17]=
> {4.531 Second,Null}
>
> Out[18]=
> {0.484 Second,Null}
>
> Out[19]=
> {0.61 Second,Null}
>
> Out[20]=
> {0.657 Second,Null}
>
> The emphasis here is with the Associated Legendre Function, since it is
> a bottel neck of my current calculations for big n's (n>50), used
> inside NIntegrate. The form LegendreP[n,m,x] takes longer to evaluate
> than in extended form (P[y])? If I try to use the extended form to do
> some numerical integration, it results in a wrong result even though
> it is faster.
>
> When plotting the functions, it seems to show some numerical
> instabilities (for low m):
>
> n = 51; m = 3;
> Plot[LegendreP[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints ->
> 100];
> Plot[P[x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100];
>
> And it is worse for m=0. Does Mathematica evaluate LegendreP
> differently for high n's, why does it take so long? Is there any way
> that I build an array of extended Associated Legendre functions, so as
> to speed up calculations and wouldn't fail numerically as above?
>
> Antonio




  • Prev by Date: Re: question about Protect
  • Next by Date: Re: how make function of solution by NDSolve depending
  • Previous by thread: Re: LegendreP Evaluation Mystery
  • Next by thread: Re: LegendreP Evaluation Mystery