Jaroslav,
Since you are listed as maintainer of the octave-forge linear-algebra
package:
Would you mind if I commit the attached updated funm and thfm functions?
Long explanation below:
Inspired by the recent discussion about matrix functions and funm.m (see
thread here:
https://mailman.cae.wisc.edu/pipermail/help-octave/2010-December/043907.html),
and because I could use a better funm.m myself, I've looked a bit into:
(1) augmenting the current, primitive (some say: naive) funm.m in
octave-forge that I contributed 11 years ago and that has been untouched
since then, with some very basic improvements;
(2) incorporating N. Higham's funm stuff from
http://www.maths.manchester.ac.uk/~higham/NAMF/ (funm_files.tar.gz).
(3) alternatively, incorporating (parts of) N. Higham's Matrix Function
Toolbox (GPL v3) at
http://www.maths.manchester.ac.uk/~higham/mftoolbox/
especially funm_simple.m
DISCLAIMER: I'm not a mathematician and I'm no expert in linear algebra.
Yet funm.m is fairly vital for some of my applications.
Anyway, as to (1):
------------------
The existing funm.m can be improved by:
- Allowing function handles as input arg
- Checking whether the function-to-be-applied is one included in
octave-forge's thfm.m and if so, invoke thfm() (which resorts to expm()
for some functions); except for logm which has recently been improved in
octave-core.
The "trick" suggested by Davies (2007), as indicated by Higham (i.e.,
adding noise, see bottom of
https://mailman.cae.wisc.edu/pipermail/help-octave/2010-December/044012.html)
is easy to incorporate too.
I've attached the resulting funm.m for review. An improved thfm.m is
attached too as the old one called funm in turn (which would imply an
endless loop).
If you don't object, I'll commit this version to octave-forge next or
following week, as a first update and place holder until (hopefully) a
much better version becomes available.
I know funm.m is still ultra-basic, but it's better than what's
presently available.
As to (2):
----------
As far as I understand, Higham's (and current Matlab's) funm.m can only
invoke functions that have Taylor expansion, or an easily coded series
of successive/recursive derivatives. Helper function files built for
this have to be supplied by the user. TMW site has an example using
(exp() + sin()).
In my case I invoke a.o. (modified) Bessel functions for which simple
relations exist between successive derivatives (and they can easily be
written as Taylor expansions anyway).
But in case of other functions, Higham's funm version may pose a problem.
Moreover, I have the impression that for small, diagonizable matrices,
Higham's (and Matlab's) funm.m may be unduly demanding and expensive.
As to (3):
----------
The funm_simple.m mentioned in there is only marginally better than the
old one in octave-forge.
Philip
%USAGE y = thfm ( x, MODE )
%
% trigonometric/hyperbolic functions of square matrix x
%
%MODE cos sin tan sec csc cot
% cosh sinh tanh sech csch coth
% acos asin atan asec acsc acot
% acosh asinh atanh asech acsch acoth
% sqrt log exp
%
%NOTE --- IMPORTANT ---
% This algorithm does NOT USE an eigensystem
% similarity transformation. It maps the MODE
% functions to functions of expm, logm and sqrtm,
% which are known to be robust with respect to
% non-diagonalizable ('defective') x
%
%EXA thfm( x ,'cos' ) calculates matrix cosinus
% EVEN IF input matrix x IS NOT DIAGONALIZABLE
%
%ASSOC expm, sqrtm, logm, funm
% Copyright (C) 2001 Rolf Fabian <fab...@tu-cottbus.de>
% 010213
% published under current GNU GENERAL PUBLIC LICENSE
% 2001-03-15 Paul Kienzle
% * extend with inverse functions and power functions
% * optimize handling of real input
% 2011-03-27 Philip Nienhuis
% * Dropped call to funm at bottom (as funm calls thfm...)
function y=thfm(x,M)
#% minimal arg check only
if nargin~=2||~ischar(M)||ischar(x)
usage ("y = thfm (x, MODE)");
endif
## look for known functions of sqrt, log, exp
I = eye(size(x));
match = 1;
len = length(M);
if len==3
if M=='cos',
if (isreal(x)) y = real( expm( i*x ) );
else y = ( expm( i*x ) + expm( -i*x ) ) / 2;
endif
elseif M=='sin',
if (isreal(x)) y = imag( expm( i*x ) );
else y = ( expm( i*x ) - expm( -i*x ) ) / (2*i);
endif
elseif M=='tan',
if (isreal(x)) t = expm( i*x ); y = imag(t)/real(t);
else t = expm( -2*i*x ); y = -i*(I-t)/(I+t);
endif
elseif M=='cot', % == cos/sin
if (isreal(x)) t = expm( i*x ); y = real(t)/imag(t);
else t = expm( -2*i*x ); y = i*(I+t)/(I-t);
endif
elseif M=='sec',
if (isreal(x)) y = inv( real(expm(i*x)) );
else y = inv( expm(i*x)+expm(-i*x) )*2 ;
endif
elseif M=='csc',
if (isreal(x)) y = inv( imag(expm(i*x)) );
else y = inv( expm(i*x)-expm(-i*x) )*2i;
endif
elseif M=='log', y = logm(x);
elseif M=='exp', y = expm(x);
else match = 0;
endif
elseif len==4
if M=='cosh', y = ( expm(x)+expm(-x) )/2;
elseif M=='sinh', y = ( expm(x)-expm(-x) )/2;
elseif M=='tanh' t = expm( -2*x ); y = (I - t)/(I + t);
elseif M=='coth', t = expm( -2*x ); y = (I + t)/(I - t);
elseif M=='sech', y = 2 * inv( expm(x) + expm(-x) );
elseif M=='csch', y = 2 * inv( expm(x) - expm(-x) );
elseif M=='asin', y = -i * logm( i*x + sqrtm(I - x*x) );
elseif M=='acos', y = i * logm( x - i*sqrtm(I - x*x) );
elseif M=='atan', y = -i/2 * logm( (I + i*x)/(I - i*x) );
elseif M=='acot', y = i/2 * logm( (I + i*x)/(i*x - I) );
elseif M=='asec', y = i * logm( ( I - sqrtm(I - x*x) ) / x );
elseif M=='acsc', y = -i * logm( i*( I + sqrtm(I - x*x) ) / x );
elseif M=='sqrt', y = sqrtm(x);
else match = 0;
end
elseif len==5
if M=='asinh', y = logm( x + sqrtm (x*x + I) );
elseif M=='acosh', y = logm( x + sqrtm (x*x - I) );
elseif M=='atanh', y = logm( (I + x)/(I - x) ) / 2;
elseif M=='acoth', y = logm( (I + x)/(x - I) ) / 2;
elseif M=='asech', y = logm( (I + sqrtm (I - x*x)) / x );
elseif M=='acsch', y = logm( (I + sqrtm (I + x*x)) / x );
else match = 0;
endif
else match = 0;
endif
## if no known function found, issue warning
if (match == 0)
warning ("thfm doesn't support function M");
endif
endfunction
## Copyright (C) 2000, 2011 P.R. Nienhuis
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} [ @var{B} ] = funm (@var{A}, @var{F})
## Compute matrix equivalent of function F; F can be a function name or
## a function handle.
##
## For trigonometric and hyperbolic functions, thfm() is automatically
## invoked as that is based on expm() and diagonalization is avoided.
## For other functions diagonalization is invoked, which implies that
## -depending on the properties of input matrix @var{A}- the results
## can be very inaccurate WITHOUT ANY WARNING. For easy diagonizable and
## stable matrices results of funm will be sufficiently accurate.
##
## Note that you should not use funm for 'sqrt', 'log' or 'exp'; instead
## use sqrtm, logm and expm as these are more robust.
##
## Examples:
##
## @example
## B = funm (A, sin);
## (Compute matrix equivalent of sin() )
## @end example
##
## @example
## function bk1 = besselk1 (x)
## bk1 = besselk(x, 1);
## endfunction
## B = funm (A, besselk1);
## (Compute matrix equivalent of bessel function K1(); a helper function
## is needed here to convey extra args for besselk() )
## @end example
##
## @seealso thfm, expm, logm, sqrtm
##
## @end deftypefn
## Author: P.R. Nienhuis <prnienh...@users.sf.net>
## Additions by P. Kienzle, .........
## 2001-03-01 Paul Kienzle
## * Many code improvements
## 2011-03-27 Philip Nienhuis
## * Function handles
## * Texinfo header
## * Fallback to thfm for trig & hyperb funcs to avoid diagonalization
function B = funm (A, name)
persistent thfuncs = {"cos", "sin", "tan", "sec", "csc", "cot", \
"cosh", "sinh", "tanh", "sech", "csch", "coth", \
"acos", "asin", "atan", "asec", "acsc", "acot", \
"acosh", "asinh", "atanh", "asech", "acsch", "acoth", \
}
# Function handle supplied?
try
ishndl = isstruct (functions (name));
fname = functions (name).function;
catch
ishdnl = 0;
fname = ' '
end_try_catch
if (nargin < 2 || (!(ischar (name) || ishndl)) || ischar (A))
usage ("B = funm (A, 'f'");
endif
if (~isempty (find (ismember (thfuncs, fname))))
# Use more robust thfm ()
if (ishndl), name = fname; endif
B = thfm (A, name);
else
# Simply invoke eigenvalues. Note: risk for repeated eigenvalues!!
# Modeled after suggestion by N. Higham (based on R. Davis, 2007)
## FIXME Do we need automatic setting of TOL?
tol = 1.e-15;
[V, D] = eig (A + tol * randn (size(A)));
D = diag (feval (name, diag(D)));
B = V * D / V;
endif
endfunction
------------------------------------------------------------------------------
Enable your software for Intel(R) Active Management Technology to meet the
growing manageability and security demands of your customers. Businesses
are taking advantage of Intel(R) vPro (TM) technology - will your software
be a part of the solution? Download the Intel(R) Manageability Checker
today! http://p.sf.net/sfu/intel-dev2devmar
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
