Thread Subject: Legendre Polynomial Evaluation

what is the fastest way to evaluate a legendre polynomial series

 yi = c0P(0,xi) + c1P(1,xi) + ... + cnP(n,xi)

for numel(x) >> 1.

That is, iv already done the least squares to get c0, c1, ... , cn. Now i want to evaluate this for very large number of points xi.

Matlab's built-in function has god awful timing.

thanks!

"Abel Brown" <brown.2179@osu.edu> wrote in message <ggh03q$ms7$1@fred.mathworks.com>...
> what is the fastest way to evaluate a legendre polynomial series
>
> yi = c0P(0,xi) + c1P(1,xi) + ... + cnP(n,xi)
>
> for numel(x) >> 1.
>
> That is, iv already done the least squares to get c0, c1, ... , cn. Now i want to evaluate this for very large number of points xi.
>
> Matlab's built-in function has god awful timing.

You could just use the three term recurrence relation
for the Legendre polynomials. Thus,

  P0 (x) = 1
  P1 (x) = x
  (n+1)P(n+1,x) = (2n+1)x P(n,x)-nP(n-1,x)

A simple loop will suffice to do it very quickly.
For example, assume the coefficients are stored in
the vector C. (Since your coefficients start at c0,
there is an offset.)

Pnm1= ones(size(x));
Pn = x;
Yi = C(1)*Pnm1 + C(2)*Pn;
if N>2
  for n = 2:N
    % 3 term recurrence
    Pnp1 = ((2*n+1)*x.*Pn - n*Pnm1)/(n+1);
    Yi = Yi + C(n+1)*Pnp1; % zero base for C index
    Pnm1= Pn;
    Pn = Pnp1;
  end
end

HTH,
John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <ggh2en$sos$1@fred.mathworks.com>...
> "Abel Brown" <brown.2179@osu.edu> wrote in message <ggh03q$ms7$1@fred.mathworks.com>...
> > what is the fastest way to evaluate a legendre polynomial series
> >
> > yi = c0P(0,xi) + c1P(1,xi) + ... + cnP(n,xi)
> >
> > for numel(x) >> 1.
> >
> > That is, iv already done the least squares to get c0, c1, ... , cn. Now i want to evaluate this for very large number of points xi.
> >
> > Matlab's built-in function has god awful timing.
>
> You could just use the three term recurrence relation
> for the Legendre polynomials. Thus,
>
> P0 (x) = 1
> P1 (x) = x
> (n+1)P(n+1,x) = (2n+1)x P(n,x)-nP(n-1,x)
>
> A simple loop will suffice to do it very quickly.
> For example, assume the coefficients are stored in
> the vector C. (Since your coefficients start at c0,
> there is an offset.)
>
> Pnm1= ones(size(x));
> Pn = x;
> Yi = C(1)*Pnm1 + C(2)*Pn;
> if N>2
> for n = 2:N
> % 3 term recurrence
> Pnp1 = ((2*n+1)*x.*Pn - n*Pnm1)/(n+1);
> Yi = Yi + C(n+1)*Pnp1; % zero base for C index
> Pnm1= Pn;
> Pn = Pnp1;
> end
> end
>
> HTH,
> John

Thanks John!

just a quick question shouldn't n go from 1:N-1? otherwise P2(x) gets skipped ... ?

"Abel Brown" <brown.2179@osu.edu> wrote in message <gghg7q$t21$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <ggh2en$sos$1@fred.mathworks.com>...
> > "Abel Brown" <brown.2179@osu.edu> wrote in message <ggh03q$ms7$1@fred.mathworks.com>...
> > > what is the fastest way to evaluate a legendre polynomial series
> > >
> > > yi = c0P(0,xi) + c1P(1,xi) + ... + cnP(n,xi)
> > >
> > > for numel(x) >> 1.
> > >
> > > That is, iv already done the least squares to get c0, c1, ... , cn. Now i want to evaluate this for very large number of points xi.
> > >
> > > Matlab's built-in function has god awful timing.
> >
> > You could just use the three term recurrence relation
> > for the Legendre polynomials. Thus,
> >
> > P0 (x) = 1
> > P1 (x) = x
> > (n+1)P(n+1,x) = (2n+1)x P(n,x)-nP(n-1,x)
> >
> > A simple loop will suffice to do it very quickly.
> > For example, assume the coefficients are stored in
> > the vector C. (Since your coefficients start at c0,
> > there is an offset.)
> >
> > Pnm1= ones(size(x));
> > Pn = x;
> > Yi = C(1)*Pnm1 + C(2)*Pn;
> > if N>2
> > for n = 2:N
> > % 3 term recurrence
> > Pnp1 = ((2*n+1)*x.*Pn - n*Pnm1)/(n+1);
> > Yi = Yi + C(n+1)*Pnp1; % zero base for C index
> > Pnm1= Pn;
> > Pn = Pnp1;
> > end
> > end
> >
> > HTH,
> > John
>
> Thanks John!
>
> just a quick question shouldn't n go from 1:N-1? otherwise P2(x) gets skipped ... ?

Yes, you are right. I typed too fast and did not verify
my results. Or maybe I was just testing you. Yeah,
thats it, and I'm sticking to the story. ;-)

John