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