1 ## Copyright (C) 2000, 2011 P.R. Nienhuis <prnienhuis@users.sf.net>
2 ## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
4 ## This program is free software; you can redistribute it and/or modify it under
5 ## the terms of the GNU General Public License as published by the Free Software
6 ## Foundation; either version 3 of the License, or (at your option) any later
9 ## This program is distributed in the hope that it will be useful, but WITHOUT
10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 ## You should have received a copy of the GNU General Public License along with
15 ## this program; if not, see <http://www.gnu.org/licenses/>.
18 ## @deftypefn {Function File} {@var{B} =} funm (@var{A}, @var{F})
19 ## Compute matrix equivalent of function F; F can be a function name or
22 ## For trigonometric and hyperbolic functions, @code{thfm} is automatically
23 ## invoked as that is based on @code{expm} and diagonalization is avoided.
24 ## For other functions diagonalization is invoked, which implies that
25 ## -depending on the properties of input matrix @var{A}- the results
26 ## can be very inaccurate @emph{without any warning}. For easy diagonizable and
27 ## stable matrices results of funm will be sufficiently accurate.
29 ## Note that you should not use funm for 'sqrt', 'log' or 'exp'; instead
30 ## use sqrtm, logm and expm as these are more robust.
36 ## (Compute matrix equivalent of sin() )
40 ## function bk1 = besselk1 (x)
41 ## bk1 = besselk(x, 1);
43 ## B = funm (A, besselk1);
44 ## (Compute matrix equivalent of bessel function K1(); a helper function
45 ## is needed here to convey extra args for besselk() )
48 ## @seealso{thfm, expm, logm, sqrtm}
51 function B = funm (A, name)
53 persistent thfuncs = {"cos", "sin", "tan", "sec", "csc", "cot", ...
54 "cosh", "sinh", "tanh", "sech", "csch", "coth", ...
55 "acos", "asin", "atan", "asec", "acsc", "acot", ...
56 "acosh", "asinh", "atanh", "asech", "acsch", "acoth", ...
59 ## Function handle supplied?
61 ishndl = isstruct (functions (name));
62 fname = functions (name).function;
68 if (nargin < 2 || (!(ischar (name) || ishndl)) || ischar (A))
69 usage ("B = funm (A, 'f' where A = square matrix and f = function name");
72 if (~isempty (find (ismember (thfuncs, fname))))
73 ## Use more robust thfm ()
74 if (ishndl); name = fname; endif
77 ## Simply invoke eigenvalues. Note: risk for repeated eigenvalues!!
78 ## Modeled after suggestion by N. Higham (based on R. Davis, 2007)
79 ## FIXME Do we need automatic setting of TOL?
81 [V, D] = eig (A + tol * randn (size(A)));
82 D = diag (feval (name, diag(D)));