octave_packages/linear-algebra-2.2.0/thfm.m

   1 ## Copyright (C) 2001 Rolf Fabian <fabian@tu-cottbus.de>
   2 ## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
   3 ## Copyright (C) 2011 Philip Nienhuis <pr.nienhuis@hccnet.nl>
   4 ## Copyright (C) 2011 Carnë Draug <carandraug+dev@gmail.com>
   5 ##
   6 ## This program is free software; you can redistribute it and/or modify it under
   7 ## the terms of the GNU General Public License as published by the Free Software
   8 ## Foundation; either version 3 of the License, or (at your option) any later
   9 ## version.
  10 ##
  11 ## This program is distributed in the hope that it will be useful, but WITHOUT
  12 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  14 ## details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License along with
  17 ## this program; if not, see <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn{Function File} {@var{y} =} thfm (@var{x}, @var{mode})
  21 ## Trigonometric/hyperbolic functions of square matrix @var{x}.
  22 ##
  23 ## @var{mode} must be the name of a function. Valid functions are 'sin', 'cos',
  24 ## 'tan', 'sec', 'csc', 'cot' and all their inverses and/or hyperbolic variants,
  25 ## and 'sqrt', 'log' and 'exp'.
  26 ##
  27 ## The code @code{thfm (x, 'cos')} calculates matrix cosinus @emph{even if} input
  28 ## matrix @var{x} is @emph{not} diagonalizable.
  29 ##
  30 ## @emph{Important note}:
  31 ## This algorithm does @emph{not} use an eigensystem similarity transformation. It
  32 ## maps the @var{mode} functions to functions of @code{expm}, @code{logm} and
  33 ## @code{sqrtm}, which are known to be robust with respect to non-diagonalizable
  34 ## ('defective') @var{x}.
  35 ##
  36 ## @seealso{funm}
  37 ## @end deftypefn
  38
  39 function y = thfm (x,M)
  40   ## minimal arg check only
  41   if ( nargin != 2 || !ischar (M) || ischar (x) )
  42     print_usage;
  43   endif
  44
  45   ## look for known functions of sqrt, log, exp
  46   I = eye (size (x));
  47
  48   switch (M)
  49     case {'cos'}
  50       if (isreal(x))  y = real( expm( i*x ) );
  51       else            y = ( expm( i*x ) + expm( -i*x ) ) / 2;
  52       endif
  53     case {'sin'}
  54       if (isreal(x))  y = imag( expm( i*x ) );
  55       else            y = ( expm( i*x ) - expm( -i*x ) ) / (2*i);
  56       endif
  57     case {'tan'}
  58       if (isreal(x))  t = expm( i*x );    y = imag(t)/real(t);
  59       else            t = expm( -2*i*x ); y = -i*(I-t)/(I+t);
  60       endif
  61     case {'cot'}
  62       if (isreal(x))  t = expm( i*x );    y = real(t)/imag(t);
  63       else            t = expm( -2*i*x ); y = i*(I+t)/(I-t);
  64       endif
  65     case {'sec'}
  66       if (isreal(x))  y = inv( real(expm(i*x)) );
  67       else            y = inv( expm(i*x)+expm(-i*x) )*2 ;
  68       endif
  69     case {'csc'}
  70       if (isreal(x))  y = inv( imag(expm(i*x)) );
  71       else            y = inv( expm(i*x)-expm(-i*x) )*2i;
  72       endif
  73     case {'log'}      y = logm(x);
  74     case {'exp'}      y = expm(x);
  75     case {'cosh'}     y = ( expm(x)+expm(-x) )/2;
  76     case {'sinh'}     y = ( expm(x)-expm(-x) )/2;
  77     case {'tanh'}     t = expm( -2*x ); y = (I - t)/(I + t);
  78     case {'coth'}     t = expm( -2*x ); y = (I + t)/(I - t);
  79     case {'sech'}     y = 2 * inv( expm(x) + expm(-x) );
  80     case {'csch'}     y = 2 * inv( expm(x) - expm(-x) );
  81     case {'asin'}     y = -i * logm( i*x + sqrtm(I - x*x) );
  82     case {'acos'}     y =  i * logm( x - i*sqrtm(I - x*x) );
  83     case {'atan'}     y = -i/2 * logm( (I + i*x)/(I - i*x) );
  84     case {'acot'}     y =  i/2 * logm( (I + i*x)/(i*x - I) );
  85     case {'asec'}     y = i * logm( ( I - sqrtm(I - x*x) ) / x );
  86     case {'acsc'}     y = -i * logm( i*( I + sqrtm(I - x*x) ) / x );
  87     case {'sqrt'}     y = sqrtm(x);
  88     case {'asinh'}    y = logm( x + sqrtm (x*x + I) );
  89     case {'acosh'}    y = logm( x + sqrtm (x*x - I) );
  90     case {'atanh'}    y = logm( (I + x)/(I - x) ) / 2;
  91     case {'acoth'}    y = logm( (I + x)/(x - I) ) / 2;
  92     case {'asech'}    y = logm( (I + sqrtm (I - x*x)) / x );
  93     case {'acsch'}    y = logm( (I + sqrtm (I + x*x)) / x );
  94     otherwise
  95       error ("thfm doesn't support function %s - try to use funm instead.", M);
  96   endswitch
  97
  98 endfunction