X-Mozilla-Status: 0001
        X-Mozilla-Status2: 00000000
               Message-ID: <3A59F43C.2C688A0E@umich.edu>
                     Date: Mon, 08 Jan 2001 12:09:16 -0500
                     From: Eric Durant <edurant@umich.edu>
                 X-Mailer: Mozilla 4.7 [en] (WinNT; U)
        X-Accept-Language: en
             MIME-Version: 1.0
               Newsgroups: comp.soft-sys.matlab
                  Subject: Re: lagrange
               References: <978892964.49463@tubarao.ip.pt>
             Content-Type: text/plain; charset=iso-8859-1
Content-Transfer-Encoding: 8bit

Guilherme Santos wrote:
>does anyone have a  .m file with the Lagrange polinomial
>interpolater?

A similar question was posted here in November, 1999.  Here's a summary
of my response...

Here's an algorithm for finding the Lagrange interpolation polynomial
(equivalently, filter) that isn't vectorized, but does eliminate the
redundant multiplications of the straightforward implementation.  Also,
it doesn't require huge matrix storage.  It allows irregularly spaced
values of the independent variable.
    
I've included a test case.  sum([polyval(p,x)-y].^2) approaches 0 as the
roundoff error approaches 0.
    
This function uses a neat, Horner-inspired factorization of the Lagrange
polynomial.

=================================================
function [p] = lagrange(x,y)
% Eric Durant <edurant@umich.edu>
% Wednesday 6 October 1999
%
% Find the interpolation polynomial for f(x) = y
% using the Lagrange interpolation method.

l = length(y);
for i=1:l,   
    alpha(i) = y(i);
    for j=1:i-1,
        alpha(i) = (alpha(i) - alpha(j)) / (x(i) - x(j));
    end
end
ip = 1;
p = zeros(1,l);   
for i=1:l-1,  
    p(l-i+1:l) = p(l-i+1:l) + alpha(i) * ip;
    ip = conv(ip, [1 -x(i)]);
end 
p = p + alpha(l) * ip;

%%% Test case:
% » x = 1:4;
% » y = [2 7 18 41];
% » p = lagrange(x,y)
%
% p =
%      1    -3     7    -3
%   
% » polyval(p,x)
%
% ans =
%      2     7    18    41
=================================================  

Eric Durant
http://edurant.com/