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 6.3 (Differentiation Based on N+1 Nodes).
% Section 6.2, Numerical Differentiation Formulas, Page 342
echo on; clc; format long; hold off; clear

% This program differentiates a Newton polynomial approximation.
% 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.

% The example investigates f(x) = cos(x) using 5 nodes.

% The derivative is approximated at X(1).

% The function f(x) is stored in the M-file f.m.
% function z = f(x)
% z = cos(x);

delete f.m
diary f.m; disp('function z = f(x)');...
disp('z = cos(x);');...
diary off;

f(0); % Remark. f.m and diffnew.m are used for Algorithm 6.3

pause % Press any key to continue.

clc;
% Place abscissa for differentiation in x0.

% Place the endpoints interval [a0,b0] about x0 in a0 and b0.

x0 = 0.8;

a0 = 0;

b0 = pi/2;

pause % Press any key to graph the function.

clc; clg;
y0 = f(x0);
hs = (b0-a0)/100;
Xs = a0:hs:b0;
Ys = f(Xs);
a = 0;
b = 1.6;
c = 0;
d = 1.1;
axis([a b c d]);...
plot(x0,y0,'o',Xs,Ys);...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('y = f(x) and the given point');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;
% Place abscissa for differentiation in x0.

% Place the number of points in n.

% Place the step size in h.

x0 = 0.8;
n = 5;
h = 0.01;
X = x0:h:x0+(n-1)*h;
Y = f(X);

[A,df] = diffnew(X,Y);

pause % Press any key to continue.

format short;
Mx0 = 'Approximating the derivative by differentiating a Newton polynomial.';
Mx1 = 'First construct the Newton approximation polynomial.';
Mx2 = 'The abscissas are:';
Mx3 = 'The ordinates are:';
Mx4 ='The coefficients for the Newton polynomial are:';
Mx5 = 'The centers for the Newton polynomial are:';
clc,echo off,diary output,...
disp(''),disp(Mx0),disp(''),disp(Mx1),disp(''),disp(Mx2),disp(X),...
disp(Mx3),disp(Y),disp(''),disp(Mx4),format short,disp(A),...
disp(Mx5),disp(X(1:n-1)),diary off, echo on,format long

pause % Press any key to see the approximate derivative.

clc;
X1 = [a b];
C = [df f(x0)];
Y1 = polyval(C,X1-x0);
axis([a b c d]);...
plot(x0,y0,'or',Xs,Ys,'-g',X1,Y1,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('The tangent line has slope m = f`(x0).');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'Numerical approximation for the derivative.';
Mx2 = ['f`(',num2str(X(1)),')'];
Mx3 = ['P`(',num2str(X(1)),')'];
Mx4 = [Mx2,' ~ ',Mx3,' = '];
clc,disp(''),disp(Mx1),...
disp(''),disp(Mx4),disp(df)

pause % Press any key to zoom in.

a = min(X)-3*h;
b = max(X)+3*h;
c = min(Y)-3*h;
d = max(Y)+3*h;
X1 = [a b];
C = [df f(x0)];
Y1 = polyval(C,X1-x0);
axis([a b c d]);...
X0 = X(2:n);...
Y0 = Y(2:n);...
plot(x0,y0,'or',X0,Y0,'og',Xs,Ys,'-g',X1,Y1,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('The tangent line has slope m = f`(x0).');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx1 = 'Numerical approximation for the derivative.';
Mx2 = ['f`(',num2str(X(1)),')'];
Mx3 = ['P`(',num2str(X(1)),')'];
Mx4 = [Mx2,' ~ ',Mx3,' = '];
clc,echo off,diary output,...
disp(''),disp(Mx1),...
disp(''),disp(Mx4),disp(df),diary off,echo on