% October 2017
%
% Written by Per-Gunnar Martinsson, University of Oxford
%
% The methods follow Lloyd N. Trefethen's 
%
%    "Approximation Theory and Approximation Practice"
%
% textbook published by SIAM in 2013. Many code snippets are taken from 
% this text. They have been modified freey, however, and any errors are 
% due to PGM. 

function lecture01

DRIVER_uniform
DRIVER_chebpoints
DRIVER_cheb

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_uniform

%rng(0)
n = 5;
f = inline('exp(-x).*sin(3.3*x)./(1+x.*x)','x');
%f = inline('abs(x).*exp(-x).*sin(3*x)./(1+x.*x)','x');
%f = inline('exp(-x).*sin(3*x)./(1+x.*x) + sign(x-0.5)','x');

xx = linspace(-1,1,n+1)';
%xx = sort(1 - 2*rand(1,n+1))';
ff = f(xx);

XX = xx*ones(1,n+1);
NN = ones(n+1,1)*(0:n);
V  = XX.^NN;
aa = V\ff;

figure(1)
yy = linspace(-1,1,1000);
pp = aa(1)*ones(size(yy));
for i = 1:n
  pp = pp + aa(i+1)*(yy.^i);
end

plot(yy,f(yy),'k',...
     yy,pp,'g',...
     xx,f(xx),'rx',...
     'LineWidth',2)
legend('Given function','Approximating polynomial','Interp. points')

ss = svd(V);
fprintf(1,'Conditioning of Vandermode matrix = %10.3e\n',ss(1)/ss(end))

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_chebpoints

%rng(0)
n = 5;
%f = inline('exp(-x).*sin(3*x)./(1+x.*x)','x');
f = inline('abs(x).*exp(-x).*sin(3*x)./(1+x.*x)','x');
%f = inline('exp(-x).*sin(3*x)./(1+x.*x) + sign(x-0.5)','x');

tt = linspace(0,pi,n+1)';
xx = cos(tt);
ff = f(xx);

XX = xx*ones(1,n+1);
NN = ones(n+1,1)*(0:n);
V  = XX.^NN;
aa = V\ff;

figure(1)
yy = linspace(-1,1,1000);
pp = aa(1)*ones(size(yy));
for i = 1:n
  pp = pp + aa(i+1)*(yy.^i);
end

plot(yy,f(yy),'k',...
     yy,pp,'g',...
     xx,f(xx),'rx',...
     'LineWidth',2)
legend('Given function','Approximating polynomial','Interp. points')

ss = svd(V);
fprintf(1,'Conditioning of Vandermode matrix = %10.3e\n',ss(1)/ss(end))

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_cheb

%rng(0)
n = 5;
%f = inline('abs(x).*exp(-x).*sin(3*x)./(1+x.*x)','x');
%f = inline('exp(-x).*sin(3*x)./(1+x.*x)','x');
f = inline('exp(-x).*sin(3*x)./(1+x.*x) + sign(x-0.5)','x');

tt = linspace(0,pi,n+1)';
xx = cos(tt);
ff = f(xx);

V = zeros(n+1);
for i = 0:n
  T = chebpoly(i);
  V(:,i+1) = T(xx);
end

aa = V\ff;

figure(1)
yy = linspace(-1,1,1000);
T = chebpoly(0);
pp = aa(1)*T(yy);
for i = 1:n
  T = chebpoly(i);
  pp = pp + aa(i+1)*T(yy);
end

plot(yy,f(yy),'k',...
     yy,pp,'g',...
     xx,f(xx),'rx',...
     'LineWidth',2)
legend('Given function','Approximating polynomial','Interp. points')

ss = svd(V);
fprintf(1,'Conditioning of Vandermode matrix = %10.3e\n',ss(1)/ss(end))

return