Mercurial > hg > octave-lyh
annotate scripts/linear-algebra/expm.m @ 9048:867d5d1aed06
swap out args in balance for M*b compat
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 27 Mar 2009 14:28:10 +0100 |
| parents | eb63fbe60fab |
| children | 8207b833557f |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2008, 2009 Jaroslav Hajek, Marco Caliari |
| 8393 | 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 | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
30 ## @ifnottex |
| 8393 | 31 ## |
| 32 ## @example | |
| 33 ## expm(a) = I + a + a^2/2! + a^3/3! + ... | |
| 34 ## @end example | |
| 35 ## | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
36 ## @end ifnottex |
| 8393 | 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 | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
46 ## @ifnottex |
| 8393 | 47 ## Pade' |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
48 ## @end ifnottex |
| 8393 | 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 | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
56 ## @ifnottex |
| 8393 | 57 ## Pade' |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
58 ## @end ifnottex |
| 8393 | 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 | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
65 ## @ifnottex |
| 8393 | 66 ## |
| 67 ## @example | |
| 68 ## -1 | |
| 69 ## D (a) N (a) | |
| 70 ## @end example | |
| 71 ## | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
72 ## @end ifnottex |
| 8393 | 73 ## whose Taylor series matches the first |
| 74 ## @iftex | |
| 75 ## @tex | |
| 76 ## $2 q + 1 $ | |
| 77 ## @end tex | |
| 78 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
79 ## @ifnottex |
| 8393 | 80 ## @code{2q+1} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
81 ## @end ifnottex |
| 8393 | 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 | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
89 ## @ifnottex |
| 8393 | 90 ## Pade' |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
91 ## @end ifnottex |
| 8393 | 92 ## approximation when |
| 93 ## @iftex | |
| 94 ## @tex | |
| 95 ## $D_q(a)$ | |
| 96 ## @end tex | |
| 97 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
98 ## @ifnottex |
| 8393 | 99 ## @code{Dq(a)} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
100 ## @end ifnottex |
| 8393 | 101 ## is ill-conditioned. |
| 102 ## @end deftypefn | |
| 103 | |
| 104 function r = expm (a) | |
| 105 | |
| 106 if (! ismatrix (a) || ! issquare (a)) | |
| 8664 | 107 error ("expm requires a square matrix"); |
| 8393 | 108 endif |
| 109 | |
| 110 n = rows (a); | |
| 8506 | 111 ## Trace reduction. |
| 8393 | 112 a(a == -Inf) = -realmax; |
| 113 trshift = trace (a) / length (a); | |
| 114 if (trshift > 0) | |
| 115 a -= trshift*eye (n); | |
| 116 endif | |
| 8506 | 117 ## Balancing. |
|
9048
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
118 [d, p, aa] = balance (a); |
| 8506 | 119 ## FIXME: can we both permute and scale at once? Or should we rather do |
| 120 ## this: | |
| 121 ## | |
|
9048
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
122 ## [d, xx, aa] = balance (a, "noperm"); |
|
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
123 ## [xx, p, aa] = balance (aa, "noscal"); |
| 8393 | 124 [f, e] = log2 (norm (aa, "inf")); |
| 125 s = max (0, e); | |
| 126 s = min (s, 1023); | |
| 127 aa *= 2^(-s); | |
| 128 | |
| 8506 | 129 ## Pade approximation for exp(A). |
| 8393 | 130 c = [5.0000000000000000e-1,... |
| 131 1.1666666666666667e-1,... | |
| 132 1.6666666666666667e-2,... | |
| 133 1.6025641025641026e-3,... | |
| 134 1.0683760683760684e-4,... | |
| 135 4.8562548562548563e-6,... | |
| 136 1.3875013875013875e-7,... | |
| 137 1.9270852604185938e-9]; | |
| 138 | |
| 139 a2 = aa^2; | |
| 140 id = eye (n); | |
| 141 x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 + id; | |
| 142 y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa; | |
| 143 | |
| 144 r = (x - y) \ (x + y); | |
| 145 | |
| 8506 | 146 ## Undo scaling by repeated squaring. |
| 147 for k = 1:s | |
| 8393 | 148 r ^= 2; |
| 149 endfor | |
| 150 | |
| 8506 | 151 ## inverse balancing. |
| 8757 | 152 d = diag (d); |
| 153 r = d * r / d; | |
| 8393 | 154 r = r(p, p); |
| 8506 | 155 ## Inverse trace reduction. |
| 8393 | 156 if (trshift >0) |
| 157 r *= exp (trshift); | |
| 158 endif | |
| 159 | |
| 160 endfunction |