function hw2
  
close all
clear all


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Question #1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Nmax = 10^5;
N = 4;
i = 1;

%I1e = 2*pi*besselI(0, 1);
I1e = besseli(0, 1);
%I2e = 12.4175160714222;
I2e = 1.97630906368989;

while ( N < Nmax )
  x = linspace(-0.5,0.5,N+1);
  x0 = x(1:N);
  dx = 1.0 / N;
  y1 = fun1( x0 );
  y2 = fun2( x0 );
  
  I1 = dx * sum( y1 );
  I2 = dx * sum( y2 );
 
  Err1(i,1) = N;
  Err1(i,2) = abs( (I1 - I1e) / I1e );
  Err1(i,3) = abs( (I2 - I2e) / I2e );

  i = i + 1;
  N = N * 2;
end


figure(1)
loglog(Err1(:,1),Err1(:,2),'-xb');
hold on
loglog(Err1(:,1),Err1(:,3),'-or');
loglog(Err1(:,1),1.0e2*Err1(:,1).^(-2),'-k');
legend('exp(sin(x))','exp(sin(x/2))','C h^{-2}');
title('Question 1')

disp('== Question #1:')
disp('   f_1 is real-analytic on the periodic interval [-1/2,1/2], ')
disp('   hence the periodic trapezoidal quadrature converges spectrally.')
disp('    ')
disp('   f_2 is only continuous on the periodic interval [-1/2,1/2], ')
disp('   hence the periodic trapezoidal quadrature converges algebraically')
disp('   (with second order).')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Question #1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


N0 = 10^5;
N = 100;

xx0 = linspace(-0.5,0.5,N0+1);
xx = xx0(1:end-1);
dx0 = 1.0 / N0;

f1  =  g1( xx );
f1s = zeros(1,N0);
f1_2 = norm( f1, 2 );
f1_inf = norm( f1, inf );

f2  =  g2( xx );
f2s = zeros(1,N0);
f2_2 = norm( f2, 2 );
f2_inf = norm( f2, inf );


figure(5)
plot(xx,f1,'-r');
title('g1')
hold on

figure(6)
plot(xx,f2,'-r');
title('g2')
hold on


for k=1:N+1
% computing the expansion coefficients via quadratures;
  if ( k == 1 ) 
    phik = ones(1,N0);
  else
    phik = sqrt(2.0) * cos( 2.0 * pi * (k-1) * xx );
  end;
  
  f1k(k) = dx0 * sum( f1 .* phik );
  f2k(k) = dx0 * sum( f2 .* phik );  
  
% constructing the series approximation;
  f1s = f1s + f1k(k) * phik;
  f2s = f2s + f2k(k) * phik;
  if ( mod(k,20) == 1 )
    figure(5)
    plot(xx(1:100:end),f1s(1:100:end),'--b');
    figure(6)
    plot(xx(1:100:end),f2s(1:100:end),'--b');
  end;
  
% computing the error norms;
  E2(k,1) = norm( f1s - f1, 2) / f1_2;
  E2(k,2) = norm( f2s - f2, 2) / f2_2;
  Einf(k,1) = norm( f1s - f1, inf) / f1_inf;  
  Einf(k,2) = norm( f2s - f2, inf) / f2_inf;  
end;


figure(2)
loglog(1:N+1,abs(f1k(:)),'-xb',1:N+1,abs(f2k(:)),'-or');
legend('f1k','f2k');
title('Question 2 - Fourier coefficients')


figure(3)
loglog(1:N+1,E2(:,1),'-xb',1:N+1,E2(:,2),'-or');
legend('g1','g2');
title('Question 2 - L^2 norm')


figure(4)
loglog(1:N+1,Einf(:,1),'-xb',1:N+1,Einf(:,2),'-or');
legend('g1','g2');
title('Question 2 - L^{inf} norm')


fprintf('\n\n\n')

disp('== Question #2:')
disp('   The L^2 errors vanish with increase of N for both functions g_1 and g_2,')
disp('   but the convergence is faster for g_2 because it is smoother.')
disp('    ')
disp('   The L^inf errors vanish with increase of N only for the function g_2.')
disp('   The function g_1 is discontinuous and therefore not sufficiently')
disp('   regular to ensure pointwise (uniform) convergence of the Fourier series')
disp('   (for g_1 convergence is achieved only in the L^2 sense)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Question #3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f1s = zeros(1,N0);
f1s2 = zeros(1,N0);
epsilon = 0.1;
figure(7), hold on
N = 1000;

for k=1:N+1
% computing the expansion coefficients via quadratures;
  if ( k == 1 ) 
    phik = ones(1,N0);
  else
    phik = sqrt(2.0) * cos( 2.0 * pi * (k-1) * xx );
  end;
  
  f1k(k) = dx0 * sum( f1 .* phik );
  if ( k == 1 )
    f1k2(k) = f1k(k);
  else    
    f1k2(k) = f1k(k) / (2.0 * pi * (k-1))^epsilon;
  end;
  
% constructing the series approximation;
  f1s = f1s + f1k(k) * phik;
  f1s2 = f1s2 + f1k2(k) * phik;
  
% $$$   if ( mod(k,10) == 1 )
% $$$     plot(xx(1:10:end),f1s(1:10:end),'--b');
% $$$   end;
  
end;

plot(xx(1:10:end),f1s(1:10:end),'-b',xx(1:10:end),f1s2(1:10:end),'-r');
legend('original (Runges phenomenon)','smoothed (no Runges phenomenon)')
title('Question 3')


fprintf('\n\n\n')

disp('== Question #3:')
disp('   To have square-integrable derivatives of order 1/2+eps,')
disp('   the function must have Fourier coefficients decaying ')
disp('   at least as O(k^(-1-eps)). We thus take the Fourier coefficients')
disp('   of g_1, which decay as O(k^(-1)), and make them vanish a bit faster')





function [f] = fun1( x )
  
  f = exp(sin(2.0*pi*x));
  
  
  
function [f] = fun2( x )
  
  f = exp(abs(sin(2.0*pi*x)));

  
function [f] = g1( x )
  
  for i=1:length(x)    
    if ( abs(x(i)) >= 0.25 )
      f(i) = 1.0;
    else
      f(i) = - 1.0;
    end;
  end;
   
  
function [f] = g2( x )
  
  for i=1:length(x)    
    if ( x(i) < 0.0 )
      x(i) = - x(i);
    end;
    
    if ( x(i) <= 0.25 )
      f(i) = -6.0 * x(i) + 1.0;
    else
      f(i) = 4.0 * x(i) - 1.5;
    end;
  end;
   
