%  DEMOFL  demonstrates some phenomena related
%  to floating point arithmetic.
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Feb 12, 2001
%
clear; clc;
disp('                                                                  ')


disp(' Hit a key to continue '); pause

clear;clc;
disp('  Example 1                                     ')
disp('  Computation of  1-cos(x)  for x near 0        ')
disp('                                                ')


x = 1;
fprintf(1, '         x             1-cos(x) \n' );
for i = 1:20
    fprintf(1, '    %12.6e     %12.6e \n', x, 1-cos(x)  );
    x = x/4;
end


disp(' Hit a key to continue '); pause

clear;clc;


x = 1;
fprintf(1, '         x           1-cos(x)      sin(x)^2/      Taylor \n' );
fprintf(1, '                                  (1+cos(x))   \n' );
for i = 1:20
    res1 = 1-cos(x);
    res2 = sin(x)^2/(1+cos(x));
    res3 = ((x^2/360 - 1/12)*x^2+1)*x^2/2;
    fprintf(1, '    %12.6e   %12.6e   %12.6e   %12.6e \n', x, res1, res2, res3 );
    x = x/4;
end


disp(' Hit a key to continue '); pause
clear;clc;

disp('  Example 2                                         ')
disp('  Solution of a*x^2 + b*x + c = 0                   ')         
disp(' ')

a = 5.e-6;
b = 100;
c = 5.e-7;

fprintf(1, '   a = %12.6e   \n', a );
fprintf(1, '   b = %12.6e   \n', b );
fprintf(1, '   c = %12.6e   \n\n', c );


x1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a);
x2 = (-b - sqrt(b^2 - 4*a*c)) / (2*a);
fprintf(1, ' x+ = (-b + sqrt(b^2 - 4*a*c)) / (2*a) = %12.6e   \n', x1 );
fprintf(1, ' x- = (-b - sqrt(b^2 - 4*a*c)) / (2*a) = %12.6e   \n', x2 );

x1 = 2*c / (-b - sqrt(b^2 - 4*a*c));
fprintf(1, ' x+ = 2*c / (-b - sqrt(b^2 - 4*a*c))   = %12.6e   \n', x1 );





disp(' Hit a key to continue '); pause
clear;clc;

disp('  Example 3                                         ')
disp('  Evaluation of  exp(x)  using Taylor series        ')
disp('                                                    ')

x_1    = -10;
fprintf(1, '   x = %12.6e   \n', x_1 );

expT_1(1) = 1;
sum_1  = 1;
expT_2(1) = 1;
sum_2  = 1;
x_2    = -x_1;

N = 50;
fprintf(1, '   N    sum_k=1^N x^k/k!  sum_k=1^N (-x)^k/k!  \n' );
for k = 1:N
    sum_1  = sum_1 * x_1 / k;
    expT_1(k+1) = expT_1(k) + sum_1;
    sum_2  = sum_2 * x_2 / k;
    expT_2(k+1) = expT_2(k) + sum_2;
    fprintf(1, ' %3d    %14.6e   %14.6e    \n', k, expT_1(k+1), 1/expT_2(k+1) );
end


figure(1)
semilogy((0:N), expT_1, 'r', (0:N), expT_2.^(-1), 'g' )
title('Taylor approximation of exp(-10)')
xlabel('N')
legend('sum_{k=1}^N x^k/k! ', '1/ sum_{k=1}^N (-x)^k/k!')

figure(2)
semilogy((0:N), abs(expT_1'-exp(x_1)*ones(N+1,1))/exp(x_1), ...
         (0:N), abs((expT_2.^(-1))'-exp(x_1)*ones(N+1,1))/exp(x_1) )
title('Relative error in Taylor approximation of exp(-10)')
xlabel('N')
legend('|e^x - sum_{k=1}^N x^k/k!|/|e^x|', '|e^x 1/ sum_{k=1}^N (-x)^k/k!|/|e^x|')




disp(' Hit a key to continue '); pause
clear;clc;

disp('  Example 4                                         ')
disp('  One-sided finite difference approximation of the  ')
disp('  derivative of exp(x) at x = 1                     ')
disp('                                                    ')


for i = 1:20
    dx(i) = 10^(-i);
    df(i) = (exp(1+dx(i))-exp(1))/dx(i);
end


loglog(dx,abs(df - exp(1)))
xlabel('\delta')
ylabel('|exp(1) - (exp(1+\delta)-exp(1))/\delta|')