%  Here is a really simple illustration of POD on
%  a linear problem   dot(x) = Ax,  x(0) = b.
%  
%  A  = Q(C circle(+) D)Q'
%  
%  Where D = diag(-[3:n])
%  
%  and      C = [ alpha  beta^2;
%                -1 alpha]
%  
%  with alpha = slightly negative  and beta >> abs(alpha)
%  
%  So there are two mildly decaying oscillatory components
%  (corresponding to nonorthogonal eigenvectors -- see 
%  the example in Spectra and Pseudospectra, p. 456)
%  and the rest are rapidly decaying
%  
%  This script illustrates the excellent reproduction
%  of the full solution from the reduced dimension solution
%  both in norm and  in components as a function of time
%  in each of the leading k  SVD co-ordinates.
%  
%  Input   n  =    dimension    (500 recommended)
%         tau =  drop tolerance  (.00001 recommended)
%  
%  ----------------------------------------------
%  D.C. Sorensen
%  21 Sep. 09
%  This test problem due to M. Embree 21 Sep. 09
% 


   n = input(' dimension ')
%
%  e.g. n = 100
%
   tau = input(' drop tolerance ')
%
%  e.g. tau = .00001 
%

   S = -diag([1:n]);
   alpha = -.0001;   beta = 100;
   
   S(1:2,1:2)  = [alpha beta^2;
                  -1 alpha];

   [Q,R] = qr(randn(n));
   A = Q*S*Q';

   h = .01;
   N = 100;

   b = randn(n,1);
   I = eye(n);
   hh = .5*h;

   [Q,R] = qr(I - hh*A);
   x = b;
   X = [x];
   t = [0];
%
%  Integrate dot(x) = Ax, x(0) = b
%  using trapezoidal rule.
%
%  Record the trajectory in X
%
   for j = 1:N,  
       x = R\(Q'*((I + hh*A)*x));
       X = [X x];
       t = [t j*h];
   end

   [U,S,V] = svd(X);
   s = diag(S);

%
%  Determine the level k of SVD truncation
%  required to meet the drop tolerance tau
%
   k = 1;
   while( s(k) > tau*s(1)),
        k = k+1; 
   end 
   Xk = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';

%
%  Get the best rank k  approximation Xk to X
%

   Xdiff_SVD = norm(X - Xk)/s(1)
   
   Xk = [];

%
%  construct the reduced dimension problem
%
   Ak = U(:,1:k)'*A*U(:,1:k);
   b1 = U(:,1:k)'*b;

 
   Ik = eye(k);
   [Qk,Rk] = qr(Ik - hh*Ak);
   xk = b1;
   Xk = [xk];
   
%
%  Integrate dot(xk) = Ak xk, xk(0) = b1
%  using trapezoidal rule.
%
%  Record the componentes trajectory in Xk
%
   for j = 1:N,  
       xk = Rk\(Qk'*((Ik + hh*Ak)*xk));
       Xk = [Xk xk];
   end
%
%  Get an approximation  XXk  to   X   
%  in the  U(:,1:k)  basis
%
   XXk = U(:,1:k)*Xk;

   Xdiff_POD = norm(X - XXk)/s(1)
   Order_of_Reduced_System = k
%
%  Display components as functions of time
%  in the truncated SVD basis  Ul(:,1:k)
%
   for i = 1:k, 
       clf
       plot(t,Xk(i,:), 'b-','linewidth', 2.5)
       hold on
       plot(t, U(:,i)'*X,'r--','linewidth', 2.5)
       legend('Reduced', 'Full')
       hold off
       disp('strike a key')
       pause
   end