% Script File: ExpTaylor2.m
% Plots, as a function of n, the relative error in the 
% Taylor approximation
%
%           1 + x + x^2/2! +...+ x^n/n! 
%
% to exp(x), together with running bounds on the error,
% for x=10 and x=-10, running both backwards and forwards.

N = 51;  % zero to 50

% Pre-compute the terms of the Taylor series, along with
% a bound for the absolute error in each term.
%
% For exp(-10), we'll just negate the odd-order terms
% (and keep the same error bound); two more columns
% will be the same thing in the reverse order.

term = zeros(N,4);  % [10fwd, -10fwd, 10rev, -10rev]
termErr = zeros(N,4);
DELTA = 0;

% Easiest to do -10fwd first and then take absolute values

thisterm = 1;
for k=1:N
  term(k,2) = thisterm;
  termErr(k,2) = abs(thisterm*DELTA); % 0 term is exact

  thisterm = - thisterm * 10;
    DELTA = DELTA + eps;    %  = 2*eps/2
  thisterm = thisterm / k;
    DELTA = DELTA + eps;
end

% Fill in terms for 10fwd
% term(:, 1) means the entire 1st column

term(:, 1) = abs(term(:, 2));
termErr(:, 1) = termErr(:, 2);

% Reverse order

term(:, 3:4) = term(N:-1:1, 1:2);
termErr(:, 3:4) = termErr(N:-1:1, 1:2);

s = zeros(N,4);     % four sets of partial sums
sumErr = zeros(N,4);  % running absolute error bounds
actErr = zeros(N,4);  % actual relative error
                      % from partial sum and roundoff

s(1, :) = term(1, :);
f = exp([10 -10 10 -10]);

for k=2:N  % 1st row already taken care of
  s(k, :) = s(k-1, :) + term(k, :);
    sumErr(k,:) = (sumErr(k-1, 1:4) + termErr(k, 1:4)) * (1+eps) ...
                  + (abs(s(k, :)) + abs(term(k, :))) * eps/2;
  actErr(k, :) = abs(s(k, :) -f)./f;
end

% Convert absolute error bounds to relative error bounds

sumErr = sumErr./ (abs(s) - sumErr);
% This is assuming abs(sum) > sumErr, which might not
% always be true.  (For example, in exp(-1), the partial
% sum 1 -1 is zero.)

titles = [' exp(10), forward '
          ' exp(-10), forward'
          ' exp(10), backward'
          'exp(-10), backward'];  % matrix of characters

for k=1:4
  subplot(2,2,k);
  
  % The 'actual error' includes the error from only summing
  % a few terms of the Taylor series, which is much larger
  % than the running error bounds on the accuracy of the
  % partial sum until the end of the summation.
  
  % In order to plot them on the same scale, here we only
  % plot the last 10 terms of the 'actual error'.
  
  semilogy(2:N,sumErr(2:N,k) ,N-9:N,actErr(N-9:N,k))
  ylabel('Relative Error in Partial Sum.')
  title(titles(k, :))        
end