% % Mfile name
%   mtl_inp_sim_ndd.m

% Language: Matlab 7.4.0 (R2007a)
% Revised: 
%   July 7, 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 the Newton Divided Difference method of interpolation
%   using Matlab.

% Keywords
%   Interpolation
%   Newton Divided Difference

% 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.

% Enter arrays of x and y data and the value of x at which y is desired.

    % Array of x-data
    x=[10 0 20 15 30 22.5];
    
    % Array of y-data
    y=[227.04 0 517.35 362.78 901.67 602.97];
    
    % Value of x at which y is desired
    xdesired = 16;
    
% *************************************************************************

disp(sprintf('\n\nNewton Divided Difference Method of Interpolation'))
disp(sprintf('\nUniversity of South Florida'))
disp(sprintf('United States of America'))
disp(sprintf('kaw@eng.usf.edu'))
disp(sprintf('\nNOTE: This worksheet illustrates the use of Matlab to illustrate'))
disp(sprintf('the concepts of the Newton Divided Difference method of interpolation.'))

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

disp(sprintf('\nThe following simulation uses Matlab to illustrate the Newton Divided'))
disp(sprintf('Difference Method of interpolation.  Given n data points of y vs x,  it'))
disp(sprintf('is required to find the value of y at a particular value of x using'))
disp(sprintf('first, second, or third order interpolation.  It is necessary to first'))
disp(sprintf('pick the needed data points and use those specific points to interpolate'))
disp(sprintf('the data.  This method differs from the Direct Method of interpolation'))
disp(sprintf('by first finding the value of the function at x0, and then solving for'))
disp(sprintf('the slope by using the other known points that bracked the value of'))
disp(sprintf('"xdesired".  If higher order interpolations are required, the coefficients'))
disp(sprintf('for the polynomial functions can be calculated recursively.  Once the'))
disp(sprintf('coefficients are found, the value of f(x) can be approximated at the'))
disp(sprintf('value of "xdesired".'))

disp(sprintf('\n\n***************************Input Data******************************'))
fprintf('\n');
disp('x array of data:')
x
disp('y array of data:')
y
disp(sprintf('The value of x at which y is desired, xdesired = %g',xdesired))

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

% Determining whether or not "xdesired" is a valid point to ask for given
% the range of the x data.
high = max(x);
low = min(x);
if xdesired > high
    disp(sprintf('\nThe value entered for "xdesired" is too high. Please pick a smaller value'))
    disp(sprintf('that falls within the range of x data.'))
    break;
elseif xdesired < low
    disp(sprintf('\nThe value entered for "xdesired" is too low.  Please pick a larger value'))
    disp(sprintf('that falls within the range of x data.'))
    break;
else
    disp(sprintf('\nThe value for "xdesired" falls within the given range of xdata.  The'))
    disp(sprintf('simulation will now commence.'));

% The following considers the x and y data and selects the two closest points to xdesired
% that also bracket it.
n = numel(x);
comp = abs(x-xdesired);
c=min(comp);
for i=1:n
    if comp(i)==c; 
    ci=i;    
    end
end

if x(ci) < xdesired
    q=1;
    for i=1:n
        if x(i) > xdesired
            ne(q)=x(i);
            q=q+1;
        end
    end
    
    b=min(ne);
    for i=1:n
        if x(i)==b
            bi=i;
        end
    end
end

if x(ci) > xdesired
    q=1;
    for i=1:n
        if x(i) < xdesired
            ne(q)=x(i);
            q=q+1;
        end
    end
    b=max(ne);
    for i=1:n
        if x(i)==b
            bi=i;
        end
    end
end

% If more than two values are needed, the following selects the subsequent values and puts
% them into a matrix, preserving the orginal data order.
for i = 1:n
    A(i,2)=i;
    A(i,1)=comp(i);
end
A=sortrows(A,1);

for i=1:n
    A(i,3)=i;
end
A=sortrows(A,2);
d=A(1:n,3);

if d(bi)~=2
    temp=d(bi);
    d(bi)=1;
    for i=1:n
        if i ~= bi & i ~= ci & d(i) <= temp
            d(i)=d(i)+1;
        end
        d(ci)=1;
    end
end

%%%%%%%%% LINEAR INTERPOLATION %%%%%%%%%

% Pick two data points
datapoints=2;
p=1;
for i=1:n
    if d(i) <= datapoints
        xdata(p)=x(i);
        ydata(p)=y(i);
        p=p+1;
    end
end

% Calculating coefficients of Newton's Divided difference polynomial
b0=sym('b0');
b1=sym('b1');
z=sym('z');
b0=ydata(1);
b1=(ydata(2)-ydata(1))/(xdata(2)-xdata(1));
fl=b0 +b1*(z-xdata(1));
fxdesired=subs(fl,z,xdesired);
fprev=fxdesired;

% Displaying the outputs:
disp(sprintf('\nLINEAR INTERPOLATION:'))
disp(sprintf('\nx data chosen: x1 = %g, x2 = %g',xdata(1),xdata(2)))
disp(sprintf('\ny data chosen: y1 = %g, y2 = %g',ydata(1),ydata(2)))

% Displaying the determined coefficients b0,b1.
disp(sprintf('\nCoefficients of the linear interpolation, b0 = %g, b1 = %g',b0,b1))

% Using concactenation to properly display the final function in the
% command window.
fprintf('\n');
str1 = ['Final Function: f1(x) = ',num2str(b0)];
if b1 > 0
    str2 = [' + ',num2str(b1),'(x'];
else
    str2 = [' ',num2str(b1),'(x'];
end    

if xdata(1) > 0
    str3 = [' - ',num2str(xdata(1)),')'];
else
    str3 = [' + ',num2str(xdata(1)),')'];
end

% Putting all the strings together.
finalstr = [str1,str2,str3];
disp(finalstr);

% Displaying the final value of the function at xdesired.
disp(sprintf('\nFinal value at xdesired = %g, f(%g) = %g',xdesired,xdesired,fxdesired))

% Creating inline function for plotting.
func = inline([num2str(b0),'+',num2str(b1),'*(z-',num2str(xdata(1)),')']);

% Figures displaying the interpolation and data points.
axis on
figure(1)
subplot(2,1,1)
fplot(func,[min(xdata),max(xdata)])
hold on
title('Linear Interpolation (Data Points Used)','Fontweight','bold')
xlabel('x data','Fontweight','bold');
ylabel('y data','Fontweight','bold');
plot(xdata,ydata,'ro','MarkerSize',8','MarkerFaceColor',[1,0,0]);
plot(xdesired,fxdesired,'kx','Linewidth',2,'MarkerSize',12');


% Creating the second figure.
axis on
subplot(2,1,2)
fplot(func,[min(xdata),max(xdata)])
hold on
title('Linear Interpolation (Full Data Set)','Fontweight','bold')
xlabel('x data','Fontweight','bold');
ylabel('y data','Fontweight','bold');
plot(x,y,'ro','MarkerSize',8','MarkerFaceColor',[1,0,0])
plot(xdesired,fxdesired,'kx','Linewidth',2,'MarkerSize',12')
xlim([min(x) max(x)])
ylim([min(y) max(y)])
hold off


%%%%%%%%% QUADRATIC INTERPOLATION %%%%%%%%%

% Pick three data points
datapoints=3;
p=1;
for i=1:n
    if d(i) <= datapoints
        xdata(p)=x(i);
        ydata(p)=y(i);
        p=p+1;
    end
end

% Calculating coefficients of Newton's Divided difference polynomial
b0=sym('b0');
b1=sym('b1');
b2=sym('b2');
z=sym('z');
b0=ydata(1);
b1=(ydata(2)-ydata(1))/(xdata(2)-xdata(1));
b2=((ydata(3)-ydata(2))/(xdata(3)-xdata(2))-(ydata(2)-ydata(1))/(xdata(2)-xdata(1)))/(xdata(3)-xdata(1));
fq=b0 +b1*(z-xdata(1)) + b2*(z-xdata(1))*(z-xdata(2));
fxdesired=subs(fq,z,xdesired);

fnew=fxdesired;
ea=abs((fnew-fprev)/fnew*100);
if ea >= 5
    sd=0;
else
    sd=floor(2-log10(abs(ea)/0.5));
end

% Displaying the outputs:
disp('---------------------------------------------------------------------')
disp(sprintf('\nQUADRATIC INTERPOLATION:'))
disp(sprintf('\nx data chosen: x1 = %g, x2 = %g, x3 = %g',xdata(1),xdata(2),xdata(3)))
disp(sprintf('\ny data chosen: y1 = %g, y2 = %g, y3 = %g',ydata(1),ydata(2),ydata(3)))

% Displaying the determined coefficients b0, b1 and b2.
disp(sprintf('\nCoefficients of the quadratic interpolation: b0 = %g, b1 = %g, b2 = %g',b0,b1,b2))

% Using concactenation to properly display the final function in the
% command window.
fprintf('\n')
str1 = ['Final Function: f2(x) = ',num2str(b0)]; 
if b1 > 0
    str2 = [' + ',num2str(b1),'(x'];
else
    str2 = [' ',num2str(b1),'(x'];
end    

if xdata(1) > 0
    str3 = [' - ',num2str(xdata(1)),')'];
    str5 = [str3,'(x'];
else
    str3 = [' + ',num2str(xdata(1)),')'];
    str5 = [str3,'(x'];
end

if b2 > 0
    str4 = [' + ',num2str(b2),'(x'];
else
    str4 = [' ',num2str(b2),'(x'];
end

if xdata(2) > 0
    str6 = [' - ',num2str(xdata(2)),')'];
else
    str6 = [' ',num2str(xdata(2)),')'];
end

% Putting all the strings together.
finalstr = [str1,str2,str3,str4,str5,str6];
disp(finalstr);

% Displaying the final value of the function at xdesired.
disp(sprintf('\nFinal value at xdesired = %g, f(%g) = %g',xdesired,xdesired,fxdesired))

% Displaying the number of significant digits that can be taken and
% absolute relative approximate error.
disp(sprintf('\nAbsolute relative approximate error, abrae = %g%%',ea))
disp(sprintf('\nNumber of significant digits at least correct, sig_dig = %g',sd))

% Creating the inline function for plotting to take place.
func1 = [num2str(b0),'+',num2str(b1),'*(z-',num2str(xdata(1)),')+',num2str(b2)];
func2 = ['*(z-',num2str(xdata(1)),')*(z-',num2str(xdata(2)),')'];
func = inline([func1,func2]);

% Creating the plots for quadratic interpolation.
axis on
figure(2)
subplot(2,1,1)
fplot(func,[min(xdata),max(xdata)])
title('Quadratic Interpolation (Data Points Used)','Fontweight','bold');
xlabel('x data','Fontweight','bold');
ylabel('y data','Fontweight','bold');
hold on
plot(xdata,ydata,'ro','MarkerSize',8','MarkerFaceColor',[1,0,0])
plot(xdesired,fxdesired,'kx','Linewidth',2,'MarkerSize',12')


% Creating the second figure for the Quadratic interpolation.
axis on
subplot(2,1,2)
fplot(func,[min(xdata),max(xdata)])
title('Quadratic Interpolation (Full data set)','Fontweight','bold')
xlabel('x data','Fontweight','bold');
ylabel('y data','Fontweight','bold');
hold on
plot(x,y,'ro','MarkerSize',8','MarkerFaceColor',[1,0,0])
plot(xdesired,fxdesired,'kx','Linewidth',2,'MarkerSize',12')
xlim([min(x) max(x)])
ylim([min(y) max(y)])
hold off
fprev=fnew;


%%%%%%%%% CUBIC INTERPOLATION %%%%%%%%%

% Pick four data points
datapoints=4;
p=1;
for i=1:n
    if d(i) <= datapoints
        xdata(p)=x(i);
        ydata(p)=y(i);
        p=p+1;
    end
end

% Calculating coefficients of Newton's Divided difference polynomial
b0=sym('b0');
b1=sym('b1');
b2=sym('b2');
b3=sym('b3');
z=sym('z');
b0=ydata(1);
b1=(ydata(2)-ydata(1))/(xdata(2)-xdata(1));
b2=((ydata(3)-ydata(2))/(xdata(3)-xdata(2))-(ydata(2)-ydata(1))/(xdata(2)-xdata(1)))/(xdata(3)-xdata(1));
b3=(((ydata(4)-ydata(3))/(xdata(4)-xdata(3))-(ydata(3)-ydata(2))/(xdata(3)-xdata(2)))/(xdata(4)-xdata(2))-b2)/(xdata(4)-xdata(1));
fc=b0 +b1*(z-xdata(1)) + b2*(z-xdata(1))*(z-xdata(2)) + b3*(z-xdata(1))*(z-xdata(2))*(z-xdata(3));
fxdesired=subs(fc,z,xdesired);

fnew=fxdesired;
ea=abs((fnew-fprev)/fnew*100);
if ea >= 5
    sd1=0;
else
    sd1=floor(2-log10(abs(ea)/0.5));
end

% Displaying the results for Cubic Interpolation.
disp('---------------------------------------------------------------------')
disp(sprintf('\nCUBIC INTERPOLATION:'))
disp(sprintf('\nx data chosen: x1 = %g, x2 = %g, x3 = %g, x4 = %g',xdata(1),xdata(2),xdata(3),xdata(4)))
disp(sprintf('\ny data chosen: y1 = %g, y2 = %g, y3 = %g, y4 = %g',ydata(1),ydata(2),ydata(3),ydata(4)))

% Displaying the determined coefficients b0,b1,b2,b3.
disp(sprintf('\nCoefficients of the cubic interpolation: b0 = %g, b1 = %g, b2 = %g, b3 = %g',b0,b1,b2,b3))

% Using concactenation to properly display the final function in the
% command window.
fprintf('\n')
str1 = ['Final Function: f3(x) = ',num2str(b0)]; 
if b1 > 0
    str2 = ['+',num2str(b1),'(x'];
else
    str2 = [' ',num2str(b1),'(x'];
end    

if xdata(1) > 0
    str3 = ['-',num2str(xdata(1)),')'];
    str5 = [str3,'(x'];
else
    str3 = ['+',num2str(xdata(1)),')'];
    str5 = [str3,'(x'];
end

if b2 > 0
    str4 = ['+',num2str(b2),'(x'];
else
    str4 = [' ',num2str(b2),'(x'];
end

if xdata(2) > 0
    str6 = ['-',num2str(xdata(2)),')'];
    str8 = [str5,str6,'(x'];
else
    str6 = [' ',num2str(xdata(2)),')'];
    str8 = [str5,str6,'(x'];
end

if b3 > 0
    str7 = ['+',num2str(b3),'(x'];
else
    str7 = [' ',num2str(b3),'(x'];
end

if xdata(3) > 0
    str9 = ['-',num2str(xdata(3)),')'];
else
    str9 = [' ',num2str(xdata(3)),')'];
end

% Putting it all together.
finalstr = [str1,str2,str3,str4,str5,str6,str7,str8,str9];
disp(finalstr);

% Displaying the value of the function at xdesired.
disp(sprintf('\nFunction value at xdesired = %g, f(%g) = %g',xdesired,xdesired,fxdesired))

% Displaying the number of significant digits that can be taken and
% absolute relative approximate error.
disp(sprintf('\nAbsolute relative approximate error, abrae = %g%%',ea))
disp(sprintf('\nNumber of significant digits at least correct, sig_dig = %g',sd1))

% Creating the inline function
func1 = [num2str(b0),'+',num2str(b1),'*(z-',num2str(xdata(1)),')+',num2str(b2),'*(z-',num2str(xdata(1)),')*(z-'];
func2 = [num2str(xdata(2)),')+',num2str(b3),'*(z-',num2str(xdata(1)),')*(z-',num2str(xdata(2)),')*(z-',num2str(xdata(3)),')'];
func = inline([func1,func2]);

% Creating the plots for the Cubic Interpolation.
axis on
figure(3)
subplot(2,1,1)
fplot(func,[min(xdata),max(xdata)])
title('Cubic Interpolation (Data Points Used)','Fontweight','bold')
xlabel('x data','Fontweight','bold');
ylabel('y data','Fontweight','bold');
hold on
plot(xdata,ydata,'ro','MarkerSize',8','MarkerFaceColor',[1,0,0])
plot(xdesired,fxdesired,'kx','Linewidth',2,'MarkerSize',12')
hold off

% Creating the second plot for Cubic Interpolation
axis on
subplot(2,1,2)
fplot(func,[min(xdata),max(xdata)])
title('Cubic Interpolation (Full Data Set)','Fontweight','bold')
xlabel('x data','Fontweight','bold');
ylabel('y data','Fontweight','bold');
hold on
plot(x,y,'ro','MarkerSize',8','MarkerFaceColor',[1,0,0])
plot(xdesired,fxdesired,'kx','Linewidth',2,'MarkerSize',12')
xlim([min(x) max(x)])
ylim([min(y) max(y)])
hold off

end