% Steepest descents and conjugate gradients, applied to FD approximation to 
%   u'' = 6*x, u(0) = 0, u(1) = 1
% Shows error convergence for several grid spacings, with initial guess u = x
% for all interior gridpoints.  

  errtolrel = 1e-6;           % Relative error tolerance norm(f-Au)/norm(f)

  a = 0;
  b = 1;

  pmax = 5;
  h = 2.^-(2 + (1:pmax));     % Vector of pmax grid spacings to investigate
  
  for p = 1:pmax
    m = 1/h(p) - 1;    
    maxiter = 100*m^2;        % Relate max number of iterations to 
                              % condition number of A, which is O(m^2), to
                              % ensure convergence of steepest descents.
    x = h(p)*(1:m)';
    f = 6*x;
    f(1) = f(1) - a/h(p)^2;   % Correct f for boundary conditions
    f(m) = f(m) - b/h(p)^2;
    fnorm = norm(f,2);
    uex = x.^3;               % Exact solution
    u0 = x;                   % Initial guess at u

    e = ones(m,1);
    A = spdiags([e -2*e e]/h(p)^2, -1:1, m, m); % A is negative definite, but
                                             % methods also work for this case.
     
    [u,rho] = steepest_descents(A,u0,f,errtolrel,maxiter);  
    resid_rel_sd = sqrt(rho)/fnorm;   % relative residuals after each iteration
    [u,rho] = CG(A,u0,f,errtolrel,maxiter);  
    resid_rel_cg = sqrt(rho)/fnorm;

    if (p==1)

   %   Plot the error convergence of both methods for h = 1/8

      subplot(2,2,1)
      semilogy(1:maxiter,resid_rel_sd,'m-',1:maxiter,resid_rel_cg,'b+');
      axis([0 100 1e-6 1])
      xlabel('Iteration')
      ylabel('Relative residual')
      legend('Steep Desc','Conj Grad',1)
      title('u"=6x, u(0)=0, u(1)=1 with h=1/8')
    end

%    Save required number of iterations for convergence

    niter_sd(p) = length(resid_rel_sd(resid_rel_sd>0));
    niter_cg(p) = length(resid_rel_cg(resid_rel_cg>0));
  end

%  Plot iterations required for convergence vs. h

  subplot(2,2,2)
  loglog(h,niter_sd,'m*',h,niter_cg,'b+')
  xlabel('h')
  ylabel('No. of iterations')
  legend('Steep Desc','Conj Grad',3)
