home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / MATFUN.DI$ / FUNM.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.0 KB  |  63 lines
  1. function [F,esterr] = funm(A,fun)
  2. %FUNM    Evaluate function of a matrix.
  3. %    F = FUNM(A,'fun') for a matrix argument A evaluates the matrix
  4. %    function specified by 'fun'.   For example, FUNM(A,'sin') is
  5. %    the matrix sine.  For matrix exponentials, logarithms and
  6. %    square roots, use EXPM(A), LOGM(A) and SQRTM(A) instead.
  7. %
  8. %    FUNM uses a potentially unstable algorithm.  If A is close to a
  9. %    matrix with multiple eigenvalues and poorly conditioned eigenvectors,
  10. %    FUNM may produce inaccurate results.  An attempt is made to detect
  11. %    this situation and print a warning message.  The error detector is
  12. %    sometimes too sensitive and a message is printed even though the
  13. %    the computed result is accurate.
  14. %
  15. %    [F,ESTERR] = FUNM(A,'fun') does not print any message, but returns
  16. %    a very rough estimate of the relative error in the computed result.
  17. %
  18. %    If A is symmetric or Hermitian, then its Schur form is diagonal and
  19. %    FUNM will be able to produce an accurate result.
  20. %
  21. %    S = SQRTM(A) and L = LOGM(A) use FUNM to do their computations,
  22. %    but they can get more reliable error estimates by comparing S*S
  23. %    and EXPM(L) with A.  E = EXPM(A) uses a completely different
  24. %    algorithm.
  25. %
  26. %    See also EXPM, SQRTM, LOGM.
  27.  
  28. %    C.B. Moler 12-2-85, 7-21-86, 7-11-92.
  29. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  30.  
  31. %   Parlett's method.  See Golub and VanLoan (1983), p. 384.
  32.  
  33. [Q,T] = schur(A);
  34. [Q,T] = rsf2csf(Q,T);
  35. F = eval(['diag(',fun,'(diag(T)))']);
  36. n = max(size(A));
  37. dmin = abs(T(1,1));
  38. for p = 1:n-1
  39.    for i = 1:n-p
  40.       j = i+p;
  41.       s = T(i,j)*(F(j,j)-F(i,i));
  42.       if p > 1
  43.          k = i+1:j-1;
  44.          s = s + T(i,k)*F(k,j) - F(i,k)*T(k,j);
  45.       end
  46.       d = T(j,j) - T(i,i);
  47.       if d ~= 0 
  48.          s = s/d;
  49.       end
  50.       F(i,j) = s;
  51.       dmin = min(dmin,abs(d));
  52.    end
  53. end
  54. F = Q*F*Q';
  55.  
  56. if dmin == 0, dmin = eps; end
  57. esterr = min(1,max(eps,(eps/dmin)*norm(triu(T,1),1)));
  58. if (nargout < 2) & (esterr > 1000*eps)
  59.    disp(' ')
  60.    disp(['WARNING: Result from FUNM may be inaccurate.' ...
  61.           ' esterr = ' num2str(esterr)])
  62. end
  63.