octave_packages/linear-algebra-2.2.0/funm.m

   1 ## Copyright (C) 2000, 2011 P.R. Nienhuis <prnienhuis@users.sf.net>
   2 ## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
   3 ##
   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
   7 ## version.
   8 ##
   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
  12 ## details.
  13 ##
  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/>.
  16
  17 ## -*- texinfo -*-
  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
  20 ## a function handle.
  21 ##
  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.
  28 ##
  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.
  31 ##
  32 ## Examples:
  33 ##
  34 ## @example
  35 ##   B = funm (A, sin);
  36 ##   (Compute matrix equivalent of sin() )
  37 ## @end example
  38 ##
  39 ## @example
  40 ##   function bk1 = besselk1 (x)
  41 ##      bk1 = besselk(x, 1);
  42 ##   endfunction
  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() )
  46 ## @end example
  47 ##
  48 ## @seealso{thfm, expm, logm, sqrtm}
  49 ## @end deftypefn
  50
  51 function B = funm (A, name)
  52
  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", ...
  57                         }
  58
  59   ## Function handle supplied?
  60   try
  61     ishndl = isstruct (functions (name));
  62     fname = functions (name).function;
  63   catch
  64     ishdnl = 0;
  65     fname = ' '
  66   end_try_catch
  67
  68   if (nargin < 2 || (!(ischar (name) || ishndl)) || ischar (A))
  69     usage ("B = funm (A, 'f'  where A = square matrix and f = function name");
  70   endif
  71
  72   if (~isempty (find (ismember (thfuncs, fname))))
  73     ## Use more robust thfm ()
  74     if (ishndl); name = fname; endif
  75     B = thfm (A, name);
  76   else
  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?
  80     tol = 1.e-15;
  81     [V, D] = eig (A + tol * randn (size(A)));
  82     D = diag (feval (name, diag(D)));
  83     B = V * D / V;
  84   endif
  85
  86 endfunction