Re: LegendreP Evaluation Mystery

*To*: mathgroup at smc.vnet.net*Subject*: [mg75134] Re: LegendreP Evaluation Mystery*From*: Roman <rschmied at gmail.com>*Date*: Wed, 18 Apr 2007 05:01:21 -0400 (EDT)*References*: <evvaqs$95r$1@smc.vnet.net>

Antonio, The evaluation of Legende polynomials is a tricky business. For high indices n and m, the explicit writing as a polynomial, as you do in P[y], is in fact numerically unstable as you have noticed. If you look at the resulting polynomial P[y], it contains very large terms, but we know that the result must be between -1 and +1, so there is a very delicate cancellation of huge terms going on. Numerically this is just asking for trouble. The function LegendreP[...] does not simply evaluate a polynomial but rather uses a numerically stable recursion relation, which then turns out to be slower but much better behaved. See, e.g., http://en.wikipedia.org/wiki/Legendre_polynomials at the section "Additional properties of Legendre polynomials" The same is probably true for Bessel functions. Maybe you can implement your own Legendre and Bessel procedures starting from the "Numerical Recipes" book, chapter "Special Functions". Or you could call the GNU Scientific Library through MathLink http://www.gnu.org/software/gsl/ But I doubt you'll do better than the Mathematica built-in numerical routines, except if you need a whole series of LegendreP[n,m,x] for fixed x but varying n,m (the recursion relations tend to produce not only the LegendreP[n,m,x] you are asking for, but all those with smaller indices as well). Or you can just define P[y] etc. through SetDelayed (:=), as Norbert Marxer points out. Or do some steps of your calculation analytically instead of numerically, if possible. Cheers! Roman On Apr 16, 10:06 am, "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