Mercurial > octave-antonio
annotate scripts/linear-algebra/expm.m @ 20170:af2b7695f1c4 draft default tip @
gallery.m: clean and vectorize part of qmult.
* scripts/special-matrix/gallery.m: make qmult() follow Octave guidelines, vectorize last for loop in qmult().
| author | Antonio Pino Robles <data.script93@gmail.com> |
|---|---|
| date | Thu, 28 May 2015 18:32:47 +0200 |
| parents | 9fc020886ae9 |
| children | 03b9d17a2d95 |
| rev | line source |
|---|---|
|
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19040
diff
changeset
|
1 ## Copyright (C) 2008-2015 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 |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
29 ## |
| 8393 | 30 ## @example |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
31 ## expm (A) = I + A + A^2/2! + A^3/3! + @dots{} |
| 8393 | 32 ## @end example |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
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 |
|
19040
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
36 ## exponential; see @nospell{Moler and Van Loan}, @cite{Nineteen Dubious Ways |
|
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
37 ## to Compute the Exponential of a Matrix}, SIAM Review, 1978. This routine |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
38 ## uses Ward's diagonal Pad@'e approximation method with three step |
|
10791
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 |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
45 ## |
| 8393 | 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 |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
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. |
|
12584
7ef7e20057fa
Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
71 ## @seealso{logm, sqrtm} |
| 8393 | 72 ## @end deftypefn |
| 73 | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
74 function r = expm (A) |
| 8393 | 75 |
| 11477 | 76 if (nargin != 1) |
| 77 print_usage (); | |
| 78 endif | |
| 79 | |
|
19700
00e31f316a3a
Fix Matlab incompatibility of "ismatrix" (bug #42422).
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
19697
diff
changeset
|
80 if (! isnumeric (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
|
81 error ("expm: A must be a square matrix"); |
| 8393 | 82 endif |
| 83 | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
84 if (isscalar (A)) |
|
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
85 r = exp (A); |
|
17312
088d014a7fe2
Use semicolon after "return" statement in core m-files.
Rik <rik@octave.org>
parents:
14363
diff
changeset
|
86 return; |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
87 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
|
88 r = diag (exp (diag (A))); |
|
17312
088d014a7fe2
Use semicolon after "return" statement in core m-files.
Rik <rik@octave.org>
parents:
14363
diff
changeset
|
89 return; |
|
10831
1646bd8e3735
special case diagonal matrices and scalars in expm
Jaroslav Hajek <highegg@gmail.com>
parents:
10791
diff
changeset
|
90 endif |
|
1646bd8e3735
special case diagonal matrices and scalars in expm
Jaroslav Hajek <highegg@gmail.com>
parents:
10791
diff
changeset
|
91 |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
92 n = rows (A); |
| 8506 | 93 ## Trace reduction. |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
94 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
|
95 trshift = trace (A) / length (A); |
| 8393 | 96 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
|
97 A -= trshift*eye (n); |
| 8393 | 98 endif |
| 8506 | 99 ## Balancing. |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
100 [d, p, aa] = balance (A); |
| 8506 | 101 ## FIXME: can we both permute and scale at once? Or should we rather do |
| 102 ## this: | |
| 103 ## | |
|
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
10831
diff
changeset
|
104 ## [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
|
105 ## [xx, p, aa] = balance (aa, "noscal"); |
| 8393 | 106 [f, e] = log2 (norm (aa, "inf")); |
| 107 s = max (0, e); | |
| 108 s = min (s, 1023); | |
| 109 aa *= 2^(-s); | |
| 110 | |
| 8506 | 111 ## Pade approximation for exp(A). |
| 8393 | 112 c = [5.0000000000000000e-1,... |
| 113 1.1666666666666667e-1,... | |
| 114 1.6666666666666667e-2,... | |
| 115 1.6025641025641026e-3,... | |
| 116 1.0683760683760684e-4,... | |
| 117 4.8562548562548563e-6,... | |
| 118 1.3875013875013875e-7,... | |
| 119 1.9270852604185938e-9]; | |
| 120 | |
| 121 a2 = aa^2; | |
| 122 id = eye (n); | |
| 123 x = (((c(8) * a2 + c(6) * id) * a2 + c(4) * id) * a2 + c(2) * id) * a2 + id; | |
| 124 y = (((c(7) * a2 + c(5) * id) * a2 + c(3) * id) * a2 + c(1) * id) * aa; | |
| 125 | |
| 126 r = (x - y) \ (x + y); | |
| 127 | |
| 8506 | 128 ## Undo scaling by repeated squaring. |
| 129 for k = 1:s | |
| 8393 | 130 r ^= 2; |
| 131 endfor | |
| 132 | |
| 8506 | 133 ## inverse balancing. |
| 8757 | 134 d = diag (d); |
| 135 r = d * r / d; | |
| 9819 | 136 r(p, p) = r; |
| 8506 | 137 ## Inverse trace reduction. |
| 8393 | 138 if (trshift >0) |
| 139 r *= exp (trshift); | |
| 140 endif | |
| 141 | |
| 142 endfunction | |
| 11477 | 143 |
| 144 | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
145 %!assert (norm (expm ([1 -1;0 1]) - [e -e; 0 e]) < 1e-5); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
146 %!assert (expm ([1 -1 -1;0 1 -1; 0 0 1]), [e -e -e/2; 0 e -e; 0 0 e], 1e-5); |
| 11477 | 147 |
| 148 %!assert (expm (10), expm (10)) | |
| 149 %!assert (full (expm (eye (3))), expm (full (eye (3)))) | |
| 150 %!assert (full (expm (10*eye (3))), expm (full (10*eye (3))), 8*eps) | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
151 |
|
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19700
diff
changeset
|
152 ## Test input validation |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
153 %!error expm () |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
154 %!error expm (1, 2) |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
155 %!error <expm: A must be a square matrix> expm ([1 0;0 1; 2 2]) |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
156 |