% compare continuous polynomial interpolation with piecewise linear interpolation
% CAAM 453/553, Rice University

 maxdeg = 40;
 xx = linspace(-5,5,1000);                  % fine grid
 fxx = (xx+sin(3*xx)).*exp(-(xx.^2)/6);     % true function value on a fine grid
 err_poly   = zeros(maxdeg+1,1); err_pwpoly = zeros(maxdeg,1);

% continuous polynomial approximation
 for n=0:maxdeg
    x = linspace(-5,5,n+1);
    fx = (x+sin(3*x)).*exp(-(x.^2)/6);
    p = polyfit(x,fx,n);               % MATLAB's polynomial interpolation routine
    err_poly(n+1) = max(abs(polyval(p,xx)-fxx));
    figure(1), clf
    plot(xx,fxx, 'b-','linewidth',2), hold on
    plot(xx,polyval(p,xx), 'r--','linewidth',2)
    plot(x,fx, 'r.','markersize',30)
    axis equal, axis([-5 5 -3 3])
    title(sprintf('Continuous Polynomial Interpolant of Degree %d',n),'fontsize',20)
    input('');
 end

% piecewise linear interpolation
 for n=1:maxdeg
    x = linspace(-5,5,n+1)';
    fx = (x+sin(3*x)).*exp(-(x.^2)/6);
    pxx = interp1(x,fx,xx,'linear');    % MATLAB's piecewise linear interpolation routine
    err_pwpoly(n) = max(abs(pxx-fxx));
    figure(2), clf
    plot(xx,fxx, 'b-','linewidth',2), hold on
    plot(xx,pxx, 'r--','linewidth',2)
    plot(x,fx, 'r.','markersize',30)
    axis equal, axis([-5 5 -3 3])
    title(sprintf('Piecewise Polynomial Interpolant with %d Intervals',n),'fontsize',20)
    input('');
 end

% compare the errors
 figure(3), clf
 semilogy([0:maxdeg], err_poly, 'b-','linewidth', 2), hold on
 semilogy([1:maxdeg], err_pwpoly, 'r--','linewidth', 2)
 legend('Continuous Polynomial', 'Piecewise Linear Polynomials',3)
 title('Maximum Error in Interpolants', 'fontsize', 20)
 xlabel('Number of Data Points, n', 'fontsize', 18)
 ylabel('Maximum Error on [-5, 5]', 'fontsize', 18)