function fto = barycentric(xfr, ffr, xto)
% fto performs polynomial interpolation using barycentric formula.
%
% Input arguments
%
% xfr: x values to fit (length n vector)
% ffr: y values to fit (length n vector)
% xto: x values at which to evaluate fit (length nto vector)
%
% Output argument
%
% yto: Interpolated values (length nto vector)
% Initialize output ...
fto = zeros(1, length(xto));
% Check for consistent length of input arguments ...
if length(xfr) ~= length(ffr)
error('Input arguments xfr and ffr must be of same length');
end
% Calculate barycentric weights ...
n = length(xfr);
w = ones(1,n);
for j = 1 : n
for k = 1 : n
if k ~= j
w(j) = w(j) * (xfr(j) - xfr(k));
end
end
w(j) = 1.0 / w(j);
end
% Perform interpolation to xto values ...
for k = 1 : length(xto)
num = 0;
den = 0;
for j = 1 : n
num = num + w(j) * ffr(j) / (xto(k) - xfr(j));
den = den + w(j) / (xto(k) - xfr(j));
end
fto(k) = num / den;
end
end
clf;
x = linspace(-0.9,0.9,101);
fxct = tpoly(x);
plot(x,fxct);
hold on;
xfr = linspace(-1,1,8);
ffr = tpoly(xfr);
xto = linspace(-0.9,0.9,101);
fto = barycentric(xfr,ffr,xto);
df = fxct - fto
plot(x,1.0e15 .* df);
title('Barycentric polynomial interpolation');
xlabel('x');
legend('Exact values', '10^{15} * Error in interpolated values');
tbarycentric.m
Observe that the error in the interpolated values has been scaled by \( 10^{15} \), which results in values that are \( O(1) \). This indicates that the interpolation error is of the order of machine precision, as we would expect given the relatively high degree of the interpolating polynomial.