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 5.3 (Non-linear Curve Fitting).
% A detailed look at the exponential curve fit.
% Section 5.2, Curve Fitting, Page 280
echo on; clc; format short; hold off; clear

% This program performs exponential curve fitting,

% by using the method of data linearization. Given a set

% of data points { (x , y ), (x , y ) ,..., (x , y ) },
% 1 1 2 2 n n

% The abscissas and ordinates are stored in X and Y, respectively.

% X = [x , x ,..., x ]; Y = [y , y ,..., y ];
% 1 2 n 1 2 n

pause % Press any key to continue.

clc;
% Place the abscissas for the points in X.

% Place the ordinates for the points in Y.

X = [0 1 2 3 4];

Y = [1.5 2.5 3.5 5.0 7.5];

pause % Press any key to graph data points.

clc; clg;
axis([-0.1 4.1 -0.2 8.2]);
plot(X,Y,'or');...
hold on;...
plot([-10000 10000],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
title('The given x-y data points');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

points = [X;Y];
format short;
clc,echo off,diary output,disp('The given x-y points:'),...
disp(' x y'),disp(points'),diary off,echo on
pause % Press any key to "linearize" the data points.

clc; clg;
XT = X;
YT = log(Y);
Ya = 1;
p = polyfit(XT,YT,1);
A = p(1);
B = p(2);
Xs = [min(X) max(X)];
Ys = polyval(p,Xs);
axis([-0.2 4.2 -0.1 2.1]);...
plot(XT,YT,'or',Xs,Ys,'-g');...
hold on;...
plot([-10000 10000],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('X');...
ylabel('Y');...
Mx1 = 'Least squares line: Y = ';...
Mx2 = [Mx1,num2str(A),' X'];...
if B > 0,
Mx2 = [Mx2,' + ',num2str(B)];
else
Mx2 = [Mx2,' - ',num2str(abs(B))];
end;...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

points = [XT;YT];
format short;

Mx1 = 'The transformed points in the X-Y plane.';
Mx3 = ' X Y';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(Mx2),...
disp(Mx3),disp(points'),diary off, echo on
pause % Press any key to fit y = C exp(Ax).

clc; clg;
A = p(1);
B = p(2);
C = exp(B);
h = (max(X)-min(X))/150;
Xs = min(X):h:max(X);
Ys = C*exp(A*Xs);
axis([-0.1 4.1 -0.2 8.2]);...
plot(X,Y,'or',Xs,Ys,'-g');...
hold on;...
plot([-10000 10000],[0 0],'b',[0 0],[-10000 10000],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = 'The fit f(x) = ';...
Mx2 = [Mx1,num2str(C),' exp(',num2str(A),' x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

points = [X;Y;C*exp(A*X);Y-C.*exp(A*X)]';
Mx3=' x(k) y(k) f(x(k)) error';
clc,echo off,diary output,disp(''),disp(Mx2),...
disp(Mx3),disp(points),diary off,echo on