Re: Re: LegendreP Evaluation Mystery
- To: mathgroup at smc.vnet.net
- Subject: [mg75194] Re: [mg75117] Re: LegendreP Evaluation Mystery
- From: Carl Woll <carlw at wolfram.com>
- Date: Fri, 20 Apr 2007 00:35:57 -0400 (EDT)
- References: <evvaqs$95r$1@smc.vnet.net><f013r8$8ih$1@smc.vnet.net> <200704180852.EAA07964@smc.vnet.net>
Antonio wrote:
>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.
>
>
If you use Set, you are *not* using Mathematica's algorithm for
computing LegendreP functions. Notice:
In[11]:=P1[x_] = LegendreP[4, 3, x]
Out[11]= -105 x (1 - x^2)^(3/2)
does not have the head LegendreP. So, when evaluating P1[.7], you are
evaluating the expression -105 .7 (1 - .7^2)^(3/2) and not the
expression LegendreP[4,3,.7]. Mathematica is smart enough to figure out
a good machine precision approximation to LegendreP[4,3,.7]. On the
other hand, Mathematica knows nothing about the expression -105 x
(1-x^2)^(3/2), so all it can do is plug in .7 and use standard machine
precision arithmetic. For n=4, m=3, plugging in .7 and using arithmetic
produces a reasonable answer. For large n, small m, standard machine
precision arithmetic will typically produce a completely incorrect
answer due to numeric instability.
Perhaps you can give a simple example of the integrals you are trying to
evaluate. Some possible remedies might include using extended precision
arithmetic with the polynomial form of the LegendreP functions, or using
FunctionInterpolation to create a substitute for the LegendreP functions.
Carl Woll
Wolfram Research
>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
>
>
>
>
>
>
>
- References:
- Re: LegendreP Evaluation Mystery
- From: Antonio <aneves@gmail.com>
- Re: LegendreP Evaluation Mystery