% % Mfile name
%   mtl_inp_sim_spechoiceofpoints.m

% Revised: 
%   August 8, 2008

% Authors:
%   Nathan Collier, Luke Snyder, Dr. Autar Kaw
%   University of South Florida
%   kaw@eng.usf.edu
%   Website: http://numericalmethods.eng.usf.edu
       
% Purpose
%   To illustrate the dangers of higher polynomial interpolation as a means
%   of approximation.

% Keywords
%   Interpolation
%   Carl Runge

% Clearing all data, variable names, and files from any other source and
% clearing the command window after each succesive run of the program.
clc
clf
clear all

% Inputs:
%    This is the only place in the program where the user makes the changes
%    based on their wishes.  The user can make three separate choices for
%    the number of equidistant points in [-1,1]. In other words, if n1=5, a
%    4th order polynomial will  be used to approximate 1/(1+25x^2) in
%    [-1,1].

    % Number of equidistant points to use:
    n1 = 5;
    n2 = 13;
    n3 = 18;
    
% *************************************************************************

disp(sprintf('\n\nHigher Order Interpolation is a Bad Idea'))
disp(sprintf('\nUniversity of South Florida'))
disp(sprintf('United States of America'))
disp(sprintf('kaw@eng.usf.edu'))
disp(sprintf('\nNOTE: This worksheet utilizes Matlab to illustrate the concept'))
disp(sprintf('that higher order interpolation is a bad idea as proven by Carl Runge'))
disp(sprintf('in 1901.'));

disp(sprintf('\n*****************************Introduction******************************'))

disp(sprintf('\nThe following simulation illustrates why higher order interpolation'))
disp(sprintf('as a means of approximation is often a bad idea.  In 1901, Carl Runge'))
disp(sprintf('published his work on the danger of higher order polynomial interpolation.'))
disp(sprintf('He used a simple looking function, f(x) = 1/(1+25x^2) on the interval of'))
disp(sprintf('[-1,1], took points equidistantly spaced in that interval and interpolated'))
disp(sprintf('between the points with polynomials.  He found that as he took more points,'))
disp(sprintf('the polynomials and the original curve differed considerably in their value.'))
disp(sprintf('This resulted in a large propagation of error as the order of the polynomials'))
disp(sprintf('increased, leading him to conclude that higher order polynomial interpolation'))
disp(sprintf('will approximate very inaccurate results.'))

disp(sprintf('\n\n******************************Input Data*********************************'))
disp(sprintf('\nThe first approximation will use, n1 = %g points',n1));
disp(sprintf('\nThe second approximation will use n2 = %g points',n2));
disp(sprintf('\nThe third approximation will use n3 = %g points',n3));

disp(sprintf('\n******************************Simulation*********************************'))    

% SOLUTION:

% Plotting Runge's Function
frunge=inline('1./(1+25.*x.^2)');
n=1;

% Vectorizing the x data across the interval [1,1].
x = -1:0.01:1;

figure(1)
plot(x,frunge(x),'k','LineWidth',3)
title('Runge^,s function and three interpolants using equidistantly spaced points','Fontweight','bold','FontSize',10)
xlabel('\bfValue of x from [-1,1]');
ylabel('\bff(x)')
hold on

% Performing each of the three polynomial fits
% n1 polynomial fit
for i=1:n1
    % This algorithm is used to create an array of equidistant points based
    % on the number of points to be used as specified by the user.
    x1(i)=(2/(n1-1))*(i-1)-1;
    y1(i)=frunge(x1(i));
    
end

p1=polyfit(x1,y1,n1-1);
x1p=-1:0.01:1;
y1p=polyval(p1,x1p);
plot(x1p,y1p,'b','LineWidth',3)

disp(sprintf('\nFor n1, %g points are used and so the first polynomial fit is of order %g.',n1,n1-1));

% n2 polynomial fit
for i=1:n2
    % This algorithm is used to create an array of equidistant points based
    % on the number of points to be used as specified by the user.
    x2(i)=(2/(n2-1))*(i-1)-1;
    y2(i)=frunge(x2(i));
    
end

p2=polyfit(x2,y2,n2-1);
x2p=-1:0.01:1;
y2p=polyval(p2,x2p);
plot(x2p,y2p,'r','LineWidth',3)

disp(sprintf('\nFor n2, %g points are used and so the second polynomial fit is of order %g.',n2,n2-1));

% n3 polynomial fit
for i=1:n3
    % This algorithm is used to create an array of equidistant points based
    % on the number of points to be used as specified by the user.
    x3(i)=(2/(n3-1))*(i-1)-1;
    y3(i)=frunge(x3(i));
    
end

p3=polyfit(x3,y3,n3-1);
x3p=-1:0.01:1;
y3p=polyval(p3,x3p);
plot(x3p,y3p,'g','LineWidth',3)

disp(sprintf('\nFor n3, %g points are used and so the third polynomial fit is of order %g.',n3,n3-1));

% Using concactenation to display the order of the polynomials used in the
% simulation.
strn1 = ['n1 = ',num2str(n1)];
strn2 = ['n2 = ',num2str(n2)];
strn3 = ['n3 = ',num2str(n3)];

legend('Runge^,s Function',strn1,strn2,strn3,4)
ylim([-1.5,1.5])
xlim([-1,1])

% Constructing and printing the three tables of values
% n1's table:
disp(sprintf('\n\n**************************Table of Values for n1**************************'));
disp('x values           Actual Value            Interpolated         %True Error');
disp('---------------------------------------------------------------------------');

for i=1:1:2*(n1-1)+1
    % Creating a new 'xi' array that will be used to interpolate the data
    % in between the specified points based on the polynomial found using
    % those same points.
    xi(i)=2/(2*(n1-1))*(i-1)-1;
    
    % Plugging in values of 'xi' into Runge's function and the interpolated
    % polynomial function.
    M(i,1)=xi(i);
    M(i,2)=frunge(xi(i));
    M(i,3)=polyval(p1,xi(i));
    M(i,4)=abs((frunge(xi(i))-polyval(p1,xi(i)))/frunge(xi(i)))*100;
    
    % If the True Error is found to be less than 10^-7, then the value is
    % considered to be 0.
    if M(i,4)<1*10^-7
        M(i,4)=0;
    end    

    % Displaying the data in the table.
    if xi(i) == 0
        disp('0                   1                        1                     0');
    else
        str = '%1.3e         %1.3e               %1.3e           %1.3e';
        disp(sprintf(str,M(i,1),M(i,2),M(i,3),M(i,4)));
    end
end

% Constructing a table of values for n2
disp(sprintf('\n\n**************************Table of Values for n2**************************'));
disp('x values           Actual Value            Interpolated         %True Error');
disp('---------------------------------------------------------------------------');

for i=1:1:2*(n2-1)+1
    % Creating a new 'xi' array that will be used to interpolate the data
    % in between the specified points based on the polynomial found using
    % those same points.
    xi(i)=2/(2*(n2-1))*(i-1)-1;
    
    M(i,1)=xi(i);
    M(i,2)=frunge(xi(i));
    M(i,3)=polyval(p2,xi(i));
    M(i,4)=abs((frunge(xi(i))-polyval(p2,xi(i)))/frunge(xi(i)))*100;
    
    % If the True Error is found to be less than 10^-7, then the value is
    % considered to be 0.
    if M(i,4)<1*10^-7
        M(i,4)=0;
    end    

    % Displaying the data in the table.
    if xi(i) == 0
        disp('0                   1                        1                     0');
    else
        str = '%1.3e         %1.3e               %1.3e           %1.3e';
        disp(sprintf(str,M(i,1),M(i,2),M(i,3),M(i,4)));
    end
end

% Constructing a table of values for n3.
disp(sprintf('\n\n**************************Table of Values for n3**************************'));
disp('x values           Actual Value            Interpolated         %True Error');
disp('---------------------------------------------------------------------------');

for i = 1 : 2*(n3-1)+1
    % Creating a new 'xi' array that will be used to interpolate the data
    % in between the specified points based on the polynomial found using
    % those same points.
    xi(i)=2/(2*(n3-1))*(i-1)-1;
    
    M(i,1)=xi(i);
    M(i,2)=frunge(xi(i));
    M(i,3)=polyval(p3,xi(i));
    M(i,4)=abs((frunge(xi(i))-polyval(p3,xi(i)))/frunge(xi(i)))*100;
    
    % If the True Error is found to be less than 10^-7, then the value is
    % considered to be 0.
    if M(i,4)<1*10^-7
        M(i,4)=0;
    end

    % Displaying the data in the table.
    if xi(i) == 0
        disp('0                   1                        1                     0');
    else
        str = '%1.3e         %1.3e               %1.3e           %1.3e';
        disp(sprintf(str,M(i,1),M(i,2),M(i,3),M(i,4)));
    end

end

