home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / LOGM2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.3 KB  |  55 lines
  1. function [Y, warning]  = logm2(A, pert)
  2. %LOGM2    LOGM2(X)  is the matrix natural logarithm of  X .    Complex
  3. %    results  are  produced  if  X  has nonpositive eigenvalues.
  4. %    LOGM2  can be thought of as being computed using eigenvalues
  5. %    D  and eigenvectors  V,  such that if  [V,D] = EIG(X)  then  
  6. %    LOGM2(X) = V*log(D)/V.  
  7.  
  8. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  9.  
  10. warning = 0;
  11. if ~length(A), Y = []; return; end
  12. [Q,T] = schur(A);
  13. [Q,T] = rsf2csf(Q,T);
  14. dT = diag(T);
  15. F = log(diag(T));
  16.  
  17. % Set log of zero to large negative number
  18. find_vec = find(~finite(F));
  19. F(find_vec) = -ones(length(find_vec),1).*(1/eps);
  20.  
  21. F = diag(F);
  22. tol = eps*norm(T,1);
  23. n = max(size(A));
  24. for p = 1:n-1
  25.    for i = 1:n-p
  26.       j = i+p;
  27.       s = T(i,j)*(F(j,j)-F(i,i));
  28.       if p > 1
  29.          k = i+1:j-1;
  30.          s = s + T(i,k)*F(k,j) - F(i,k)*T(k,j);
  31.       end
  32.       d = T(j,j) - T(i,i);
  33.       if abs(d) <= tol,
  34.      warning = 1;
  35.          d = tol;
  36.       end
  37.       F(i,j) = s/d;
  38.    end
  39. end
  40. Y = Q*F*Q';
  41.  
  42. % If diagonal elements are the same then check accuracy against expm.
  43. % If accuracy is poor then perturb matrix and call logm2 again.
  44. if (warning)
  45.     tol = tol + eps;
  46.     if  norm(expm(Y) - A, 'inf') > 1e5*tol 
  47.         if nargin ==1, 
  48.             pert = tol; 
  49.         else
  50.             pert = pert * 100
  51.         end
  52.         Y = logm2(A + pert*randn(n,n), pert);
  53.     end
  54. end 
  55.