Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: LegendreP Evaluation Mystery

• To: mathgroup at smc.vnet.net
• Subject: [mg75117] Re: LegendreP Evaluation Mystery
• From: Antonio <aneves at gmail.com>
• Date: Wed, 18 Apr 2007 04:52:38 -0400 (EDT)
• References: <evvaqs$95r$1@smc.vnet.net><f013r8$8ih$1@smc.vnet.net>

On 17 abr, 01:19, "Norbert Marxer" <mar... at mec.li> wrote:
> On 16 Apr., 10:06, "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
>
> Hello
>
> If you use SetDelayed (:=) instead of Set (=) in your defintions for
> m, x, P and J then all your mysteries and questions will (probably)
> evaporate.
>
> With your definition the random numbers are evaluated only once and
> therefore m is constant.
>
> Similarly your associated Legendre polynomial (P) is evaluated only
> once (for specific values of n and m) and set to a specific constant
> polynomial with the independent variable y. The evaluation of this
> polynomial takes obviously less time than the calculation and
> evaluation of the associated Legendre polynomial.
>
> Best Regards
> Norbert Marxer- Ocultar texto entre aspas -
>
> - Mostrar texto entre aspas -

Dear Norbert and list members,

The Set instead of SetDelayed, was purposely placed so it will be
evaluated only once.
The problem with the example above is that the extended function
numerical result isn't identical with the one for LegendreP[n,m,x]
function. So does Mathematica evaluate LegendreP with a different
algorithm for low and high n's?

Any help with this question is appreciated.

Anyhow continuing to determine accuratelly an extended form for
associated LegendreP of high n's, I wrote:

P1 = Table[Sqrt[((2*n + 1)*(n - m)!)/((4*Pi)*(n + m)!)]*LegendreP[n,
m, x], {n, 1, 52}, {m, 0, n}];
P2 = Table[(-1)^m*((Sqrt[((2*n + 1)*(n - m)!)/((4*Pi)*(n + m)!)]*(1 -
x^2)^(m/2)*D[(x^2 - 1)^n, {x, n + m}])/(n!*2^n)), {n, 1, 52},
    {m, 0, n}];
A[n_, m_, y_] := P1[[n,m + 1]] /. x -> y;
B[n_, m_, y_] := P2[[n,m + 1]] /. x -> y;

This defines P1 using Mathematicas algoritm for LegendreP, and P2
using Rodrigues.
When plotting, for a particular n (big),m (Note that I have included a
normalization factor to avoid big numbers, Sqrt[...]):

n = 49; m = 12;
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x],
{x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
   PlotLabel -> "LegendreP"];
Plot[A[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
PlotLabel -> "P1"];
Plot[B[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
PlotLabel -> "P2"];
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x] -
A[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
   PlotLabel -> "LegendreP-P1"];
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x] -
B[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
   PlotLabel -> "LegendreP-P2"];
Plot[A[n, m, x] - B[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -
> 100, PlotLabel -> "P1-P2"];

>From this output we can see that there is a problem for Mathematicas
definition of LegendreP[n,m,x] for x around 0.6.  While the output of
P1 and P2 seems identical. If we change n and m again.

n = 51; m = 3;
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x],
{x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
   PlotLabel -> "LegendreP"];
Plot[A[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
PlotLabel -> "P1"];
Plot[B[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
PlotLabel -> "P2"];
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x] -
A[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
   PlotLabel -> "LegendreP-P1"];
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x] -
B[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 100,
   PlotLabel -> "LegendreP-P2"];
Plot[A[n, m, x] - B[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints -
> 100, PlotLabel -> "P1-P2"];

We see that the definition P1 (from Mathematica's LegendreP) shows
some numerical instabilities for x near 1. While for P2 from Rodrigues
this doesn't happen. The problem with LegendreP[n,m,x] still occurs in
this case for x around 0.85, but smaller than in the previous example.
Is Rodrigues the best way to build an array (n,m) to represent the
associated Legendre function?

Thanks,
Antonio

• Follow-Ups:
  • Re: Re: LegendreP Evaluation Mystery
    • From: Carl Woll <carlw@wolfram.com>