Matlab Lecture Examples from Week 3

Random Cubic Interpolants of sin(x)

close all
x0 = linspace(0,2*pi,100)';
y0 = sin(x0);
for eg=1:10
   x = 2*pi*sort(rand(4,1));
   y = sin(x);
   c = InterpN(x,y);
   pvals = HornerN(c,x,x0);
   plot(x0,y0)
   axis([0 2*pi -2 2])
   hold on
   plot(x,y,'*r')
   pause
   plot(x0,pvals,'r')
   pause
   hold off
end

close all
x0 = linspace(0,2*pi,100)';
y0 = sin(3*x0.^1.1).*exp(-.3*x0);
for n= [2:20 30 40 50]
   x = linspace(0,2*pi,n);
   y = sin(3*x.^1.1).*exp(-.3*x);
   c = InterpN(x,y);
   pvals = HornerN(c,x,x0);
   plot(x0,y0,[-1 8],[0 0],'-')
   axis([-1 8 -1 1])
   hold on
   plot(x,y,'*r')
   plot(x0,pvals,'r')
   text(2,.8,sprintf('n = %1d',n),'fontsize',24)
   text(2,.6,sprintf('Max Error = %7.4f',max(abs(pvals-y0))),'fontsize',24)
   pause
   hold off
end

  function c = InterpN(x,y)
% c = InterpN(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector. c is  a column n-vector with the property that if
%
%      p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
%      p(x(i)) = y(i), i=1:n.

n = length(x);
for k = 1:n-1
   y(k+1:n) = (y(k+1:n)-y(k)) ./ (x(k+1:n) - x(k));
end
c = y;

  function pVal = HornerN(c,x,z)
% pVal = HornerN(c,x,z)
% Evaluates the Newton interpolant on z where
% c and x are n-vectors and z is an m-vector.
%
% pVal is a vector the same size as z with the property that if
%
%         p(x) = c(1) +  c(2)(x-x(1))+ ... + c(n)(x-x(1))...(x-x(n-1))
% then 
%         pval(i) = p(z(i)) , i=1:m.

n = length(c); 
pVal = c(n)*ones(size(z));
for k=n-1:-1:1
   pVal = (z-x(k)).*pVal + c(k);
end