%
%                   Polynomial Interpolation
%
%     An M-file is provided to carry out interpolation on a specified
%     interval using evenly spaced points or Chebyshev nodes. The initial
%     settings give the polynomial of degree 8 for 1/(1+25*x^2) on [-1,1],
%     for evenly spaced nodes. Menus allow you to alter the function, the
%     degree of interpolation, the interpolation interval or the type of
%     nodes, and also give you control over the plots produced.
%
%     *****   Important note on entering functions   *****
%
%     When entering a new function you need to know a peculiarity of the
%     MATLAB notation. The operators *, ^ and / denote matrix (rather than
%     scalar) multiplication, exponentiation and (suitably defined) division.
%     To enforce such operations on individual elements we have to use .*,
%     .^ and ./ i.e. put a dot before each operator. To avoid confusion with
%     a decimal point, put a space before the dot if it follows a digit.
%     For example, the default function is entered as 1 ./(1+25 .*x.^2).
%
%                                                      J.F.Blowey
%                                                  University of Durham
%     Press any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
help interp
pause
fct='1 ./(1+25 .*x.^2)'; a=-1; b=1; deg=8; nodetype=1; changed=1;hold on;
clc;
choise = 0; mlen = 6;
ms = 'Polynomial Interpolation Menu';
ms2 = 'function to be approximated';
ms3 = 'interpolation interval';
ms4 = 'degree of interpolation';
ms5 = 'node type - equal intervals or Chebyshev nodes';
ms6 = 'begin calculation';
ms_mlen = 'Quit';
while (choise ~= mlen)
  choise = menu(ms,ms2,ms3,ms4,ms5,ms6,ms_mlen); 
  fprintf(' \n') 
%
  if (choise == 1)        % function 
      disp('Enter function and press RETURN ');
      fct=input(' ','s'); changed=1;
  elseif (choise == 2)      
      a=input('Enter left-hand limit  ');
      b=input('Enter right-hand limit '); changed=1;
  elseif (choise == 3)      
      deg=input('Enter degree of interpolating polynomial (<50) ');
      if (deg>=50) disp('The numerical results you see may suffer from rounding error');end
  elseif (choise == 4)      
      nodetype=input('Enter 1 for even spacing, 2 for Chebyshev nodes  ');
  elseif (choise == 5)      
clc;
      npts=deg+1; nsampl=101;

      if nodetype==1
        xnodes=a+(0:npts-1)*((b-a)/(npts-1));
      else
        nodetype=2;
      xnodes=(a+b)/2+((b-a)/2)*cos(pi/(deg+1)*((0:deg)'+.5));
      end

      xsampl=a+(0:nsampl-1)*((b-a)/(nsampl-1));
      x=xnodes; ydata=eval(fct);
      % Recomputation avoided when neither function nor interval changes
      if changed >0
         x=xsampl; yexact=eval(fct); changed=0;
      end;
      for k=1:(length(xnodes)-1), for i=1:(length(xnodes)-k)
      ydata(i)=(ydata(i+1)-ydata(i))./(xnodes(i+k)-xnodes(i));
      end, end, cof=ydata;
      yinterp=(cof(1)*ones(length(xsampl),1))';
      for i=2:length(xnodes); yinterp=cof(i)+(xsampl-xnodes(i)).*yinterp; end

      disp(['The function interpolated is ',fct]);
      disp(['Interpolation interval [',num2str(a),',',num2str(b),']']);
      disp(['Polynomial degree = ',num2str(deg)]);
      if (nodetype==1) 
        disp(['Evenly spaced nodes']);
      else
       disp(['The Cheybshev nodes']);
      end
      disp('');

choice = 0; menlen = 5;
m1 = 'Plot menu for interpolation';
m2 = 'plot interpolating polynomial';
m3 = 'plot exact function values';
m4 = 'plot error in current interpolation';
m5 = 'clear graphics screen';
m_menlen = 'Back';
while (choice ~= menlen)
  choice = menu(m1,m2,m3,m4,m5,m_menlen); 
%  fprintf(' \n') 
  hue(1)='r'; hue(2)='g'; % plotting colours
  if choice==1
    hold on;plot(xsampl,yinterp,hue(nodetype)); 
  elseif choice==2
   hold on;plot(xsampl,yexact,'b');
  elseif choice==3
    yerror=yexact-yinterp;
    hold on;plot(xsampl,yerror,['--',hue(nodetype)]);  
  elseif choice==4
    hold off;clf;
  end
end
end                  % if, elseif
%
end                    % while choise ~= mlen

clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% File interp.m                                 January 5th 1999
%
% Author:   J.F.Blowey
%           Dept. of Mathematical Sciences,
%           University of Durham,
%           South Road, Durham  DH1 3LE, England.
%
%   e-mail: J.F.Blowey @ durham.ac.uk
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
