Re: Re: LegendreP Evaluation Mystery
- To: mathgroup at smc.vnet.net
- Subject: [mg75205] Re: [mg75117] Re: LegendreP Evaluation Mystery
- From: "Antonio Neves" <aneves at gmail.com>
- Date: Fri, 20 Apr 2007 00:41:35 -0400 (EDT)
- References: <evvaqs$95r$1@smc.vnet.net> <f013r8$8ih$1@smc.vnet.net>
Dear Carl and list members,
Thank you for your suggestions. I believe that maybe the extended
precision arithmatic would be too time consuming, but the posibility
of using Function Interpolation seems interesting.
It is now clear to me that with the "Set" instead of Set Delayed the
polynomial and not the LegendreP is passed on. But still, when doing a
plot, why is there a problem with tha evaluation of LegendreP? As in,
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 -> 200,
PlotLabel -> "LegendreP[n,m,x]"];
I suppose that LegendreP[n,m,x] will behave the same way inside Plot
as in NIntegrate, true? Using "Set", as you said in your post, I
obtain a smooth curve, without apart numerical singularity around
x=0.6.
n = 49; m = 12;
NewFunc[n_, m_,
x_] = Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x];
Plot[NewFunc[n, m, x], {x, 0, 1}, PlotRange -> All,
PlotPoints -> 200, PlotLabel -> "NewFunc[n,m,x]"];
This I do not know how to explain and it seems to be a problem of how
LegendreP[n,m,x] is evaluated. But lets say, okay, lets have a lower
m, such as m=1. Then,
n = 49; m = 1;
NewFunc[n_, m_,
x_] = Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m, x];
Plot[NewFunc[n, m,
x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 200, PlotLabel -> "
NewFunc[n,m,x]"];
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[n, m,
x], {x, 0, 1}, PlotRange -> All, PlotPoints -> 200, PlotLabel -> "
LegendreP[n,m,x]"];
Yes we see some numerical instabilities, probably due to standard machine
precision, when we get nearer to x=1 (due to the (1-x)^(1/2) terms in
the polinomial). But then again, if we use,
Pnm = 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, 51}, {m, 0, n}];
PnmO = Map[Experimental`OptimizeExpression, Pnm, {2}];
PnmC = Map[(Compile[{{x, _Real}}, #]) &, PnmO, {2}];
P[n_, m_, x_?NumericQ] := PnmC[[n, m + 1]][x];
n = 49; m = 1;
Plot[Sqrt[((2*n + 1)/(4*Pi))*((n - m)!/(n + m)!)]*LegendreP[
n, m, x], {x, 0, 1}, PlotRange -> All,
PlotPoints -> 200, PlotLabel -> "LegendreP[n,m,x]"];
Plot[P[n, m, x], {x, 0, 1}, PlotRange -> All, PlotPoints ->
200, PlotLabel -> "P[n,m,x]"];
We obtain P[n,m,x] as a curve without the numerical instabilities as
the NewFunc[n,m,x] above, and without the singular type behaviour as
LegendreP[n,m,x]. I am using Mathematica 5.2.0.0 on Windows XP 32-bit
system, I do not know if the singularities/numerical instabilities
appears for 64bit system.
And finally, my problematic integral that resulted in this query.
(* Some typical values *)
nGlass = 1.55;
nWater = 1.33;
Crit = ArcSin[nWater/nGlass];
ko = 7.98;
zo = 0.03;
d = 50.0;
fw = 0.68;
(*The Integral*)
IntFunc[(n_Integer)?Positive, (m_Integer)?NonNegative, (r_Real)?Positive] :=
Block[{func}, func[(a_)?NumericQ] := I^(n - m)*Cos[a]*Sqrt[((2*n +
1)/(n*(n + 1)))*((n - m)!/(n + m)!)]*Sqrt[Cos[a]]*Exp[-(fw*Sin[a])^2]*
Exp[(-I)*ko*nWater*zo*Sqrt[1 -
(nGlass*(Sin[a]/nWater))^2]]*Exp[(-I)*ko*d*(nGlass*Cos[a] -
nWater*Sqrt[1 - (nGlass*(Sin[a]/nWater))^2])]*
(m^2*LegendreP[n, m, Sqrt[1 -
(nGlass*(Sin[a]/nWater))^2]]*(BesselJ[m,
ko*nGlass*r*Sin[a]]/(ko*nGlass*r*Sin[a]*
(nGlass*Cos[a] + nWater*Sqrt[1 - (nGlass*(Sin[a]/nWater))^2]))) -
((1/2)*(BesselJ[-1 + m, ko*nGlass*r*Sin[a]] - BesselJ[1 + m,
ko*nGlass*r*Sin[a]]))*
(((-1 - n)*Sqrt[1 - (nGlass*(Sin[a]/nWater))^2]*LegendreP[n,
m, Sqrt[1 - (nGlass*(Sin[a]/nWater))^2]] +
(1 - m + n)*LegendreP[1 + n, m, Sqrt[1 -
(nGlass*(Sin[a]/nWater))^2]])/(nWater*Cos[a] + nGlass*Sqrt[1 -
(nGlass*(Sin[a]/nWater))^2])));
NIntegrate[func[a], {a, 0, Crit}, AccuracyGoal -> Infinity,
PrecisionGoal -> 7, MaxRecursion -> 30, MinRecursion -> 6, Method ->
DoubleExponential,
Compiled -> False, SingularityDepth -> 1000]];
Trying,
n = 50; m = 4; r = 3.21;
IntFunc[n, m, r]
Takes too long to evaluate, I believe is LegendreP[n,m,x] is the
bottleneck, and I need to determine for each element of a grid as in,
Grid = Table[{n, m, r}, {0, 70}, {m, 0, n}, {r, 0.001, 5, .1}];
Carl, how do you suggest to use Function Interpolation?
Regards,
Antonio
On 4/19/07, Carl Woll <carlw at wolfram.com> wrote:
> 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
> >
> >
> >
> >
> >
> >
> >
>
>