home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / MATFUN.DI$ / EXPM1.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  823 b   |  40 lines
  1. function E = expm1(A)
  2. %EXPM1    Matrix exponential via Pade approximation.
  3. %    E = EXPM1(A) is an M-file implementation of the built-in
  4. %    algorithm used by MATLAB for the matrix exponential.
  5. %    See Golub and Van Loan, Matrix Computations, Algorithm 11.3-1.
  6. %
  7. %    See also EXPM, EXPM2, EXPM3.
  8.  
  9. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  10.  
  11. % Scale A by power of 2 so that its norm is < 1/2 .
  12. [f,e] = log2(norm(A,'inf'));
  13. s = max(0,e+1);
  14. A = A/2^s;
  15.  
  16. % Pade approximation for exp(A)
  17. X = A; 
  18. c = 1/2;
  19. E = eye(size(A)) + c*A;
  20. D = eye(size(A)) - c*A;
  21. q = 6;
  22. p = 1;
  23. for k = 2:q
  24.    c = c * (q-k+1) / (k*(2*q-k+1));
  25.    X = A*X;
  26.    cX = c*X;
  27.    E = E + cX;
  28.    if p
  29.      D = D + cX;
  30.    else
  31.      D = D - cX;
  32.    end
  33.    p = ~p;
  34. end
  35. E = D\E;
  36.  
  37. % Undo scaling by repeated squaring
  38. for k=1:s, E = E*E; end
  39.  
  40.