Re: LegendreP Evaluation Mystery
- To: mathgroup at smc.vnet.net
- Subject: [mg75247] Re: LegendreP Evaluation Mystery
- From: Roman <rschmied at gmail.com>
- Date: Sat, 21 Apr 2007 23:17:58 -0400 (EDT)
- References: <evvaqs$95r$1@smc.vnet.net><f013r8$8ih$1@smc.vnet.net>
Antonio:
> Well if I try to evaluate LegendreP numerically by recursion relation, no=
t only will it be slow but also carry on any round offs up to the high n's =
I am interested in.
I think you misunderstand here: all numbers will always carry a
roundoff on the order of the machine precision, but this is not a
problem unless this error exponentially (or otherwise badly) increases
during each step of a recursion relation. If you have a stable
recursion, this is not a problem, and is indeed how things are done
properly.
> Since I need all {n,m}'s up to the big ones none, I think it is best to h=
ave the extended function (polynimial) determined by Rodrigues.
Again, not really. Evaluating high-order polynomials numerically is
not a simple issue, as you can see in Numerical Recipes, section on
"evaluation of functions" - "polynomials and rational functions".
Abramowitz and Stegun list some useful recursion relations, and I've
assembled a procedure that will generate all the LegendreP[n,m,x] for
all {n,m} up to a given value:
first, define a list of coefficients to be used later, to speed things
up a little:
coefflist[n_] := coefflist[n] = N[Table[{1/a-1, 2-1/a}, {a,1,n}]]
make a list of LegendreP[a,z] for a=0..n by recursion of the formula
P[n,z] = (2-1/n)*z*P[n-1,z]-(1-1/n)*P[n-2,z]
P0list[n_, z_] := Transpose[FoldList[{#1[[2]], ({1, z}#1).#2} &,
{0,1}, coefflist[n]]][[2]]
lastly, from these P[n,z] compute the associated Legendre polynomials
from the formula
P[n,m+1,z] = ((n-m)*z*P[n,m,z]-(n+m)*P[n-1,m,z])/sqrt(1-z^2)
Plist[n_, z_] := Module[{zz,nextm},
zz = 1/Sqrt[1-z^2];
nextm[m_, p_] := MapIndexed[zz(#2[[1]] z #1[[2]]-(2m
+#2[[1]])#1[[1]])&,
Transpose[{Most[p],Rest[p]}]];
FoldList[nextm[#2,#1]&, P0list[n,z], Range[0,n-1]]]
Now if you call Plist[10,0.3] it will calculate all the P[n,m,0.3] for
n=0..10 and m=0..n
This should be fast and numerically stable. Maybe you can make it even
faster by using Compile[]; please ask if you need help with that.
Roman.
On Apr 20, 6:52 am, "Antonio Neves" <ane... at gmail.com> wrote:
> 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 <c... 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 illustra=
te
> > >>>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 arithm=
etic
> > 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 function=
s=2E
>
> > 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] -
>
> ...
>
> read more =BB