|
8393
|
1 ## Copyright (C) 2008 Jaroslav Hajek, Marco Caliari |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
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. |
|
|
9 ## |
|
|
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. |
|
|
14 ## |
|
|
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/>. |
|
|
18 |
|
|
19 ##-*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {} expm (@var{a}) |
|
|
21 ## Return the exponential of a matrix, defined as the infinite Taylor |
|
|
22 ## series |
|
|
23 ## @iftex |
|
|
24 ## @tex |
|
|
25 ## $$ |
|
|
26 ## \exp (A) = I + A + {A^2 \over 2!} + {A^3 \over 3!} + \cdots |
|
|
27 ## $$ |
|
|
28 ## @end tex |
|
|
29 ## @end iftex |
|
|
30 ## @ifinfo |
|
|
31 ## |
|
|
32 ## @example |
|
|
33 ## expm(a) = I + a + a^2/2! + a^3/3! + ... |
|
|
34 ## @end example |
|
|
35 ## |
|
|
36 ## @end ifinfo |
|
|
37 ## The Taylor series is @emph{not} the way to compute the matrix |
|
|
38 ## exponential; see Moler and Van Loan, @cite{Nineteen Dubious Ways to |
|
|
39 ## Compute the Exponential of a Matrix}, SIAM Review, 1978. This routine |
|
|
40 ## uses Ward's diagonal |
|
|
41 ## @iftex |
|
|
42 ## @tex |
|
|
43 ## Pad\'e |
|
|
44 ## @end tex |
|
|
45 ## @end iftex |
|
|
46 ## @ifinfo |
|
|
47 ## Pade' |
|
|
48 ## @end ifinfo |
|
|
49 ## approximation method with three step preconditioning (SIAM Journal on |
|
|
50 ## Numerical Analysis, 1977). Diagonal |
|
|
51 ## @iftex |
|
|
52 ## @tex |
|
|
53 ## Pad\'e |
|
|
54 ## @end tex |
|
|
55 ## @end iftex |
|
|
56 ## @ifinfo |
|
|
57 ## Pade' |
|
|
58 ## @end ifinfo |
|
|
59 ## approximations are rational polynomials of matrices |
|
|
60 ## @iftex |
|
|
61 ## @tex |
|
|
62 ## $D_q(a)^{-1}N_q(a)$ |
|
|
63 ## @end tex |
|
|
64 ## @end iftex |
|
|
65 ## @ifinfo |
|
|
66 ## |
|
|
67 ## @example |
|
|
68 ## -1 |
|
|
69 ## D (a) N (a) |
|
|
70 ## @end example |
|
|
71 ## |
|
|
72 ## @end ifinfo |
|
|
73 ## whose Taylor series matches the first |
|
|
74 ## @iftex |
|
|
75 ## @tex |
|
|
76 ## $2 q + 1 $ |
|
|
77 ## @end tex |
|
|
78 ## @end iftex |
|
|
79 ## @ifinfo |
|
|
80 ## @code{2q+1} |
|
|
81 ## @end ifinfo |
|
|
82 ## terms of the Taylor series above; direct evaluation of the Taylor series |
|
|
83 ## (with the same preconditioning steps) may be desirable in lieu of the |
|
|
84 ## @iftex |
|
|
85 ## @tex |
|
|
86 ## Pad\'e |
|
|
87 ## @end tex |
|
|
88 ## @end iftex |
|
|
89 ## @ifinfo |
|
|
90 ## Pade' |
|
|
91 ## @end ifinfo |
|
|
92 ## approximation when |
|
|
93 ## @iftex |
|
|
94 ## @tex |
|
|
95 ## $D_q(a)$ |
|
|
96 ## @end tex |
|
|
97 ## @end iftex |
|
|
98 ## @ifinfo |
|
|
99 ## @code{Dq(a)} |
|
|
100 ## @end ifinfo |
|
|
101 ## is ill-conditioned. |
|
|
102 ## @end deftypefn |
|
|
103 |
|
|
104 function r = expm (a) |
|
|
105 |
|
|
106 if (! ismatrix (a) || ! issquare (a)) |
|
|
107 error ("expm requires a square matrix") |
|
|
108 endif |
|
|
109 |
|
|
110 n = rows (a); |
|
|
111 # trace reduction |
|
|
112 a(a == -Inf) = -realmax; |
|
|
113 trshift = trace (a) / length (a); |
|
|
114 if (trshift > 0) |
|
|
115 a -= trshift*eye (n); |
|
|
116 endif |
|
|
117 # balancing |
|
|
118 [p, d, aa] = balance (a); |
|
|
119 # FIXME: can we both permute and scale at once? Or should we rather do this: |
|
|
120 # [p, xx, aa] = balance (a, "noscal"); |
|
|
121 # [xx, d, aa] = balance (aa, "noperm"); |
|
|
122 [f, e] = log2 (norm (aa, "inf")); |
|
|
123 s = max (0, e); |
|
|
124 s = min (s, 1023); |
|
|
125 aa *= 2^(-s); |
|
|
126 |
|
|
127 # Pade approximation for exp(A) |
|
|
128 c = [5.0000000000000000e-1,... |
|
|
129 1.1666666666666667e-1,... |
|
|
130 1.6666666666666667e-2,... |
|
|
131 1.6025641025641026e-3,... |
|
|
132 1.0683760683760684e-4,... |
|
|
133 4.8562548562548563e-6,... |
|
|
134 1.3875013875013875e-7,... |
|
|
135 1.9270852604185938e-9]; |
|
|
136 |
|
|
137 a2 = aa^2; |
|
|
138 id = eye (n); |
|
|
139 x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 + id; |
|
|
140 y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa; |
|
|
141 |
|
|
142 r = (x - y) \ (x + y); |
|
|
143 |
|
|
144 # Undo scaling by repeated squaring |
|
|
145 for k = 1:s |
|
|
146 r ^= 2; |
|
|
147 endfor |
|
|
148 |
|
|
149 # inverse balancing |
|
|
150 ds = diag (s); |
|
|
151 r = ds * r / ds; |
|
|
152 r = r(p, p); |
|
|
153 # inverse trace reduction |
|
|
154 if (trshift >0) |
|
|
155 r *= exp (trshift); |
|
|
156 endif |
|
|
157 |
|
|
158 endfunction |
