% CLASS18OCT  Short code to determine the first degree n for
%   which a certain polynomial approximates a certain function
%   within a certain tolerance.
%
%   given:  f(x) = 2^x
%           interval  [0,3]
%           P(x) is equally-spaced poly interpolant of degree n
%
%   seek: first n so that  |f(x) - P(x)| <= 10^-8  for x in [0,3]
%
%   We have already used Lagrange Remainder theorem to get the 
%   estimate
%                         (ln 2)^(n+1) 2^3 3^(n+1)
%     |f(x) - P(x)|  <=  --------------------------
%                                  (n+1)!

n = 1:40;  % enough, it turns out; I used "n=1:100" in class

%% uncomment to show estimate:
%8 * (log(2)*3).^(n+1) ./ factorial(n+1)

%% uncomment to show list which is 0 where inequality is false, and 1 where true:
%find(8 * (log(2)*3).^(n+1) ./ factorial(n+1) < 1.0e-8)

% in class I did this:
f = find(8 * (log(2)*3).^(n+1) ./ factorial(n+1) < 1.0e-8);
f(1)   % = 16

% simpler but equivalent way:
n = find((8 * (log(2)*3).^(n+1) ./ factorial(n+1) < 1.0e-8), 1)

% now CHECK our answer by finding P and computing approximation error

% compute P;  internally polyfit uses the Vandermonde method
x = linspace(0,3,n+1);
y = 2.^x; 
P = polyfit(x,y,n);

% plot f(x) and P(x);  polyval uses result from polyfit to evaluate P(x) [Horner's]
figure(1)
xplot=0:.001:3;  plot(xplot,2.^xplot,xplot,polyval(P,xplot))
grid on,  xlabel x,  legend('f(x)','P(x)')
title('this is not a useful plot because P(x) matches f(x) to screen resolution!')

% compute actual max approx error    max_{0<=x<=3} |f(x) - P(x)|
max(abs(2.^xplot - polyval(P,xplot)))

% plot difference  |f(x) - P(x)|
figure(2)
plot(xplot,abs(2.^xplot - polyval(P,xplot)))
grid on,  xlabel x,  title('difference  |f(x)-P(x)| when n=16')

% redo with much smaller n; seems to suffice!;  this n found by trial-and-error
n = 10;
x = linspace(0,3,n+1);
y = 2.^x; 
P = polyfit(x,y,n);
figure(3)
plot(xplot,abs(2.^xplot - polyval(P,xplot)))
grid on,  xlabel x,  title(['difference  |f(x)-P(x)| when n=' num2str(n)])

