echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 4.5 (Newton Interpolation Polynomial).
% Section 4.4, Newton Polynomials, Page 234
echo on; clc; format short; hold off; clear
% Investigation a Newton polynomial approximation.

% The n+1 are points needed to construct Pn(x).

% The abscissas are stored in X.

% The ordinates are stored in Y.

% The points are counted k=1,2,...,n+1.

pause % Press any key to continue.

% Remark. newpoly.m is used for Algorithm 4.5

clc;
% Newton polynomial approximation Pn(x) of f(x) over [a,b].

% Place the degree of approximation in n.

% Place the left endpoint in a.

% Place the right endpoint in b.

% Place the 'string expression' for f(x) in fun.

n = 4;
a = 0;
b = 4;

fun = 'cos(x)';

pause % Press any key to form the coefficient polynomials.

clc;
% The Newton interpolating polynomial is being constructed.

h = (b-a)/n;

X = a:h:b;

x = X;
Y = eval(fun);

[C,D] = newpoly(X,Y);

pause % Press any key to continue.

format short;
Mx1 = 'Construction of a Newton approximation polynomial.';
Mx2 = 'The abscissas are:';
Mx3 = 'The ordinates are:';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(X),disp(Mx3),disp(Y),...
diary off, echo on

pause % Press any key for the divided difference table.

clc;
Mx = 'The divided difference table is:';
clc,echo off,diary output,...
disp(' '),disp(Mx),disp(D),diary off,echo on

pause % Press any key to graph the Newton approximation.

clc; clg;
a = min(X);
b = max(X);
h = (b-a)/150;
X1 = a:h:b;
Y1 = polyval(C,X1);
x = X1;
Y2 = eval(fun);
plot(X,Y,'or',X1,Y1,'-r',X1,Y2,'-g');...
hold on;...
plot([a-0.2 b+0.2],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['Comparison of ',fun,' and P'];...
Mx2 = [Mx1,num2str(n),'(x).'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx1='The Newton polynomial has been rearranged in ordinary polynomial form.';
Mx2='This ordinary polynomial looks like:';
Mx3='Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
Mx4 = 'The degree is n = ';
Mx5 = ', and the coefficients list C is:';
clc,echo off, diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(''),...
disp(Mx3),disp(''),disp([Mx4,num2str(n-1),Mx5]),disp(''),...
for i=1:5:n+1, disp(C([i:min(i+4,n)])); end,...
diary off, echo on

pause % Press any key to view f(x) - Pn(x).

clc;
Z = zeros(1,length(X));
plot(X,Z,'or',X1,Y2-Y1,'-r');...
hold on;...
plot([a-0.2 b+0.2],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['The error: ',fun,' - P'];...
Mx2 = [Mx1,num2str(n),'(x).'];...
title(Mx2);...
grid;...
hold off;...
shg; pause % Press any key to continue.

pause % Press any key for a list of numerical computations.

clc; format long;
X = a:0.25:b;
x = X;
Y = eval(fun);
P = polyval(C,X);
points = [X;Y;P;Y-P];
Mx1=['Newton polynomial approximation of f(x) = ',fun];
Mx2=' x(k) f(x(k)) Pn(x(k)) error';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
diary off,echo on