home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / EIGMOVIE.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.6 KB  |  73 lines
  1. function eigmovie(A)
  2. %EIGMOVIE Symmetric eigenvalue movie.
  3. %    EIGMOVIE(A) shows a "movie" that depicts the computation of
  4. %    the eigenvalues of a symmetric matrix.
  5. %    EIGMOVIE, by itself, uses A = pascal(6)/2.
  6.  
  7. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  8.  
  9. if nargin < 1
  10.    A = pascal(6)/2;
  11. end
  12.  
  13. n = max(size(A));
  14. A = (A + A')/2;
  15. more off
  16. clc
  17. format short
  18. home
  19. disp('Symmetrized input matrix')
  20. A
  21. disp('Pause. Press any key to continue.'), pause
  22. I = eye(size(A));
  23.  
  24. for k = 2:n-1
  25.    home
  26.    disp(['Householder similarity transformation ' int2str(k-1)]) 
  27.    u = A(:,k-1);
  28.    u(1:k-1) = zeros(k-1,1);
  29.    u = u/norm(u);
  30.    if u(k) < 0, u = -u; end
  31.    u(k) = u(k) + 1;
  32.    beta = -u(k);  % -(norm(u)^2)/2
  33.    v = A*u/beta;
  34.    gamma = v'*u/(2*beta);
  35.    v = v + gamma*u;
  36.    A  =  A + u*v' + v*u';
  37.    j = (k+1):n;  z = zeros(n-k,1);
  38.    A(j,k-1) = z;  A(k-1,j) = z';
  39.    A
  40.    disp('Pause                            '), pause
  41. end
  42. home
  43.  
  44. clc
  45. it = 0;
  46. while n > 0
  47.    home
  48.    disp(['Symmetric tridiagonal QR iteration ' int2str(it) '   '])
  49.    d = diag(A); e = diag(A,-1);
  50.    A = diag(d) + diag(e,-1) + diag(e,1)
  51.    if n > 1, k = n-1:n; else k = 1:2; end
  52.    S = A(k,k);
  53.    disp('Lower 2-by-2 = ')
  54.    disp(' ')
  55.    format long e
  56.    disp(S)
  57.    format short
  58.    disp('Pause'), pause
  59.    if n == 1, break, end
  60.    if abs(S(2,1)) < eps*(abs(S(1,1)) + abs(S(2,2)))
  61.       A(n,n-1) = 0;  A(n-1,n) = 0;
  62.       n = n-1;
  63.    else
  64.       shift = eig(S);
  65.       if abs(shift(1)-A(n,n)) < abs(shift(2)-A(n,n))
  66.          shift = shift(1); else, shift = shift(2); end
  67.       [Q,R] = qr(A-shift*I);
  68.       A = R*Q+shift*I;
  69.    end
  70.    it = it + 1;
  71. end
  72. end
  73.