Mercurial > matrix-functions
view funm_files/funm_atom.m @ 0:8f23314345f4 draft
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 |
line wrap: on
line source
function [F,n_terms] = fun_atom(T,fun,tol,prnt) %FUNM_ATOM Function of triangular matrix with nearly constant diagonal. % [F, N_TERMS] = FUNM_ATOM(T, FUN, TOL, PRNT) evaluates function % FUN at the upper triangular matrix T, where T has nearly constant % diagonal. A Taylor series is used. % FUN(X,K) must return the K'th derivative of % the function represented by FUN evaluated at the vector X. % TOL is a convergence tolerance for the Taylor series, % defaulting to EPS. % If PRNT ~= 0 trace information is printed. % N_TERMS is the number of terms taken in the Taylor series. % N_TERMS = -1 signals lack of convergence. if nargin < 3 | isempty(tol), tol = eps; end if nargin < 4, prnt = 0; end if isequal(fun,@fun_log) % LOG is special case. [F,n_terms] = logm_isst(T,prnt); return end itmax = 500; n = length(T); if n == 1, F = feval(fun,T,0); n_terms = 1; return, end lambda = sum(diag(T))/n; F = eye(n)*feval(fun,lambda,0); f_deriv_max = zeros(itmax+n-1,1); N = T - lambda*eye(n); mu = norm( (eye(n)-abs(triu(T,1)))\ones(n,1),inf ); P = N; max_d = 1; for k = 1:itmax f = feval(fun,lambda,k); F_old = F; F = F + P*f; rel_diff = norm(F - F_old,inf)/(tol+norm(F_old,inf)); if prnt fprintf('%3.0f: coef = %5.0e', k, abs(f)/factorial(k)); fprintf(' N^k/k! = %7.1e', norm(P,inf)); fprintf(' rel_d = %5.0e',rel_diff); fprintf(' abs_d = %5.0e',norm(F - F_old,inf)); end P = P*N/(k+1); if rel_diff <= tol % Approximate the maximum of derivatives in convex set containing % eigenvalues by maximum of derivatives at eigenvalues. for j = max_d:k+n-1 f_deriv_max(j) = norm(feval(fun,diag(T),j),inf); end max_d = k+n; omega = 0; for j = 0:n-1 omega = max(omega,f_deriv_max(k+j)/factorial(j)); end trunc = norm(P,inf)*mu*omega; % norm(F) moved to RHS to avoid / 0. if prnt fprintf(' [trunc,test] = [%5.0e %5.0e]', ... trunc, tol*norm(F,inf)) end if prnt == 5, trunc = 0; end % Force simple stopping test. if trunc <= tol*norm(F,inf) n_terms = k; if prnt, fprintf('\n'), end, return end end if prnt, fprintf('\n'), end end n_terms = -1;