%% MATLAB Approximation and Errors Demo 
% The purpose of this file is to demonstrate basic limitations of 
% numerical approximations, including data, rounding, truncation, 
% and discretization errors. 

%% Set up Domain

% Grid Spacing:
dx=1/100;
% Grid Vectors:
x=[-1:dx:1];

%% Demonstrate Data Approximation Errors 

figure(1);
for m=1:10
    %Pie=pi to m significant figures
    pie=round(pi,m);
    accurate=sin(pi.*x./2);
    approximate=sin(pie.*x./2);
    %Plot sin(pi*x) solution
    plot(x,accurate);
    hold on;
    %Plot sin(pie*x) solution
    plot(x,approximate);
    %Plot relative error %
    plot(x,100.*abs((approximate-accurate)./accurate));
    hold off;
    legend('Accurate','Approximate','% Error','Location','Northwest');
    axis([x(1) x(end) -1 1]);
    xlabel('x (arb.)');
    ylabel('y (arb.)');
    title(['Pi to ',int2str(m),' significant figures']);
    pause(1);
end;

%% Demonstrate Precision and Rounding Errors

a=eps('single');

% Demonstrate Precision Limitation (about 1):
figure(2);
for m=1:20
    n=m;
    %Evaluate small oscillation about "1" in single and double precision
    spvalue=1+n.*a.*sin(pi.*x./2);
    dpvalue=1+double(n.*a).*sin(pi.*x./2);
    %Plot single-precision solution
    plot(x,spvalue);
    hold on;
    %Plot double-precision solution
    plot(x,dpvalue);
    %Plot absolute error + 1
    plot(x,dpvalue-double(spvalue)+1);
    hold off;
    legend('Single','Double','Error','Location','Northwest');
    axis([x(1) x(end) (1-n*a) (1+n*a)]);
    xlabel('x (arb.)');
    ylabel('y (arb.)');
    title(['Single Precision 1+',int2str(n),'*eps*sin(pi.*x./2)']);
    pause(0.5);
end;

%% Demonstrate Truncation Errors 

figure(3);
for m=0:5
    accurate=sin(pi.*x./2);
    %First term in Taylor Series
    approximate=0;
    for n=0:m
       nextterm=((-1).^n).*((x.*pi./2).^(2.*n+1))./factorial(2.*n+1);
       approximate=approximate+nextterm;
    end;
    %Plot sine solution
    plot(x,accurate);
    hold on;
    %Plot series solution
    plot(x,approximate);
    %Plot relative error
    plot(x,abs((approximate-accurate)./accurate));
    hold off;
    legend('Accurate','Approximate','Relative Error','Location','Northwest');
    axis([x(1) x(end) -1 1]);
    xlabel('x (arb.)');
    ylabel('y (arb.)');
    title(['Sine Approximated to ',int2str(m+1),' terms in Taylor Series']);
    pause;
end;

%% Demonstrate Discretization Errors 

figure(3);
for m=1:100
    % Grid Spacing:
    dxA=1/(m);
    % Grid Vectors:
    xA=[-1:dxA:1];
    %Sine on a higher-res grid
    accurate=sin(pi.*x./2);
    %Sine on a coarse grid
    approximate=sin(pi.*xA./2);
    %Nearest-neighbor interpolation to higher grid resolution
    approximateI=interp1(xA,approximate,x,'nearest');
    %Plot sine solution
    plot(x,accurate);
    hold on;
    %Plot low-res sine solution
    plot(x,approximateI);
    %Plot relative error
    plot(x,abs((approximateI-accurate)./accurate));
    hold off;
    legend('Accurate','Approximate','Relative Error','Location','Northwest');
    axis([x(1) x(end) -1 1]);
    xlabel('x (arb.)');
    ylabel('y (arb.)');
    title(['Sine Evaluated on ',int2str(m.*2+1),' grid points']);
    pause(0.25);
end;


%% Summary
% Here we have demonstrated some errors.
%
% (J. B. Snively, for EP501, 9/6/2016.)