function func = intrunge ( iknot, n );
% IntRunge
%
% This m-file is part of the material of the lectures on
%   Numerical Mathematics
% at the Math Department of Halle University, Germany.
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  Apr 5, 2005
%
% Polynomial interpolation, Function of Runge
%
% Parameters:
%   iknot  (input)  : control flag "knots"
%                       =1 ... equidistant
%                       =2 ... Tschebyscheff
%   n      (input)  : degree of interpolation polynomial
%
% Examples:
%   intrunge ( 1, 2 );   print -dpng ../png/intrunge2102
%   intrunge ( 1, 6 );   print -dpng ../png/intrunge2106
%   intrunge ( 1, 8 );   print -dpng ../png/intrunge2108
%   intrunge ( 2, 6 );   print -dpng ../png/intrunge2206
%   intrunge ( 2, 8 );   print -dpng ../png/intrunge2208

% -> data grid
xmin = -5.0;
xmax =  5.0;
if iknot==1,
    x = linspace ( xmin, xmax, n+1 );
    kstr = sprintf ( 'Knoten: %2i (aequidistant)', n+1 );
elseif iknot==2,
    x = ( xmin + xmax + ( xmax - xmin ) * cos(pi*(2*(1:(n+1))-1)/(2*(n+1))) ) / 2;
    kstr = sprintf ( 'Knoten: %2i (Tschebyscheff)', n+1 );
end;

% -> data points
f = 1 ./ ( 1 + x.^2 );

% -> interpolation polynomial in Newton form: difference scheme
a = zeros ( n+1, 1 ); 
for i=n:-1:0, 
    a(i+1) = f(i+1); 
    for j=i+1:n, 
        a(j+1) = (a(j+1)-a(j)) / (x(j+1)-x(i+1));
    end;
end;

% -> output data
nxx = 1001;
xx  = linspace ( -5, 5, nxx );
ff = 1 ./ ( 1 + xx.^2 );

pp  = zeros ( size(xx) );
for ixx=1:nxx, 
    pp(ixx) = a(n+1); 
    for i=n-1:-1:0, 
        pp(ixx) = pp(ixx) * (xx(ixx)-x(i+1)) + a(i+1);
    end;
end;

% -> plot
subplot ( 1, 2, 1 );
plot ( x, f, '*b', xx, ff, 'g', xx, pp, 'r' )
hold on;
plot ( [ xmin xmax ], [ 0 0 ], ':k' );
plot ( [ 0 0 ], [ -1.0 2.0 ], ':k' );
hold off;
axis ( [ -5.0 5.0 -1.0 2.0 ] );
xlabel ( 'x' );
ylabel ( 'f' );
title ( [ 'f(x) = 1/(1+x^2),',  kstr ] );

subplot ( 1, 2, 2 );
semilogy ( xx, abs(pp-ff), 'b' );
axis ( [ -5.0 5.0 1.0e-6 10.0 ] );
xlabel ( 'x' );
ylabel ( 'Interpolationsfehler' );
title ( [ 'f(x) = 1/(1+x^2),',  kstr ] );
