Mercurial > hg > octave-lyh
annotate scripts/linear-algebra/expm.m @ 9065:8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Grammarcheck input .txi files
Spellcheck .texi files
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Sat, 28 Mar 2009 17:43:22 -0700 |
| parents | 867d5d1aed06 |
| children | f0c3d3fc4903 |
| 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 ## @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 | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
33 ## expm(a) = I + a + a^2/2! + a^3/3! + @dots{} |
| 8393 | 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 | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
68 ## @group |
| 8393 | 69 ## -1 |
| 70 ## D (a) N (a) | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
71 ## @end group |
| 8393 | 72 ## @end example |
| 73 ## | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
74 ## @end ifnottex |
| 8393 | 75 ## whose Taylor series matches the first |
| 76 ## @iftex | |
| 77 ## @tex | |
| 78 ## $2 q + 1 $ | |
| 79 ## @end tex | |
| 80 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
81 ## @ifnottex |
| 8393 | 82 ## @code{2q+1} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
83 ## @end ifnottex |
| 8393 | 84 ## terms of the Taylor series above; direct evaluation of the Taylor series |
| 85 ## (with the same preconditioning steps) may be desirable in lieu of the | |
| 86 ## @iftex | |
| 87 ## @tex | |
| 88 ## Pad\'e | |
| 89 ## @end tex | |
| 90 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
91 ## @ifnottex |
| 8393 | 92 ## Pade' |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
93 ## @end ifnottex |
| 8393 | 94 ## approximation when |
| 95 ## @iftex | |
| 96 ## @tex | |
| 97 ## $D_q(a)$ | |
| 98 ## @end tex | |
| 99 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
100 ## @ifnottex |
| 8393 | 101 ## @code{Dq(a)} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
102 ## @end ifnottex |
| 8393 | 103 ## is ill-conditioned. |
| 104 ## @end deftypefn | |
| 105 | |
| 106 function r = expm (a) | |
| 107 | |
| 108 if (! ismatrix (a) || ! issquare (a)) | |
| 8664 | 109 error ("expm requires a square matrix"); |
| 8393 | 110 endif |
| 111 | |
| 112 n = rows (a); | |
| 8506 | 113 ## Trace reduction. |
| 8393 | 114 a(a == -Inf) = -realmax; |
| 115 trshift = trace (a) / length (a); | |
| 116 if (trshift > 0) | |
| 117 a -= trshift*eye (n); | |
| 118 endif | |
| 8506 | 119 ## Balancing. |
|
9048
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
120 [d, p, aa] = balance (a); |
| 8506 | 121 ## FIXME: can we both permute and scale at once? Or should we rather do |
| 122 ## this: | |
| 123 ## | |
|
9048
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
124 ## [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
|
125 ## [xx, p, aa] = balance (aa, "noscal"); |
| 8393 | 126 [f, e] = log2 (norm (aa, "inf")); |
| 127 s = max (0, e); | |
| 128 s = min (s, 1023); | |
| 129 aa *= 2^(-s); | |
| 130 | |
| 8506 | 131 ## Pade approximation for exp(A). |
| 8393 | 132 c = [5.0000000000000000e-1,... |
| 133 1.1666666666666667e-1,... | |
| 134 1.6666666666666667e-2,... | |
| 135 1.6025641025641026e-3,... | |
| 136 1.0683760683760684e-4,... | |
| 137 4.8562548562548563e-6,... | |
| 138 1.3875013875013875e-7,... | |
| 139 1.9270852604185938e-9]; | |
| 140 | |
| 141 a2 = aa^2; | |
| 142 id = eye (n); | |
| 143 x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 + id; | |
| 144 y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa; | |
| 145 | |
| 146 r = (x - y) \ (x + y); | |
| 147 | |
| 8506 | 148 ## Undo scaling by repeated squaring. |
| 149 for k = 1:s | |
| 8393 | 150 r ^= 2; |
| 151 endfor | |
| 152 | |
| 8506 | 153 ## inverse balancing. |
| 8757 | 154 d = diag (d); |
| 155 r = d * r / d; | |
| 8393 | 156 r = r(p, p); |
| 8506 | 157 ## Inverse trace reduction. |
| 8393 | 158 if (trshift >0) |
| 159 r *= exp (trshift); | |
| 160 endif | |
| 161 | |
| 162 endfunction |