Mercurial > octave
annotate scripts/linear-algebra/logm.m @ 10827:228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Wed, 28 Jul 2010 11:57:39 -0700 |
parents | cace99cb01ab |
children | d9c8916bb9dd |
rev | line source |
---|---|
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
1 ## Copyright (C) 2010 Richard T. Guy <guyrt7@wfu.edu> |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
2 ## Copyright (C) 2010 Marco Caliari <marco.caliari@univr.it> |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
3 ## Copyright (C) 2008 N.J. Higham |
4334 | 4 ## |
5 ## This file is part of Octave. | |
6 ## | |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
7 ## Octave is free software; you can redistribute it and/or modify it |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
8 ## under the terms of the GNU General Public License as published by |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
9 ## the Free Software Foundation; either version 3 of the License, or (at |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
10 ## your option) any later version. |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
11 ## |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
12 ## Octave is distributed in the hope that it will be useful, but |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
15 ## General Public License for more details. |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
16 ## |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
17 ## You should have received a copy of the GNU General Public License |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
18 ## along with Octave; see the file COPYING. If not, see |
7016 | 19 ## <http://www.gnu.org/licenses/>. |
4334 | 20 |
21 ## -*- texinfo -*- | |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
22 ## @deftypefn {Function File} {@var{s} =} logm (@var{a}) |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
23 ## @deftypefnx {Function File} {@var{s} =} logm (@var{a}, @var{opt_iters}) |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
24 ## @deftypefnx {Function File} {[@var{s}, @var{iters}] =} logm (@dots{}) |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
25 ## Compute the matrix logarithm of the square matrix @var{a}. The |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
26 ## implementation utilizes a Pad@'e approximant and the identity |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
27 ## |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
28 ## @example |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
29 ## logm(@var{a}) = 2^k * logm(@var{a}^(1 / 2^k)) |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
30 ## @end example |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
31 ## |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
32 ## The optional argument @var{opt_iters} is the maximum number of square roots |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
33 ## to compute and defaults to 100. The optional output @var{iters} is the |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
34 ## number of square roots actually computed. |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
35 ## |
4334 | 36 ## @end deftypefn |
37 | |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
38 ## Reference: N. J. Higham, Functions of Matrices: Theory and Computation |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
39 ## (SIAM, 2008.) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
40 ## |
4334 | 41 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
42 function [s, iters] = logm (a, opt_iters) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
43 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
44 if (nargin == 0) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
45 print_usage (); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
46 elseif (nargin < 2) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
47 opt_iters = 100; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
48 elseif (nargin > 2) |
6046 | 49 print_usage (); |
4334 | 50 endif |
51 | |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
52 if (! issquare (a)) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
53 error ("logm: argument must be a square matrix."); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
54 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
55 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
56 [u, s] = schur (a); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
57 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
58 if (isreal (a)) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
59 [u, s] = rsf2csf (u, s); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
60 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
61 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
62 if (any (diag (s) < 0)) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
63 warning ("Octave:logm:non-principal", |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
64 ["logm: Matrix has negative eigenvalues.", ... |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
65 " Principal matrix logarithm is not defined.", ... |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
66 " Computing non-principal logarithm."]); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
67 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
68 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
69 k = 0; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
70 ## Algorithm 11.9 in "Function of matrices", by N. Higham |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
71 theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1]; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
72 p = 0; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
73 m = 7; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
74 while (k < opt_iters) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
75 tau = norm (s - eye (size (s)),1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
76 if (tau <= theta (7)) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
77 p = p + 1; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
78 j(1) = find (tau <= theta, 1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
79 j(2) = find (tau / 2 <= theta, 1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
80 if (j(1) - j(2) <= 1 || p == 2) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
81 m = j(1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
82 break |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
83 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
84 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
85 k = k + 1; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
86 s = sqrtm (s); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
87 endwhile |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
88 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
89 if (k >= opt_iters) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
90 warning ("logm: Maximum number of square roots exceeded. Results may still be accurate."); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
91 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
92 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
93 s = logm_pade_pf (s - eye (size (s)), m); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
94 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
95 s = 2^k * u * s * u'; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
96 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
97 if (nargout == 2) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
98 iters = k; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
99 endif |
4334 | 100 |
101 endfunction | |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
102 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
103 ################## ANCILLARY FUNCTIONS ################################ |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
104 ###### Taken from the mfttoolbox (GPL 3) by D. Higham. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
105 ###### Reference: |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
106 ###### D. Higham, Functions of Matrices: Theory and Computation |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
107 ###### (SIAM, 2008.). |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
108 ####################################################################### |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
109 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
110 ##LOGM_PADE_PF Evaluate Pade approximant to matrix log by partial fractions. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
111 ## Y = LOGM_PADE_PF(a,M) evaluates the [M/M] Pade approximation to |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
112 ## LOG(EYE(SIZE(a))+a) using a partial fraction expansion. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
113 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
114 function s = logm_pade_pf(a,m) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
115 [nodes,wts] = gauss_legendre(m); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
116 ## Convert from [-1,1] to [0,1]. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
117 nodes = (nodes + 1)/2; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
118 wts = wts/2; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
119 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
120 n = length(a); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
121 s = zeros(n); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
122 for j=1:m |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
123 s = s + wts(j)*(a/(eye(n) + nodes(j)*a)); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
124 endfor |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
125 endfunction |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
126 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
127 ###################################################################### |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
128 ## GAUSS_LEGENDRE Nodes and weights for Gauss-Legendre quadrature. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
129 ## [X,W] = GAUSS_LEGENDRE(N) computes the nodes X and weights W |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
130 ## for N-point Gauss-Legendre quadrature. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
131 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
132 ## Reference: |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
133 ## G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
134 ## rules, Math. Comp., 23(106):221-230, 1969. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
135 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
136 function [x,w] = gauss_legendre(n) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
137 i = 1:n-1; |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
138 v = i./sqrt ((2*i).^2-1); |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
139 [V,D] = eig ( diag(v,-1)+diag(v,1) ); |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
140 x = diag (D); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
141 w = 2*(V(1,:)'.^2); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
142 endfunction |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
143 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
144 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
145 %!assert(norm(logm([1 -1;0 1]) - [0 -1; 0 0]) < 1e-5); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
146 %!assert(norm(expm(logm([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5); |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
147 %!assert(logm([1 -1 -1;0 1 -1; 0 0 1]), [0 -1 -1.5; 0 0 -1; 0 0 0], 1e-5); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
148 %! |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
149 %% Test input validation |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
150 %!error logm (); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
151 %!error logm (1, 2, 3); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
152 %!error <logm: argument must be a square matrix.> logm([1 0;0 1; 2 2]); |