Mercurial > matrix-functions
comparison funm_files/funm_atom.m @ 8:a587712dcf5f draft default tip
funm_atom.m: rename fun_atom to funm_atom
* funm_atom.m: rename fun_atom to funm_atom.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Fri, 29 May 2015 09:48:36 +0200 |
parents | funm_files/fun_atom.m@e0817ccb2834 |
children |
comparison
equal
deleted
inserted
replaced
7:e0817ccb2834 | 8:a587712dcf5f |
---|---|
1 function [F,n_terms] = fun_atom(T,fun,tol,prnt) | |
2 %FUNM_ATOM Function of triangular matrix with nearly constant diagonal. | |
3 % [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function | |
4 % FUN at the upper triangular matrix T, where T has nearly constant | |
5 % diagonal. A Taylor series is used. | |
6 % FUN(X,K) must return the K'th derivative of | |
7 % the function represented by FUN evaluated at the vector X. | |
8 % TOL is a convergence tolerance for the Taylor series, | |
9 % defaulting to EPS. | |
10 % If PRNT ~= 0 trace information is printed. | |
11 % N_TERMS is the number of terms taken in the Taylor series. | |
12 % N_TERMS = -1 signals lack of convergence. | |
13 | |
14 if nargin < 3 | isempty(tol), tol = eps; end | |
15 if nargin < 4, prnt = 0; end | |
16 | |
17 if isequal(fun,@fun_log) % LOG is special case. | |
18 [F,n_terms] = logm_isst(T,prnt); | |
19 return | |
20 end | |
21 | |
22 itmax = 500; | |
23 | |
24 n = length(T); | |
25 if n == 1, F = feval(fun,T,0); n_terms = 1; return, end | |
26 | |
27 lambda = sum(diag(T))/n; | |
28 F = eye(n)*feval(fun,lambda,0); | |
29 f_deriv_max = zeros(itmax+n-1,1); | |
30 N = T - lambda*eye(n); | |
31 mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf ); | |
32 | |
33 P = N; | |
34 max_d = 1; | |
35 | |
36 for k = 1:itmax | |
37 | |
38 f = feval(fun,lambda,k); | |
39 F_old = F; | |
40 F = F + P*f; | |
41 rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf)); | |
42 if prnt | |
43 fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k)); | |
44 fprintf(' N^k/k! = %7.1e', norm(P,inf)); | |
45 fprintf(' rel_d = %5.0e',rel_diff); | |
46 fprintf(' abs_d = %5.0e',norm(F - F_old,inf)); | |
47 end | |
48 P = P*N/(k+1); | |
49 | |
50 if rel_diff <= tol | |
51 | |
52 % Approximate the maximum of derivatives in convex set containing | |
53 % eigenvalues by maximum of derivatives at eigenvalues. | |
54 for j = max_d:k+n-1 | |
55 f_deriv_max(j) = norm(feval(fun,diag(T),j),inf); | |
56 end | |
57 max_d = k+n; | |
58 omega = 0; | |
59 for j = 0:n-1 | |
60 omega = max(omega,f_deriv_max(k+j)/factorial(j)); | |
61 end | |
62 | |
63 trunc = norm(P,inf)*mu*omega; % norm(F) moved to RHS to avoid / 0. | |
64 if prnt | |
65 fprintf(' [trunc,test] = [%5.0e %5.0e]', ... | |
66 trunc, tol*norm(F,inf)) | |
67 end | |
68 if prnt == 5, trunc = 0; end % Force simple stopping test. | |
69 if trunc <= tol*norm(F,inf) | |
70 n_terms = k; | |
71 if prnt, fprintf('\n'), end, return | |
72 end | |
73 end | |
74 | |
75 if prnt, fprintf('\n'), end | |
76 | |
77 end | |
78 n_terms = -1; |