Physics 410 Tutorial 2 - Polynomial Interpolation: Solution

1 Solution from previous tutorial

1.1 barycentric.m

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

1.2 tbarycentric.m

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');

1.3 Plotting output from 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.