Mercurial > hg > octave-lyh
annotate scripts/linear-algebra/expm.m @ 11477:a02d00dd3d5f
expm.m: new tests
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 10 Jan 2011 14:50:33 -0500 |
| parents | 1740012184f9 |
| children | fd0a3ac60b0e |
| 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 -*- |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
20 ## @deftypefn {Function File} {} expm (@var{A}) |
| 8393 | 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 | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
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 | |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
38 ## uses Ward's diagonal Pad@'e approximation method with three step |
|
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
39 ## preconditioning (SIAM Journal on Numerical Analysis, 1977). Diagonal |
|
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
40 ## Pad@'e approximations are rational polynomials of matrices |
| 8393 | 41 ## @tex |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
42 ## $D_q(A)^{-1}N_q(A)$ |
| 8393 | 43 ## @end tex |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
44 ## @ifnottex |
| 8393 | 45 ## |
| 46 ## @example | |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
47 ## @group |
| 8393 | 48 ## -1 |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
49 ## D (A) N (A) |
|
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9048
diff
changeset
|
50 ## @end group |
| 8393 | 51 ## @end example |
| 52 ## | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
53 ## @end ifnottex |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
54 ## whose Taylor series matches the first |
| 8393 | 55 ## @tex |
| 56 ## $2 q + 1 $ | |
| 57 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
58 ## @ifnottex |
| 8393 | 59 ## @code{2q+1} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
60 ## @end ifnottex |
| 8393 | 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 | |
|
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10635
diff
changeset
|
63 ## Pad@'e approximation when |
| 8393 | 64 ## @tex |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
65 ## $D_q(A)$ |
| 8393 | 66 ## @end tex |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
67 ## @ifnottex |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
68 ## @code{Dq(A)} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8506
diff
changeset
|
69 ## @end ifnottex |
| 8393 | 70 ## is ill-conditioned. |
| 71 ## @end deftypefn | |
| 72 | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
73 function r = expm (A) |
| 8393 | 74 |
| 11477 | 75 if (nargin != 1) |
| 76 print_usage (); | |
| 77 endif | |
| 78 | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
79 if (! ismatrix (A) || ! issquare (A)) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11471
diff
changeset
|
80 error ("expm: A must be a square matrix"); |
| 8393 | 81 endif |
| 82 | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
83 if (isscalar (A)) |
|
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
84 r = exp (A); |
|
10831
1646bd8e3735
special case diagonal matrices and scalars in expm
Jaroslav Hajek <highegg@gmail.com>
parents:
10791
diff
changeset
|
85 return |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
86 elseif (strfind (typeinfo (A), "diagonal matrix")) |
|
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
87 r = diag (exp (diag (A))); |
|
10831
1646bd8e3735
special case diagonal matrices and scalars in expm
Jaroslav Hajek <highegg@gmail.com>
parents:
10791
diff
changeset
|
88 return |
|
1646bd8e3735
special case diagonal matrices and scalars in expm
Jaroslav Hajek <highegg@gmail.com>
parents:
10791
diff
changeset
|
89 endif |
|
1646bd8e3735
special case diagonal matrices and scalars in expm
Jaroslav Hajek <highegg@gmail.com>
parents:
10791
diff
changeset
|
90 |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
91 n = rows (A); |
| 8506 | 92 ## Trace reduction. |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
93 A(A == -Inf) = -realmax; |
|
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
94 trshift = trace (A) / length (A); |
| 8393 | 95 if (trshift > 0) |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
96 A -= trshift*eye (n); |
| 8393 | 97 endif |
| 8506 | 98 ## Balancing. |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
99 [d, p, aa] = balance (A); |
| 8506 | 100 ## FIXME: can we both permute and scale at once? Or should we rather do |
| 101 ## this: | |
| 102 ## | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
103 ## [d, xx, aa] = balance (A, "noperm"); |
|
9048
867d5d1aed06
swap out args in balance for M*b compat
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
104 ## [xx, p, aa] = balance (aa, "noscal"); |
| 8393 | 105 [f, e] = log2 (norm (aa, "inf")); |
| 106 s = max (0, e); | |
| 107 s = min (s, 1023); | |
| 108 aa *= 2^(-s); | |
| 109 | |
| 8506 | 110 ## Pade approximation for exp(A). |
| 8393 | 111 c = [5.0000000000000000e-1,... |
| 112 1.1666666666666667e-1,... | |
| 113 1.6666666666666667e-2,... | |
| 114 1.6025641025641026e-3,... | |
| 115 1.0683760683760684e-4,... | |
| 116 4.8562548562548563e-6,... | |
| 117 1.3875013875013875e-7,... | |
| 118 1.9270852604185938e-9]; | |
| 119 | |
| 120 a2 = aa^2; | |
| 121 id = eye (n); | |
| 122 x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 + id; | |
| 123 y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa; | |
| 124 | |
| 125 r = (x - y) \ (x + y); | |
| 126 | |
| 8506 | 127 ## Undo scaling by repeated squaring. |
| 128 for k = 1:s | |
| 8393 | 129 r ^= 2; |
| 130 endfor | |
| 131 | |
| 8506 | 132 ## inverse balancing. |
| 8757 | 133 d = diag (d); |
| 134 r = d * r / d; | |
| 9819 | 135 r(p, p) = r; |
| 8506 | 136 ## Inverse trace reduction. |
| 8393 | 137 if (trshift >0) |
| 138 r *= exp (trshift); | |
| 139 endif | |
| 140 | |
| 141 endfunction | |
| 11477 | 142 |
| 143 %!assert(norm(expm([1 -1;0 1]) - [e -e; 0 e]) < 1e-5); | |
| 144 %!assert(expm([1 -1 -1;0 1 -1; 0 0 1]), [e -e -e/2; 0 e -e; 0 0 e], 1e-5); | |
| 145 | |
| 146 %% Test input validation | |
| 147 %!error expm (); | |
| 148 %!error expm (1, 2); | |
| 149 %!error <expm: A must be a square matrix> expm([1 0;0 1; 2 2]); | |
| 150 | |
| 151 %!assert (expm (10), expm (10)) | |
| 152 %!assert (full (expm (eye (3))), expm (full (eye (3)))) | |
| 153 %!assert (full (expm (10*eye (3))), expm (full (10*eye (3))), 8*eps) |