%  scripts and functions for MAT 331, Tuesday, Sept 19, 2017


%-------------------------------------------------
% script   compare a function  to it  Taylor series 
N=20;
f=@(x) exp(x); a=-2; b=2; 
series=1./factorial(0:N);

%f =@(x) sin(x); a=-2*pi; b=2*pi;
%series= imag(i.^[0:N])./factorial(0:N);

%f = @(x)  1./(1-x); a=-1;b=1;
%series = ones(1,20);

figure;
hold on;
grid on;
title('function and taylor approximations');
x=[a:.01:b];
for k=1:10
    plot(x,polyval(series(k:-1:1),x),'LineWidth',2);
end
plot(x,f(x),'k','LineWidth',3);


%-------------------------------------------------


% function to compute polynomial interpolation
% input is a vector x of distinct real numbers
% and a equal length vector y.  Finds polynomial
% so that p(x_k) = y_k

function  p = poly_inter(x,y)
   L =length(x);
   t=[1:L]
   p=zeros(1,L);     % initialize output to zero
   for k=1:L
      s=t(t~=k);     % remove kth index
      L_k=poly(x(s)); % make polynomial with roots at x_j,  j~=k
      L_k = L_k / polyval(L_k,x(k));  % normalize so L_k(x(k))=1
      p = p + y(k) *L_k;
   end % for k
return;
   figure; hold on; grid on;
   title('Polynomial interpolation');
   scatter(x,y,'r');
   a=min(x);
   b=max(x);
   delta=(b-a)/100;
   t=[a:delta:b];
   plot(t,polyval(p,t),'b');
return;

%-------------------------------------------------

% function to compute polynomial interpolation
% of a given function at evenly  points
% so that p(x_k) = y_k
% try things like exp, sin, abs(x), 1/(1+x^2)

function  p = even_inter(f,a,b,n)
   delta=(b-a)/n;
   x= [a:delta:b];
   poly_inter(x,f(x));
return;


%-------------------------------------------------

% function to compute polynomial interpolation
% of a given function at the Chebyshev points

function  p = cheb_inter(f,a,b,n)
   x=cos(pi * [0:n]/n)
   x= a+ (b-a)*(x+1)/2;
   poly_inter(x,f(x));
return;


%-------------------------------------------------
%   function to use polynomial interpolation to integrate
%   interpolation points are evenly spaced

function out =  even_integrate(f,a,b,n)
   % x=cos(pi * [0:n]/n);  % define Chebushev points on [-1,1]
   % x= a+ (b-a)*(x+1)/2;  % rescale x to fit [a,b]
   delta =(b-a)/n;
   x=[a:delta:b];
   L=length(x);
   t=[1:L];
   y=f(x);
   for k=1:L
      s=find(t~=k);     % remove kth index from t
      L_k=poly(x(s)); % make polynomial with roots at x_j,  j~=k
      L_k = L_k / polyval(L_k,x(k));  % normalize so L_k(x(k))=1
      Q_k=polyint(L_k);
      w(k) = diff(polyval(Q_k,[a,b]));
   end % for k
   out = sum(w.*f(x));
   % now plot the functions (not needed for integral)
   p=poly_inter(x,y)
   figure; hold on; grid;
   title('Plots of function and polynomial with interpolation  points');
   t=[a:(b-a)/100:b];
   plot(t,polyval(p,t),'r','LineWidth',3); % plot poly
   plot(t,f(t),'--k','LineWidth',2);       % plot function
   scatter(x,polyval(p,x),'k','filled')    % plot interpolation poi9nts
return

%-------------------------------------------------

%   plot function  and its polynomial interpolation 

function out =  plot_inter(f,x)
   a=min(x);
   b=max(x);
   delta=(b-a)/100;
   p=poly_inter(x,f(x))
   figure; hold on; grid;
   title('Plots of function and polynomial at Chebyshev points');
   t=[a:(b-a)/100:b];
   plot(t,polyval(p,t),'r','LineWidth',3); % plot poly
   plot(t,f(t),'--k','LineWidth',2);       % plot function
   scatter(x,y,'k','filled')    % plot interpolation poi9nts
return


%-------------------------------------------------
%   function to use polynomial interpolation to integrate
%   interpolation points are Chebyshev points 

function out =  cheb_integrate(f,a,b,n)
   x=cos(pi * [0:n]/n);  % define Chebushev points on [-1,1]
   x= a+ (b-a)*(x+1)/2;  % rescale x to fit [a,b]
   % delta =(b-a)/n;     % used for evenly spaced points
   % x=[a:delta:b];      
   L=length(x);
   t=[1:L];
   y=f(x);
   for k=1:L
      s=find(t~=k);     % remove kth index from t
      L_k=poly(x(s)); % make polynomial with roots at x_j,  j~=k
      L_k = L_k / polyval(L_k,x(k));  % normalize so L_k(x(k))=1
      Q_k=polyint(L_k);
      w(k) = diff(polyval(Q_k,[a,b]));
   end % for k
   out = sum(w.*f(x));
   % plot function and polynomials (not needed)
   plot_inter(f,x);
return


%-------------------------------------------------




%%
% Chris Bishop
% Homework 1, MAT 331


%%
% Set initial variables 
x=1;
N=20;

%% 
% Compute Formula (1) in homework 

for k=1:N
    formula2(k)=(1+x/k)^k;
end

%%
% Compute Formula 2: 

for k=1:N
    y(k)=x^k/factorial(k);
end
formula1=1+cumsum(y);

%%
% Nest we plot both fomulas to show how they both approach e. Note that 
% Formula 1 appoaches the limit faster

figure;
hold on;
plot(formula1)
plot(formula2)

figure;
hold on;
plot(abs(formula1 - exp(x)))
plot(abs(formula2 - exp(x)))
legend('Formula 1', 'Formula 2');