Mercurial > hg > octave-lyh

annotate scripts/linear-algebra/logm.m @ 11471:994e2a93a8e2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use uppercase 'A' to refer to matrix inputs in m-files.
author Rik <octave@nomad.inbox5.com>
date Sun, 09 Jan 2011 16:01:05 -0800
parents fb98284fcc20
children 1740012184f9
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
4 ##
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
5 ## This file is part of Octave.
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
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
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
19 ## <http://www.gnu.org/licenses/>.
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
20
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
21 ## -*- texinfo -*-
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
22 ## @deftypefn {Function File} {@var{s} =} logm (@var{A})
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
23 ## @deftypefnx {Function File} {@var{s} =} logm (@var{A}, @var{opt_iters})
10827
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{})
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
25 ## Compute the matrix logarithm of the square matrix @var{A}. The
10827
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
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
29 ## logm(@var{A}) = 2^k * logm(@var{A}^(1 / 2^k))
10827
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
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
36 ## @end deftypefn
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
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
11033
d9c8916bb9dd Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents: 10827
diff changeset
39 ## (SIAM, 2008.)
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
40 ##
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
41
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
42 function [s, iters] = logm (A, opt_iters = 100)
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
43
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
44 if (nargin == 0 || nargin > 2)
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5307
diff changeset
45 print_usage ();
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
46 endif
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
47
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
48 if (! issquare (A))
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
49 error ("logm: argument must be a square matrix");
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
50 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
51
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
52 [u, s] = schur (A);
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
53
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
54 if (isreal (A))
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
55 [u, s] = rsf2csf (u, s);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
56 endif
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 (any (diag (s) < 0))
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
59 warning ("Octave:logm:non-principal",
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
60 "logm: principal matrix logarithm is not defined for matrices with negative eigenvalues; computing non-principal logarithm");
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
61 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
62
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
63 k = 0;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
64 ## 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
65 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
66 p = 0;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
67 m = 7;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
68 while (k < opt_iters)
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
69 tau = norm (s - eye (size (s)),1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
70 if (tau <= theta (7))
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
71 p = p + 1;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
72 j(1) = find (tau <= theta, 1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
73 j(2) = find (tau / 2 <= theta, 1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
74 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
75 m = j(1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
76 break
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
77 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
78 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
79 k = k + 1;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
80 s = sqrtm (s);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
81 endwhile
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
82
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
83 if (k >= opt_iters)
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
84 warning ("logm: maximum number of square roots exceeded; results may still be accurate");
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
85 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
86
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
87 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
88
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
89 s = 2^k * u * s * u';
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
90
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
91 if (nargout == 2)
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
92 iters = k;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
93 endif
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
94
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
95 endfunction
10825
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 ################## ANCILLARY FUNCTIONS ################################
11033
d9c8916bb9dd Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents: 10827
diff changeset
98 ###### Taken from the mfttoolbox (GPL 3) by D. Higham.
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
99 ###### Reference:
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
100 ###### D. Higham, Functions of Matrices: Theory and Computation
11033
d9c8916bb9dd Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents: 10827
diff changeset
101 ###### (SIAM, 2008.).
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
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
104 ##LOGM_PADE_PF Evaluate Pade approximant to matrix log by partial fractions.
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
105 ## Y = LOGM_PADE_PF(A,M) evaluates the [M/M] Pade approximation to
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
106 ## LOG(EYE(SIZE(A))+A) using a partial fraction expansion.
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
107
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
108 function s = logm_pade_pf (A, m)
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
109 [nodes, wts] = gauss_legendre (m);
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
110 ## Convert from [-1,1] to [0,1].
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
111 nodes = (nodes+1)/2;
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
112 wts = wts/2;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
113
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
114 n = length (A);
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
115 s = zeros (n);
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
116 for j = 1:m
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
117 s += wts(j)*(A/(eye (n) + nodes(j)*A));
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
118 endfor
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
119 endfunction
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
120
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
121 ######################################################################
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
122 ## 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
123 ## [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
124 ## for N-point Gauss-Legendre quadrature.
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
125
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
126 ## Reference:
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
127 ## 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
128 ## 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
129
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
130 function [x, w] = gauss_legendre (n)
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
131 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
132 v = i./sqrt ((2*i).^2-1);
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
133 [V, D] = eig (diag (v, -1) + diag (v, 1));
10827
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
134 x = diag (D);
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
135 w = 2*(V(1,:)'.^2);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
136 endfunction
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
137
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
138
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
139 %!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
140 %!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
141 %!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
142 %!
10827
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
143 %% Test input validation
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
144 %!error logm ();
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
145 %!error logm (1, 2, 3);
11465
fb98284fcc20 logm: fix test
John W. Eaton <jwe@octave.org>
parents: 11458
diff changeset
146 %!error <logm: argument must be a square matrix> logm([1 0;0 1; 2 2]);