% % Mfile name
%   mtl_inp_sim_spechoiceofpoints.m

% Language: Matlab 7.4.0 (R2007a)
% Revised: 
%   August 8, 2008

% Authors:
%   Nathan Collier, Luke Snyder, Autar Kaw
%   University of South Florida
%   kaw@eng.usf.edu
%   Website: http://numericalmethods.eng.usf.edu
       
% Purpose
%   To illustrate that choice of points can dramatically change the outcome
%   of a higher order polynomial interpolant.

% Keywords
%   Higher order Interpolation
%   Carl Runge

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


% INPUTS: This is the only section in the program where the user can make
% changes to the program.

    % Enter the number of points chosen for interpolation in [-1,1]
    n=13;
    % NOTE: For the number scheme to work properly n needs to be an odd 
    % number.

% *************************************************************************

disp(sprintf('\n\nAdditional Interpolation Topics: The Effect of Choice of Points'))
disp(sprintf('on Interpolation.'))
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 choice of points is a critical factor in higher order'))
disp(sprintf('interpolation.'))

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

disp(sprintf('\nIn 1901, Carl Runge published his work on dangers of higher order'))
disp(sprintf('interpolation. He took a simple looking function f(x)=1/(1+25x^2) on'))
disp(sprintf('the interval [-1,1].  He took points equidistantly spaced in [-1,1]'))
disp(sprintf('and interpolated the points with polynomials.  He found that as he')) 
disp(sprintf('took more points, the polynomials and the original curve differed')) 
disp(sprintf('considerably. However, if he took data points close to the ends of')) 
disp(sprintf('the interval [-1,1], the problem of large differences between '))
disp(sprintf('interpolated and actual values was less pronounced.  This simulation '))
disp(sprintf('demonstrates this phenomena.'))

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

% SOLUTION:

disp(sprintf('\nFigure 1:'))
disp(sprintf('The first plot shows Runge''s function and the corresponding'))
disp(sprintf('polynomial interpolant that is found using equidistant points.'))
disp(sprintf('Note that in this case, the interpolated function is more accurate'))
disp(sprintf('around the middle of the function, but the accuracy decreases'))
disp(sprintf('rapidly as the data points taken move away from the function center.'))
disp(sprintf('Percent relative true error values would be highest at towards the'))
disp(sprintf('ends of the interval.'))

% Creating the first plot which shows a polynomial interpolation with
% equidistantly spaced points.
figure(1)
subplot(2,1,1)
frunge=inline('1/(1+25*x^2)');
fplot(frunge,[-1 1])
title('\bfRunge^,s function and interpolant with equidistantly spaced points','FontSize',10)
xlabel('\bfx values from [-1,1]')
ylabel('\bff(x)')
hold on
xlim([-1,1])
ylim([-1.1,1.1])

% Using a loop and algorithm to calculate equidistant points for x and then
% plugging those values into the inline function "frunge" to obtain the
% actual values for the function.
for i=1:n
    X(i)=2/(n-1)*(i-1)-1;
    Y(i)=frunge(X(i));
end

% Using polynomial interpolation to create a higher order function based on
% the equidistant points.
plot(X,Y,'ro','MarkerSize',6,'MarkerFaceColor',[1,0,0])
p=polyfit(X,Y,numel(X)-1);
x1=-1:0.01:1;
y1=polyval(p,x1);
plot(x1,y1,'k')

% Creating a legend for the plot.
legend('= Runge^,s function','= Equidistant points','= Interpolated Function',3,'Location','South')

hold off

% Calculating the error in this interpolant versus the original function.
f = 1./(1+25.*x1.^2);

% Calculating the difference between the two functions.
F = abs(f-y1);

% Plotting the results.
subplot(2,1,2)
plot(x1,F);
title('\bfError in Equidistant Points versus original function')
xlabel('\bfx values from [-1,1]')
ylabel('\bfError calculated');

% Creating a figure for the second case: Points taken towards ends of the x
% interval.

disp(sprintf('\nFigure 2:'))
disp(sprintf('Note here that in contrast to the case with equidistant points,'))
disp(sprintf('the function is more accurate next to the end points of the interval'))
disp(sprintf('as opposed to the center as in the previous case.  Percent relative')) 
disp(sprintf('true error values would be highest towards the function center versus'))
disp(sprintf('its outer values.'));

figure(2)
subplot(2,1,1)
fplot(frunge,[-1,1])
title('\bfRunge^,s function and interpolant with more points near the end of the interval [-1,1]','FontSize',10)
xlabel('\bfx from -1 to 1')
ylabel('\bff(x)')
hold on
xlim([-1,1])
ylim([-1.1,1.1])

% This loop creates the factor by which the points will be divided based on
% the algorithm shown.
d=0;
for i = 2:(n-1)/2+1
    d=d+2^(i-2);
end
l=1/d;
Xb(1)=-1;

% This loop starts from -1, and uses the factor 'l' to create a new array
% Xb which houses the new data points. Data points are created from [-1,0].
for i = 2:(n-1)/2+1
    Xb(i)=Xb(i-1)+2^(i-2)*l;
end

% This is the other half of the previous loop which creates the remaining
% data points from [0,1] based on the same 'l' factor and algorithm.
for i = (n-1)/2+2:n
    Xb(i)=Xb(i-1)+2^(n-i)*l;
end

% This loop inputs the newly created values from the "Xb" array and enters
% them into the inline function "frunge".
for i=1:n
    Yb(i)=frunge(Xb(i));
end

n = numel(X);
% Plotting the newly created points and using higher order interpolation as
% previously seen to create a polynomial fit for the data.
plot(Xb,Yb,'ro','MarkerSize',6,'MarkerFaceColor',[1,0,0])
p=polyfit(Xb,Yb,n-1);
x1=-1:0.01:1;
y1=polyval(p,x1);
plot(x1,y1,'k')

% Creating a legend for the figure.
legend('= Runge^,s Function','= Created Points','= Interpolated Function',3,'Location','South')

hold off

% Calculating the error between the two functions
F = abs(f-y1);

% Plotting the results
subplot(2,1,2)
plot(x1,F)
title('\bfError in Polynomial with points near the ends versus original function')
xlabel('\bfx values from [-1,1]')
ylabel('\bfError calculated')




