% template.m
% This is a template for electronic submission of homework
% assignments for the course MATH 3Q03; use of this template is 
% obligatory
% Email the complete file to the instructor before the
% due date; before submission make sure that the code executes
% correctly in MATLAB

% Written by Bartosz Protas, September 2016

function hw3_sol

close all, clear all

% Student information
  
disp('  ------------   Student Information  ------------  ')  

% Please enter your information here
Name = 'Johny';
Surname = 'Good';
ID = 123456789;

fprintf('  Student:  %s %s (ID: %d) \n', Name, Surname, ID);



disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')


disp('  ----------   Solution to Question #1  ----------  ')  
disp('  ')
% Insert here your solution to Question #1
% (if you define your own functions, they may appear 
%  at the end of the file)


% part (a) - plotting the function;
x = [-5:0.1:5];
y = [-5:0.1:5];
[X,Y] = meshgrid( x, y );
Z = F(X,Y);

figure(1); hold on
contourf(X,Y,Z);
axis equal
colorbar
xlabel('x'), ylabel('y')
%surf(X,Y,Z);

tol = 1.0e-6; total = 100; 

% the two initial guesses;
x0(:,1) = [2;1];
x0(:,2) = [5;3];


% part (b) - finding the critical points and determining their type;
for i=1:2
  k = 0; epsilon = 1;
  x = x0(:,i); % initial guess 
  clear J;
  J = [];
  while ((epsilon > tol) & (k < total)) % Newton method    
  
    if (i == 1)
      plot(x(1),x(2),'or','MarkerFaceColor','r');
    else
      plot(x(1),x(2),'oy','MarkerFaceColor','y');
    end;

    
    xnew = x - H(x(1),x(2)) \ dF( x(1), x(2) ); 
    epsilon = norm(xnew-x);
    x=xnew; 
    k = k+1;
    J = [ J; [F(x(1), x(2)), norm( dF(x(1), x(2)), 2)] ];
    
    
  end
  fprintf('Critical point x = %12.10f, y=%12.10f found after %2.0f iterations \n',...
	  x(1),x(2),k);
  
  disp('Checking for definiteness of the Hessian')
  E = eig(H(x(1),x(2)))
  if ( E(1) >= 0.0 & E(2) >= 0.0 )
    fprintf('The critical point is a local maximum');
  elseif ( E(1) <= 0.0 & E(2) <= 0.0 )
    fprintf('The critical point is a local minimum');
  else
    fprintf('The cirtical point is a saddle');
  end;
  fprintf('\n\n');
  

end;


disp('  ')
disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')


disp('  ----------   Solution to Question #2  ----------  ')  
disp('  ')
% Insert here your solution to Question #2
% (if you define your own functions, they may appear 
%  at the end of the file)


disp('  ')
disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')



xc = chebfun(@(x) x,[0.0,2*pi]);
fc = chebfun(@(x) sin(x),[0.0,2*pi]);

% maximum number of terms
Nc = 10;

figure(2), hold on
plot(fc,'-b','LineWidth',4)
axis([0.0 2*pi -2 2]);

wc = 0.0 * xc;
for n=0:Nc-1
  wc = wc + (-1)^n * xc.^(2*n+1) / factorial(2*n+1);
  plot(wc,'-r');
  
  Err(n+1) = norm(fc-wc,inf);
  
end;
xlim([0 2*pi]), ylim([-2 2])
title('sin(x) and its polynomial approximation')

figure(3)
semilogy(1:Nc,Err,'-x')
xlabel('N'), ylabel('Error')
title('Error of polynomial interpolation');


disp('  ----------   Solution to Question #3  ----------  ')  
disp('  ')
% Insert here your solution to Question #3
% (if you define your own functions, they may appear 
%  at the end of the file)

% part (a) - condition number of the Vandermonde interpolation
% matrix;

E = 1.0e6;
kappa = 1.0;
n = 2;
while ( kappa <= E )
  Xn = linspace(-pi,pi,n+1);  
  kappa  = cond( interp_matrix( Xn, n ) );
  fprintf('  n=%2d    kappa=%e \n', n, kappa );
  n = n + 1;
end;
fprintf('  The condition number of the Vandermonde matrix exceeds 10^6 when n=%d \n\n\n', n-1)


% part (b) - plotting Lagrange polynomials;
N = 4;
X0 = linspace(-pi,pi,101);
X4 = linspace(-pi,pi,N+1);
figure(4); plot ( X0, Lag( X0, X4, 4, 1), 'r-',...
		  X0, Lag( X0, X4, 4, 2), 'g-',...
		  X0, Lag( X0, X4, 4, 3), 'm-',...
		  X0, Lag( X0, X4, 4, 4), 'b-',...
		  X0, Lag( X0, X4, 4, 4), 'k-');
title('Lagrange polynomials');


% part (5) - plotting f(x) and the interpolagin polynomials
F4 = f(X4);
I4 = interp_Lag( X0, X4, F4, N );
figure(5); plot( X0, f(X0), '-b', X4, F4, 'ok', X0, I4, '--r' );
title('Function and Interpolants, N=4');

N = 16;
X16 = linspace(-pi,pi,N+1);
F16 = f(X16);
I16 = interp_Lag( X0, X16, F16, N );
figure(6); plot( X0, f(X0), '-b', X16, F16, 'ok', X0, I16, '--r' );
title('Function and Interpolants, N=16');

% part (d) - examining the errors;

fprintf('\n\n\n');

xx = 0.01;
for n=2:2:20;
  Xn = linspace(-pi,pi,n+1);  
  Yn1 = f( Xn );
  yl1 = interp_Lag( xx, Xn, Yn1, n );
  Yn2 = abs(f( Xn ));
  yl2 = interp_Lag( xx, Xn, Yn2, n );
    
  buf = sprintf('  N=%2d  Err1= %10f  Err12= %10f', ...
		n, (yl1 - f(xx)),  (yl2 - abs(f(xx))));
  disp(buf)
end;
disp(['  The errors decay much slower for g(x), because it is not' ...
      ' differentiable at x=0 '])




disp('  ')
disp('  ------------------------------------------------  ')  
disp('  ')
disp('  ')
disp('  ')
disp('  ')


% Please copy as many headers as necessary

% The function to be minimized
function [z] = F( x, y )
  
  z = cos(y).^2 .* (x - 1).^4 + 100 * (1 - cos(y)).^2;
%  z = x.^2 + 2.0*y.^2;
  
% Gradient of the function
function [df] = dF( x, y )
  
  df = [4 * cos(y).^2 .* (x - 1).^3; ...
	-2 * cos(y) .* (x - 1).^4 .* sin(y) + 200 * (1-cos(y)) .* sin(y)];

%  df = [ 2 * x; 4 * y ];


% Hessian of the function
function [h] = H( x, y )
  
  hxx = 12 * cos(y).^2 .* (x - 1).^2;
  hxy = -8 * cos(y) .* (x - 1).^3 .* sin(y);
  hyy =  2 * sin(y).^2 .* (x - 1).^4 - 2 * cos(y).^2 .* (x-1).^4 + ...
	 200 * sin(y).^2 + 200 .* (1-cos(y)) .* cos(y);
 
%  hxx = 2.0;
%  hyy = 4.0;
%  hxy = 0.0;
  
  h = [ hxx  hxy;
        hxy  hyy ];
  
       
% the first function to be interpolated;

function [y] = f(x) 
  
  y = x .* cos( x.^2 );
 
% Vandermonde interpolation matrix 
function [A] = interp_matrix( x, n )
  
  for i=1:n+1
    for j=1:n+1
      A(i,j) = x(i)^(j-1);
    end;
  end
  
% Lagrange polynomial
function [y] = Lag( x, x0, n, k )
  
  y = 1.0;
  for i=1:n+1
    if ( i ~= k )
      y = y .* ( x - x0(i) ) ./ ( x0(k) - x0(i) );
    end;
  end;
 
% interpolation polynomial
function [y] = interp_Lag( x, x0, b, n );
  
  y = 0.0;
  for i=1:n+1
    y = y + b(i) * Lag( x, x0, n, i);
  end;
