function interpolateRunge(max_degree)
% interpolateRunge: Plots some polynomial interpolants of Runge's function.
%
%   interpolateRunge(max_degree)
%
% This function constructs some polynomial interpolants of Runge's famous
% example (taken from Heath section 7.3)
%
%        f(x) = 1 / (1 + 25x^2)
%
% The monomial basis is used, although any other basis could have been used.
% The results are plotted (nothing is returned).  The interpolants start
% with degree zero (a constant) and work their way up to degree max_degree.
% Interpolation points are chosen by linspace in the range [ -1, +1 ] for
% each n.
%
% Also plotted is the error in the interpolant.
%
% Input Parameters:
%
%   max_degree - Maximum degree of polynomial interpolant.  Default = 10.
%
% Output Parameters:
%
%   None.
%
% Ian Mitchell for CPSC 303 (January 24, 2007).

  % Assign default inputs.
  if(nargin < 1)
    max_degree = 10;
  end
  
  % What is our function?  Use element-wise divide and multiply because we
  % want to compute f for multiple x points all in one function call.
  f = inline('1 ./ (1 + 25 * x.^2)', 'x');

  % Over what range will be choose abscissae?
  x_range = [ -1, +1 ];
  
  % Over what range do we want to display the results.  Choose something
  % a little larger than the interpolation range, and a little larger
  % than the function's range.
  plot_range_x = x_range + [ -0.25, +0.25 ];
  % We should really evaluate f to figure out this value, but I'll
  % hardcode it for simplicity.
  plot_range_y = [ -0.1, +2.1 ];
  
  % Get a new figure window.
  figure;
  
  % Loop over the polynomial degrees.  
  for n = 0 : max_degree

    % Interpolation abscissae.
    % Apply transpose to get column vector.
    if(n == 0)
      x = 0;
    else
      x = linspace(x_range(1), x_range(2), n + 1)';
      % Wait for the user to digest the last plot.
      pause;
    end
    
    % Interpolation data values.
    y = feval(f, x);

    % Basis function matrix.  Note that Matlab's vander() function
    % returns the columns in the opposite order than that used in class.
    A = vander(x);
    
    % Find coefficients.
    c = A \ y;
    
    % Range of interest for the plot.
    plot_x = linspace(plot_range_x(1), plot_range_x(2), 401);
    
    % Clear the figure window.
    clf;
    
    % Plot the function we are trying to match.
    subplot(2,1,1);
    f_y = feval(f, plot_x);
    hf = plot(plot_x, f_y, 'b-');
    hold on;
    
    % Plot the interpolation points.
    hp = plot(x, y, 'r*');
    
    % Also plot the interpolant between interpolation points, to compare
    % results.
    plot_y = polyval(c, plot_x);
    hi = plot(plot_x, plot_y, 'r:');

    % Label the plot.
    ylabel('f(x)');
    title('Monomial interpolation example - Ian Mitchell for CPSC 303');
    legend('Runge function', [ 'n = ' num2str(n) ' interpolant' ]);
    
    % Easier to see lines.
    set([ hf; hi ], 'LineWidth', 2);
    set(hp, 'MarkerSize', 9);

    % Interpolant can get really crazy, so we'll tighten the plot bounds.
    axis([ plot_range_x, plot_range_y ]);
    
    % Plot the error.  We will only look at magnitude of error, and we will use
    % a separate subplot so that we can examine the error with logarithmic
    % scaling.  Note how error gets much smaller near the interpolation
    % points (it goes to exactly zero at the interpolation points).
    subplot(2,1,2);
    plot_error = abs(feval(f, plot_x) - plot_y);
    he = semilogy(plot_x, plot_error);

    % Matlab seems to use weird axis bounds for this plot, so we will
    % tighten them up.  First, get the current bounds so as to keep the
    % range in the y-axis.
    axis_error = axis;
    % Now adjust the range in the x-axis to agree with the range for the
    % plot above.
    axis_error(1:2) = plot_range_x;
    % And finally, reset the axis.
    axis(axis_error);
    
    % Label the plot.  Both axes share the same label in the x-axis.
    % Only one line in this plot, so y-axis label substitutes for a legend.
    ylabel('error');
    xlabel('x');

    set(he, 'LineWidth', 2);
    
  end