On 18 abr, 10:33, AES <sieg... at stanford.edu> wrote:
> In article <evvaqs$95... at smc.vnet.net>, "Antonio" <ane... at gmail.com>
> wrote:
>
> > 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
>
> Don't know if this will be of any help, but some years ago I was
> evaluating numerical values of Legendre polynomials for high n (in the
> 50s to 60s) and encountered round-off errors (or what looked like
> round-off errors any way) with real arguments, that disappeared when I
> used rational integer or exact arguments.
>
> Example: I wanted to vary the argument by very small steps around Pi.
> Using Pi, 1.001 Pi, 1.002 Pi, etc gave an obviously noisy sequence of
> values. Using Pi, (1001 Pi/1000), (1002 Pi/1000) gave a much smoother
> sequence. I assumed it had something to do with this being a
> polynomial, and Mathematica switching to some kind of integer arithmetic (???).
Well if I try to evaluate LegendreP numerically by recursion relation,
not only will it be slow but also carry on any round offs up to the
high n's I am interested in. Since I need all {n,m}'s up to the big
ones none, I think it is best to have the extended function
(polynimial) determined by Rodrigues. In this case the polynomials
will carry big integer numbers as pointed out by Roman. But in this
case I guess that Mathematica much better in crunching big numbers
than determining LegendreP[n,m,x] for big n's. In this case, is it
better for Mathematica numerically to have a polynomial with Integers
or Reals?
>From AES post it seems that a rational integer argument would be
better than exact arguments. But unfortunatelly this is not the case
for me. I need to use the associated Legendre polynomial inside a
numerical integral, so I can't control the points sampled by
NIntegrate, or can I?
Anyway up to now I will share below what is a good alternative to the
LegendreP[n,m,x] when all functions up to Nmax is needed.
(*Maximum number of precomiled functions*)
Nmax = 70;
(*Build a Table of associated functions through Rodrigues *)
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, Nmax}, {m, 0,
n}];
(*Optimize the obtained expression, to reduce the number of operations
*)
PnmO = Map[Experimental`OptimizeExpression, Pnm, {2}];
(*Compile so that it executes faster*)
PnmC = Map[(Compile[{{x, _Real}}, #]) &, PnmO, {2}];
(*A wrap - up to define P in a function style *)
P[n_, m_, x_?NumericQ] := PnmC[[n, m + 1]][x];
Anyone with a oneliner solution to the above?
Thanks to all,
Antonio