Mercurial > hg > octave-lyh

comparison scripts/linear-algebra/logm.m @ 11458:93a039fe681e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
logm: style fixes
author John W. Eaton <jwe@octave.org>
date Fri, 07 Jan 2011 14:14:35 -0500
parents d9c8916bb9dd
children fb98284fcc20
comparison
equal deleted inserted replaced
11457:33f6384d2b78 11458:93a039fe681e
37 37
38 ## Reference: N. J. Higham, Functions of Matrices: Theory and Computation 38 ## Reference: N. J. Higham, Functions of Matrices: Theory and Computation
39 ## (SIAM, 2008.) 39 ## (SIAM, 2008.)
40 ## 40 ##
41 41
42 function [s, iters] = logm (a, opt_iters) 42 function [s, iters] = logm (a, opt_iters = 100)
43 43
44 if (nargin == 0) 44 if (nargin == 0 || nargin > 2)
45 print_usage ();
46 elseif (nargin < 2)
47 opt_iters = 100;
48 elseif (nargin > 2)
49 print_usage (); 45 print_usage ();
50 endif 46 endif
51 47
52 if (! issquare (a)) 48 if (! issquare (a))
53 error ("logm: argument must be a square matrix."); 49 error ("logm: argument must be a square matrix");
54 endif 50 endif
55 51
56 [u, s] = schur (a); 52 [u, s] = schur (a);
57 53
58 if (isreal (a)) 54 if (isreal (a))
59 [u, s] = rsf2csf (u, s); 55 [u, s] = rsf2csf (u, s);
60 endif 56 endif
61 57
62 if (any (diag (s) < 0)) 58 if (any (diag (s) < 0))
63 warning ("Octave:logm:non-principal", 59 warning ("Octave:logm:non-principal",
64 ["logm: Matrix has negative eigenvalues.", ... 60 "logm: principal matrix logarithm is not defined for matrices with negative eigenvalues; computing non-principal logarithm");
65 " Principal matrix logarithm is not defined.", ...
66 " Computing non-principal logarithm."]);
67 endif 61 endif
68 62
69 k = 0; 63 k = 0;
70 ## Algorithm 11.9 in "Function of matrices", by N. Higham 64 ## Algorithm 11.9 in "Function of matrices", by N. Higham
71 theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1]; 65 theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1];
85 k = k + 1; 79 k = k + 1;
86 s = sqrtm (s); 80 s = sqrtm (s);
87 endwhile 81 endwhile
88 82
89 if (k >= opt_iters) 83 if (k >= opt_iters)
90 warning ("logm: Maximum number of square roots exceeded. Results may still be accurate."); 84 warning ("logm: maximum number of square roots exceeded; results may still be accurate");
91 endif 85 endif
92 86
93 s = logm_pade_pf (s - eye (size (s)), m); 87 s = logm_pade_pf (s - eye (size (s)), m);
94 88
95 s = 2^k * u * s * u'; 89 s = 2^k * u * s * u';
109 103
110 ##LOGM_PADE_PF Evaluate Pade approximant to matrix log by partial fractions. 104 ##LOGM_PADE_PF Evaluate Pade approximant to matrix log by partial fractions.
111 ## Y = LOGM_PADE_PF(a,M) evaluates the [M/M] Pade approximation to 105 ## Y = LOGM_PADE_PF(a,M) evaluates the [M/M] Pade approximation to
112 ## LOG(EYE(SIZE(a))+a) using a partial fraction expansion. 106 ## LOG(EYE(SIZE(a))+a) using a partial fraction expansion.
113 107
114 function s = logm_pade_pf(a,m) 108 function s = logm_pade_pf (a, m)
115 [nodes,wts] = gauss_legendre(m); 109 [nodes, wts] = gauss_legendre (m);
116 ## Convert from [-1,1] to [0,1]. 110 ## Convert from [-1,1] to [0,1].
117 nodes = (nodes + 1)/2; 111 nodes = (nodes+1)/2;
118 wts = wts/2; 112 wts = wts/2;
119 113
120 n = length(a); 114 n = length (a);
121 s = zeros(n); 115 s = zeros (n);
122 for j=1:m 116 for j = 1:m
123 s = s + wts(j)*(a/(eye(n) + nodes(j)*a)); 117 s += wts(j)*(a/(eye (n) + nodes(j)*a));
124 endfor 118 endfor
125 endfunction 119 endfunction
126 120
127 ###################################################################### 121 ######################################################################
128 ## GAUSS_LEGENDRE Nodes and weights for Gauss-Legendre quadrature. 122 ## GAUSS_LEGENDRE Nodes and weights for Gauss-Legendre quadrature.
131 125
132 ## Reference: 126 ## Reference:
133 ## G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature 127 ## G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature
134 ## rules, Math. Comp., 23(106):221-230, 1969. 128 ## rules, Math. Comp., 23(106):221-230, 1969.
135 129
136 function [x,w] = gauss_legendre(n) 130 function [x, w] = gauss_legendre (n)
137 i = 1:n-1; 131 i = 1:n-1;
138 v = i./sqrt ((2*i).^2-1); 132 v = i./sqrt ((2*i).^2-1);
139 [V,D] = eig ( diag(v,-1)+diag(v,1) ); 133 [V, D] = eig (diag (v, -1) + diag (v, 1));
140 x = diag (D); 134 x = diag (D);
141 w = 2*(V(1,:)'.^2); 135 w = 2*(V(1,:)'.^2);
142 endfunction 136 endfunction
143 137
144 138