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.e (Error Analysis for Numerical Differentiation).
% Section 6.1, Approximating the Derivative, Page 326
echo on; clc; format long; hold off; clear

% This program investigates the error E(f,h) when
% computing numerical approximation for f'(x). Two central
% Two central difference formulas are explored:

% f(x+h) - f(x-h)
% f'(x) = --------------- + E(f,h)
% 2h

% - f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)
% f'(x) = --------------------------------------- + E(f,h)
% 12h

% 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 is used for Algorithm 6.e

pause % Press any key to continue.

clc;
% f(x+h) - f(x-h)
% For f'(x) Å ---------------
% 2h

% The error bounds is known to have the general form:

% |E(f,h)| ² eps/h + m h^2/6

% and the optimum step size is

% h = (3 eps/m)^(1/3)

% where eps is machine epsilon (round off).

pause % Press any key to continue.

clc;
% Place the differentiation point in x

% Place a bound for |f'''(c)| in m

x = 0.8;

m = 1;

hopt = (3*eps/m)^(1/3);

emin = eps/hopt + m*hopt^2/6;

pause % Press any key to graph the error bound.

clc;
a = hopt/10000;
b = 3.44*hopt; %b = 0.00003;
c = 0;
d = 5.25*emin; %d = 0.0000000002;
h = (b-a)/150;
H = a:h:b;
B = eps./H + m*H.^2/6;
E = abs(- sin(x) - (f(x+H) - f(x-H))./(2*H));
axis([a b c d]);...
plot(H,E,'-r',H,B,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('h');...
ylabel('e');...
title('|E(f,h)| ² eps/h + m h^2/6');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;
Mx1 = 'Graph of the error bound B = eps/h + m h^2/6,';
Mx2 = 'and error in the approx. E = |E(f,h)|.';
Mx3 = 'The two curves verify the relationship:';
Mx4 = ' |E(f,h)| ² eps/h + m h^2/6';
Mx5 = 'The optimum step size from the theory is:';
Mx6 = ' hopt = (3 eps/m)^(1/3) = ';
Mx7 = [Mx6,num2str(hopt)];
Mx8 = 'The minimum error from the theory is:';
Mx9 = ' emin = eps/hopt + m hopt^2/6 = ';
Mx10 = [Mx9,num2str(emin)];
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),...
disp(''),disp(Mx3),disp(''),disp(Mx4),...
disp(''),disp(Mx5),disp(''),disp(Mx7),...
disp(''),disp(Mx8),disp(''),disp(Mx10),diary off,echo on

pause % Press any key to continue.

clear; clc;
% - f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)
% for f'(x) Å ---------------------------------------
% 12h

% The error bounds is known to have the general form:

% |E(f,h)| ² 3/2 eps/h + m h^4/30 and the

% Optimum step size

% h = (45/4 eps/m)^(1/5)

% where eps is machine epsilon (round off).

pause % Press any key to continue.

clc;
% Place the differentiation point in x

% Place a bound for |f'''''(c)| in m

x = 0.8;

m = 1;

hopt = ((45/4)*eps/m)^(1/5);

emin = (3/2)*eps/hopt + m*hopt^4/30;

pause % Press any key to graph the error bound.

clc;
a = hopt/10000;
b = 2.51*hopt; %b = 0.003;
c = 0;
d = 5.78*emin; %d = 0.000000000002;
h = (b-a)/150;
H = a:h:b;
B = (3/2)*eps./H + m*H.^4/30;
E = abs(- sin(x)-(-f(x+2*H)+8*f(x+H)- 8*f(x-H)+ f(x-2*H))./(12*H));
axis([a b c d]);...
plot(H,E,'-r',H,B,'-g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('h');...
ylabel('e');...
title('|E(f,h)| ² 3/2 eps/h + m h^4/30');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;
Mx1 = 'Graph of the error bound B = 3/2 eps/h + m h^4/30,';
Mx2 = 'and error in the approx. E = |E(f,h)|.';
Mx3 = 'The two curves verify the relationship:';
Mx4 = ' |E(f,h)| ² 3/2 eps/h + m h^4/30';
Mx5 = 'The optimum step size from the theory is:';
Mx6 = ' hopt = (45/4 eps/m)^(1/5) = ';
Mx7 = [Mx6,num2str(hopt)];
Mx8 = 'The minimum error from the theory is:';
Mx9 = ' emin = 3/2 eps/hopt + m hopt^4/30 = ';
Mx10 = [Mx9,num2str(emin)];
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),...
disp(''),disp(Mx3),disp(''),disp(Mx4),...
disp(''),disp(Mx5),disp(''),disp(Mx7),...
disp(''),disp(Mx8),disp(''),disp(Mx10),diary off,echo on