% Chebfun code to compute the eigenvalues of the advection-diffusion operator
% L u = -u'' + c u' on L^2[0,pi] with homogeneous Dirichlet boundary conditions
% by applying the Arnoldi (Lanczos) algorithm directly to the inverse operator,
% implemented via variation of parameters.  The exact eigenvalues are n^2 + c^2/4.

% Rice Eigenvalue Clinic, August 2009

 c = 4;                          % advection coefficient; c <= 4 is ok
 m = 10;                         % maximum number of iterations
 v0 = chebfun('x.^0',[0 pi]);    % starting vector
 xfun = chebfun('x',[0 pi]);     % the function "x"
 V = v0/norm(v0);                % normalize the starting vector
 H = zeros(m+1,m);               % H will contain the upper Hessenberg matrix

 for k=1:m                       % Arnoldi process
    w = cumsum(V(:,k)/c)-cumsum(V(:,k)./(c*exp(c*xfun))).*exp(c*xfun);
    c2 = (w(pi)-w(0))/(exp(c*pi)-1);
    c1 = w(0)-c2;
    w = w-c1-c2*exp(c*xfun);
    for j=1:k                    %    Gram-Schmidt
        H(j,k) = V(:,j)'*w;
        w = w-H(j,k)*V(:,j);
    end
    H(k+1,k) = norm(w);
    V(:,k+1) = w/H(k+1,k);       %    append new basis vector to V

    fprintf('\n-----------------------------------------\n')
    fprintf('step %d\n', k);
    fprintf('-----------------------------------------\n')
    fprintf('      Ritz value       relative error \n')
    fprintf('-----------------------------------------\n')
    ew = sort(1./eig(H(1:k,1:k)));
    for j=1:k
       fprintf('%20.15f    %12.5e\n', ew(j), ((j^2+c^2/4)-ew(j))/(j^2))
    end
 end
 
 fprintf('\nDeparture from orthogonality of Arnoldi basis functions: %12.5e\n\n', ...
     norm(V'*V-eye(m+1)))