Code covered by the BSD License

# Compute the value of Lagrange polynomial at given points

### Hongkai Dai (view profile)

04 May 2012 (Updated )

The function computes the value of a Lagrange polynomial with mesh points x and mesh values y.

lagrangePolyVal(t,x,y)
```function p = lagrangePolyVal(t,x,y)
% Use Barycentric Formula to compute the value of the lagrange polynomial,
% which is a polynomial interpolation at points x with values y. The code is
% vectorized, please refer to the paper 'Barycentric Lagrange
% Interpolation' by Jean-Paul Berrut and Lloyd N Trefethen for more
% information
% t     -- a row vector, compute the value of the polynomial at t(i)
% x     -- a row vector, the mesh points
% y     -- a matrix, each row of y, y(i,:), represents i'th polynomial.
%          y(i,j) is the value of the i'th polynomial at the mesh point x(j)
% p     -- a matrix, p(i,j) is the value of the i'th lagrange polynomial at
%          t(j)
% For example, if we have two polynomials p1 = x^3+x+1 and p2 = x^2-x, the
% mesh points are x = [1 3 5 7], the mesh values for p1 is y(1,:) = [3 31 131 351], the
% mesh values for p2 is y(2,:) = [0 6 20 42], we want to compute the value of p1 and
% p2 at the points t = [1 2 3 4], then p = lagrangePolyVal(t,x,y) will
% return as p(1,:) = [3 11 31 69] and p(2,:) = [0 2 6 12], exactly the
% values of p1 and p2 computed at t.
t = t(:)';
x = x(:)';
nX = length(x);
nY = size(y,1);
if(size(y,2)~= nX)
error('The size of the mesh points x and the mesh values y do not match')
end
nt = length(t);
% check if some t resides at the mesh points. If so, compute seperately
tx = bsxfun(@times,t',ones(1,nX));
xt = bsxfun(@times,x,ones(nt,1));
[iEq,jEq] = find(abs(tx-xt)<eps);
p = zeros(nY,nt);
p(:,iEq) = y(:,jEq);
if(length(iEq) == nt)
% If all t are mesh points, return the corresponding mesh values;
return;
else
ind = 1:nt;
iNeq = ind(logical(prod(tx-xt,2)'));
tn = t(iNeq);
ntn = length(tn);
xx = bsxfun(@times,x',ones(1,nX));
% omega is the barycentric weights
omega = ones(1,nX)./prod(xx-xx'+eye(nX),2)';
% den is the denominator of the barycentric formula
den = bsxfun(@times,omega',ones(1,ntn))./(bsxfun(@times,tn,ones(nX,1))-bsxfun(@times,x',ones(1,ntn)));
% To multiply with y, repmat the denominator [den(1,:);den(2,:);...] to
% [den(1,:);den(1,:);...den(1,:);den(2,:);den(2,:);...]
if(nY>1)
den = reshape(permute(reshape(repmat(reshape(permute(reshape(den',ntn,1,nX),[2,1,3]),1,[]),nY,1),nY,ntn,nX),[2,1,3]),ntn,nX*nY)';
end
nom = den.*bsxfun(@times,y(:),ones(1,ntn));
p(:,iNeq) = reshape((sum(permute(reshape(nom,nY,nX,ntn),[2 3 1]),1)./sum(permute(reshape(den,nY,nX,ntn),[2 3 1]),1)),ntn,nY)';
end
end
```

Contact us