Mercurial > hg > octave-lyh
annotate scripts/linear-algebra/expm.m @ 9819:84398271118c
fix typo in expm
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Sat, 14 Nov 2009 07:02:30 +0100 |
| parents | f0c3d3fc4903 |
| children | d1978e7364ad |
| 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 | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
19 ## -*- texinfo -*- |
| 8393 | 20 ## @deftypefn {Function File} {} expm (@var{a}) |
| 21 ## Return the exponential of a matrix, defined as the infinite Taylor | |
| 22 ## series | |
| 23 ## @tex | |
| 24 ## $$ | |
| 25 ## \exp (A) = I + A + {A^2 \over 2!} + {A^3 \over 3!} + \cdots | |
| 26 ## $$ | |
| 27 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
28 ## @ifnottex |
| 8393 | 29 ## |
| 30 ## @example | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
31 ## expm(a) = I + a + a^2/2! + a^3/3! + @dots{} |
| 8393 | 32 ## @end example |
| 33 ## | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
34 ## @end ifnottex |
| 8393 | 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 | |
| 39 ## @tex | |
| 40 ## Pad\'e | |
| 41 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
42 ## @ifnottex |
| 8393 | 43 ## Pade' |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
44 ## @end ifnottex |
| 8393 | 45 ## approximation method with three step preconditioning (SIAM Journal on |
| 46 ## Numerical Analysis, 1977). Diagonal | |
| 47 ## @tex | |
| 48 ## Pad\'e | |
| 49 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
50 ## @ifnottex |
| 8393 | 51 ## Pade' |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
52 ## @end ifnottex |
| 8393 | 53 ## approximations are rational polynomials of matrices |
| 54 ## @tex | |
| 55 ## $D_q(a)^{-1}N_q(a)$ | |
| 56 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
57 ## @ifnottex |
| 8393 | 58 ## |
| 59 ## @example | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
60 ## @group |
| 8393 | 61 ## -1 |
| 62 ## D (a) N (a) | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
63 ## @end group |
| 8393 | 64 ## @end example |
| 65 ## | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
66 ## @end ifnottex |
| 8393 | 67 ## whose Taylor series matches the first |
| 68 ## @tex | |
| 69 ## $2 q + 1 $ | |
| 70 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
71 ## @ifnottex |
| 8393 | 72 ## @code{2q+1} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
73 ## @end ifnottex |
| 8393 | 74 ## terms of the Taylor series above; direct evaluation of the Taylor series |
| 75 ## (with the same preconditioning steps) may be desirable in lieu of the | |
| 76 ## @tex | |
| 77 ## Pad\'e | |
| 78 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
79 ## @ifnottex |
| 8393 | 80 ## Pade' |
|
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 ## approximation when |
| 83 ## @tex | |
| 84 ## $D_q(a)$ | |
| 85 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
86 ## @ifnottex |
| 8393 | 87 ## @code{Dq(a)} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
88 ## @end ifnottex |
| 8393 | 89 ## is ill-conditioned. |
| 90 ## @end deftypefn | |
| 91 | |
| 92 function r = expm (a) | |
| 93 | |
| 94 if (! ismatrix (a) || ! issquare (a)) | |
| 8664 | 95 error ("expm requires a square matrix"); |
| 8393 | 96 endif |
| 97 | |
| 98 n = rows (a); | |
| 8506 | 99 ## Trace reduction. |
| 8393 | 100 a(a == -Inf) = -realmax; |
| 101 trshift = trace (a) / length (a); | |
| 102 if (trshift > 0) | |
| 103 a -= trshift*eye (n); | |
| 104 endif | |
| 8506 | 105 ## Balancing. |
|
9048
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
106 [d, p, aa] = balance (a); |
| 8506 | 107 ## FIXME: can we both permute and scale at once? Or should we rather do |
| 108 ## this: | |
| 109 ## | |
|
9048
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
110 ## [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
|
111 ## [xx, p, aa] = balance (aa, "noscal"); |
| 8393 | 112 [f, e] = log2 (norm (aa, "inf")); |
| 113 s = max (0, e); | |
| 114 s = min (s, 1023); | |
| 115 aa *= 2^(-s); | |
| 116 | |
| 8506 | 117 ## Pade approximation for exp(A). |
| 8393 | 118 c = [5.0000000000000000e-1,... |
| 119 1.1666666666666667e-1,... | |
| 120 1.6666666666666667e-2,... | |
| 121 1.6025641025641026e-3,... | |
| 122 1.0683760683760684e-4,... | |
| 123 4.8562548562548563e-6,... | |
| 124 1.3875013875013875e-7,... | |
| 125 1.9270852604185938e-9]; | |
| 126 | |
| 127 a2 = aa^2; | |
| 128 id = eye (n); | |
| 129 x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 + id; | |
| 130 y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa; | |
| 131 | |
| 132 r = (x - y) \ (x + y); | |
| 133 | |
| 8506 | 134 ## Undo scaling by repeated squaring. |
| 135 for k = 1:s | |
| 8393 | 136 r ^= 2; |
| 137 endfor | |
| 138 | |
| 8506 | 139 ## inverse balancing. |
| 8757 | 140 d = diag (d); |
| 141 r = d * r / d; | |
| 9819 | 142 r(p, p) = r; |
| 8506 | 143 ## Inverse trace reduction. |
| 8393 | 144 if (trshift >0) |
| 145 r *= exp (trshift); | |
| 146 endif | |
| 147 | |
| 148 endfunction |