1 ## Copyright (C) 2008-2012 Jaroslav Hajek, Marco Caliari
3 ## This file is part of Octave.
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
20 ## @deftypefn {Function File} {} expm (@var{A})
21 ## Return the exponential of a matrix, defined as the infinite Taylor
25 ## \exp (A) = I + A + {A^2 \over 2!} + {A^3 \over 3!} + \cdots
31 ## expm (A) = I + A + A^2/2! + A^3/3! + @dots{}
35 ## The Taylor series is @emph{not} the way to compute the matrix
36 ## exponential; see Moler and Van Loan, @cite{Nineteen Dubious Ways to
37 ## Compute the Exponential of a Matrix}, SIAM Review, 1978. This routine
38 ## uses Ward's diagonal Pad@'e approximation method with three step
39 ## preconditioning (SIAM Journal on Numerical Analysis, 1977). Diagonal
40 ## Pad@'e approximations are rational polynomials of matrices
42 ## $D_q(A)^{-1}N_q(A)$
54 ## whose Taylor series matches the first
61 ## terms of the Taylor series above; direct evaluation of the Taylor series
62 ## (with the same preconditioning steps) may be desirable in lieu of the
63 ## Pad@'e approximation when
70 ## is ill-conditioned.
71 ## @seealso{logm, sqrtm}
80 if (! ismatrix (A) || ! issquare (A))
81 error ("expm: A must be a square matrix");
87 elseif (strfind (typeinfo (A), "diagonal matrix"))
88 r = diag (exp (diag (A)));
94 A(A == -Inf) = -realmax;
95 trshift = trace (A) / length (A);
100 [d, p, aa] = balance (A);
101 ## FIXME: can we both permute and scale at once? Or should we rather do
104 ## [d, xx, aa] = balance (A, "noperm");
105 ## [xx, p, aa] = balance (aa, "noscal");
106 [f, e] = log2 (norm (aa, "inf"));
111 ## Pade approximation for exp(A).
112 c = [5.0000000000000000e-1,...
113 1.1666666666666667e-1,...
114 1.6666666666666667e-2,...
115 1.6025641025641026e-3,...
116 1.0683760683760684e-4,...
117 4.8562548562548563e-6,...
118 1.3875013875013875e-7,...
119 1.9270852604185938e-9];
123 x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 + id;
124 y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa;
126 r = (x - y) \ (x + y);
128 ## Undo scaling by repeated squaring.
133 ## inverse balancing.
137 ## Inverse trace reduction.
144 %!assert(norm(expm([1 -1;0 1]) - [e -e; 0 e]) < 1e-5);
145 %!assert(expm([1 -1 -1;0 1 -1; 0 0 1]), [e -e -e/2; 0 e -e; 0 0 e], 1e-5);
147 %% Test input validation
150 %!error <expm: A must be a square matrix> expm([1 0;0 1; 2 2]);
152 %!assert (expm (10), expm (10))
153 %!assert (full (expm (eye (3))), expm (full (eye (3))))
154 %!assert (full (expm (10*eye (3))), expm (full (10*eye (3))), 8*eps)