%% MATLAB Precision Demo 
% The purpose of this file is to demonstrate basic limitations of 
% floating point arithmetic. 

%% Set up Domain

% Grid Spacing:
dx=pi/100;
dy=pi/100;
% Grid Vectors:
x=[0:dx:2*pi];
y=[0:dy:2*pi];
% 2D Meshgrid Arrays:
[X,Y]=meshgrid(x,y);


%% Demonstrate Precision (small numbers)

a=realmin('single').*eps('single');

% Demonstrate Precision Limitation:
figure(1);
for m=1:20
    n=21-m;
    %Plot single-precision solution
    plot(x,n.*a.*sin(x));
    hold on;
    %Plot double-precision solution
    plot(x,double(n.*a).*sin(x));
    %Plot absolute error
    spvalue=n.*a.*sin(x);
    plot(x,double(n.*a).*sin(x)-double(spvalue));
    hold off;
    axis([x(1) x(end) -n*a n*a]);
    xlabel('x (rad.)');
    ylabel('y (arb.)');
    title(['Single Precision ',int2str(n),'*realmin*eps*sin(x)']);
    pause(0.5);
end;
    
% Demonstrate Precision Limitation:
figure(2);
plot(x,a.*sin(x),x,2.*a.*sin(x),x,3.4.*a.*sin(x));
hold on;
plot(x,double(a).*sin(x),x,double(2.*a).*sin(x),x,double(3.*a).*sin(x));
hold off;
xlabel('x (rad.)');
ylabel('y (arb.)');
title('Single Precision Sinusoids, Near Machine Precision');
axis([x(1) x(end) -4*a 4*a]);

%% Demonstrate Precision (deviations from 1)

a=eps('single');

% Demonstrate Precision Limitation (about 1):
figure(3);
for m=1:20
    n=21-m;
    %Plot single-precision solution
    plot(x,1+n.*a.*sin(x));
    hold on;
    %Plot double-precision solution
    plot(x,1+double(n.*a).*sin(x));
    %Plot absolute error + 1
    spvalue=(1+n.*a.*sin(x));
    plot(x,1+double(n.*a).*sin(x)-double(spvalue)+1);
    hold off;
    axis([x(1) x(end) (1-n*a) (1+n*a)]);
    xlabel('x (rad.)');
    ylabel('y (arb.)');
    title(['Single Precision 1+',int2str(n),'*eps*sin(x)']);
    pause(0.5);
end;

% Can you explain the asymmetry? 

%% Summary
% Here we have demonstrated some limits of precision.
%
% (J. B. Snively, for EP501, 9/1/2016.)