Mercurial > matrix-functions

comparison funm_files/logm_isst.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function [X, iter] = logm_isst(T, prnt)
2 %LOGM_ISST Log of triangular matrix by Schur-Pade method with scaling.
3 % X = LOGM_ISST(A) computes the logarithm of an upper triangular
4 % matrix A, for a matrix with no nonpositive real eigenvalues,
5 % using the inverse scaling and squaring method with Pade
6 % approximation. TOL is an error tolerance.
7 % [X, ITER] = LOGM_ISST(A, PRNT) returns the number ITER of square
8 % roots computed and prints this information if PRNT is nonzero.
9
10 % References:
11 % S. H. Cheng, N. J. Higham, C. S. Kenney, and A. J. Laub, Approximating the
12 % logarithm of a matrix to specified accuracy, SIAM J. Matrix Anal. Appl.,
13 % 22(4):1112-1125, 2001.
14 % N. J. Higham, Evaluating Pade approximants of the matrix logarithm,
15 % SIAM J. Matrix Anal. Appl., 22(4):1126-1135, 2001.
16
17 if nargin < 2, prnt = 0; end
18 n = length(T);
19
20 if any( imag(diag(T)) == 0 & real(diag(T)) <= 0 )
21 error('A must not have nonpositive real eigenvalues!')
22 end
23
24 if n == 1, X = log(T); iter = 0; return, end
25
26 R = T;
27 maxlogiter = 50;
28
29 for iter = 0:maxlogiter
30
31 phi = norm(T-eye(n),'fro');
32
33 if phi <= 0.25
34 if prnt, fprintf('LOGM_ISST computed %g square roots. \n', iter), end
35 break
36 end
37 if iter == maxlogiter, error('Too many square roots in LOGM_ISST.\n'), end
38
39 % Compute upper triangular square root R of T, a column at a time.
40 for j=1:n
41 R(j,j) = sqrt(T(j,j));
42 for i=j-1:-1:1
43 R(i,j) = (T(i,j) - R(i,i+1:j-1)*R(i+1:j-1,j))/(R(i,i) + R(j,j));
44 end
45 end
46 T = R;
47 end
48
49 X = 2^(iter)*logm_pf(T-eye(n),8);
50
51 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 function S = logm_pf(A,m)
53 %LOGM_PF Pade approximation to matrix log by partial fraction expansion.
54 % Y = LOGM_PF(A,m) approximates LOG(I+A).
55
56 [nodes,wts] = gauss_legendre(m);
57 % Convert from [-1,1] to [0,1].
58 nodes = (nodes + 1)/2;
59 wts = wts/2;
60
61 n = length(A);
62 S = zeros(n);
63
64 for j=1:m
65 S = S + wts(j)*(A/(eye(n) + nodes(j)*A));
66 end
67
68 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69 function [x,w] = gauss_legendre(n)
70 %GAUSS_LEGENDRE Nodes and weights for Gauss-Legendre quadrature.
71
72 % Reference:
73 % G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature
74 % rules, Math. Comp., 23(106):221-230, 1969.
75
76 i = 1:n-1;
77 v = i./sqrt((2*i).^2-1);
78 [V,D] = eig( diag(v,-1)+diag(v,1) );
79 x = diag(D);
80 w = 2*(V(1,:)'.^2);