% 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 lecture06

DRIVER_analytic
DRIVER_gibbs_interp
DRIVER_gibbs_proj
DRIVER_gibbs_general
DRIVER_remez1
DRIVER_remez2

return

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

function DRIVER_analytic

x    = chebfun('x'); 
beta = 1/20;
f    = 1/(1 + x*x/(beta*beta));
nn   = 10:10:200;

err = 0*nn;

for j = 1:length(nn)
  n      = nn(j); 
  fn     = chebfun(f,n+1); 
  err(j) = norm(f-fn,inf);
end

figure(1)
hold off
nmax = max(nn);
fn = chebfun(f,2*nmax);
kk = 0:nmax;
aa = chebcoeffs(fn);
loglog(kk,abs(aa(1:(nmax+1))),'b.',....
       'MarkerSize',30,...
       'LineWidth',3)
title('Magnitude of Chebyshev coefficient: f(x) = 1/(1 + 25*x*x)')
xlabel('k')
ylabel('|a_k|')

figure(2)
hold off
loglog(nn,err,'b.',....
       'MarkerSize',35,...
       'LineWidth',3)
title('Errors of Chebyshev interpolants in sup norm: f(x) = 1/(1 + 25*x*x)')
xlabel('n')
ylabel('||f - f_n||')
keyboard

figure(3)
hold off
semilogy(nn,err,'b.',....
         nn,(1/(1+beta)).^nn,'g--',...
       'MarkerSize',35,...
       'LineWidth',3)
title('Errors of Chebyshev interpolants in sup norm: f(x) = 1/(1 + 25*x*x)')
xlabel('n')
ylabel('||f - f_n||')
keyboard

return

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

function DRIVER_gibbs_interp

x   = chebfun('x');
f   = sign(x);
n   = 50;
fn  = chebfun(f,n+1);
xx  = chebpts(n+1);
xxx = linspace(-1,1,1000);

plot(xxx,f(xxx),'b',...
     xxx,fn(xxx),'r',...
     xx,fn(xx),'g.',...
     'MarkerSize',20)
keyboard

nn = 10*(2.^(1:10))+1;
for i = 1:length(nn)
  n   = nn(i);
  fn  = chebfun(f,n);
  xx  = chebpts(n);
  xxx = linspace(-1,1,1000);
  plot(xxx,f(xxx),'b',...
       xxx,fn(xxx),'r',...
       'MarkerSize',20)
  fprintf(1,'n = %4d   max|p_n| = %9.7f\n',n,max(fn))
  keyboard
end

return

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

function DRIVER_gibbs_proj

x   = chebfun('x');
f   = sign(x);
n   = 21;
fn  = chebfun(f,'trunc',n+1);
xx  = chebpts(n+1);
xxx = linspace(-1,1,1000);

plot(xxx,f(xxx),'b',...
     xxx,fn(xxx),'r',...
     xx,f(xx),'g.',...
     'MarkerSize',20)
keyboard

nn = 10*(2.^(1:10))+1;
for i = 1:length(nn)
  n   = nn(i);
  fn  = chebfun(f,'trunc',n+1);
  xx  = chebpts(n+1);
  xxx = linspace(-1,1,1000);
  plot(xxx,f(xxx),'b',...
       xxx,fn(xxx),'r',...
       'MarkerSize',20)
  fprintf(1,'n = %4d   max|p_n| = %9.7f\n',n,max(fn))
  keyboard
end

return

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

function DRIVER_gibbs_general

x   = chebfun('x');
f   = cos(3*x) - 0.7*sign(x-0.3) + 2.0*abs(x-0.3);
n   = 101;
fn  = chebfun(f,n+1);
xx  = chebpts(n+1);
xxx = linspace(-1,1,1000);

plot(xxx,f(xxx),'b',...
     xxx,fn(xxx),'r',...
     xx,fn(xx),'g.')
keyboard

nn = 10*(2.^(1:10))+1;
for i = 1:length(nn)
  n   = nn(i);
  fn  = chebfun(f,n+1);
  xx  = chebpts(n+1);
  xxx = linspace(-1,1,1000);
  plot(xxx,f(xxx),'b',...
       xxx,fn(xxx),'r',...
       'MarkerSize',20)
  fprintf(1,'n = %4d   max|p_n| = %9.7f\n',n,max(fn))
  keyboard
end

return

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

function DRIVER_remez1

x     = chebfun('x');
f     = min(abs(x+0.5),2*abs(x-0.5));
n     = 6;
pstar = minimax(f,n);
err   = max(abs(f-pstar));

xx = linspace(-1,1,1000);
figure(1)
subplot(1,2,1)
hold off
plot(xx,f(xx),'b',...
     xx,pstar(xx),'r',...
     'LineWidth',2)
legend('f','p^{*}','Location','NorthWest')
title(sprintf('Best approximant for n=%d',n))
subplot(1,2,2)
hold off
plot(xx,f(xx)-pstar(xx),'r',...
     xx, err*ones(size(xx)),'g',...
     xx,-err*ones(size(xx)),'g',...
     'LineWidth',2)
title(sprintf('Error "f - p^{*}r" for n=%d',n))
keyboard

return

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

function DRIVER_remez2

x     = chebfun('x');
f     = sin(exp(3*x));
n     = 28;
pstar = minimax(f,n);
err   = max(abs(f-pstar));

xx = linspace(-1,1,1000);
figure(1)
subplot(1,2,1)
hold off
plot(xx,f(xx),'b',...
     xx,pstar(xx),'r',...
     'LineWidth',2)
legend('f','p^{*}','Location','NorthWest')
title(sprintf('Best approximant for n=%d',n))
subplot(1,2,2)
hold off
plot(xx,f(xx)-pstar(xx),'r',...
     xx, err*ones(size(xx)),'g',...
     xx,-err*ones(size(xx)),'g',...
     'LineWidth',2)
title(sprintf('Error "f - p^{*}" for n=%d',n))
keyboard

return

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