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>
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
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
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/>.
20 ## @deftypefn{Function File} {@var{y} =} thfm (@var{x}, @var{mode})
21 ## Trigonometric/hyperbolic functions of square matrix @var{x}.
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'.
27 ## The code @code{thfm (x, 'cos')} calculates matrix cosinus @emph{even if} input
28 ## matrix @var{x} is @emph{not} diagonalizable.
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}.
39 function y = thfm (x,M)
40 ## minimal arg check only
41 if ( nargin != 2 || !ischar (M) || ischar (x) )
45 ## look for known functions of sqrt, log, exp
50 if (isreal(x)) y = real( expm( i*x ) );
51 else y = ( expm( i*x ) + expm( -i*x ) ) / 2;
54 if (isreal(x)) y = imag( expm( i*x ) );
55 else y = ( expm( i*x ) - expm( -i*x ) ) / (2*i);
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);
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);
66 if (isreal(x)) y = inv( real(expm(i*x)) );
67 else y = inv( expm(i*x)+expm(-i*x) )*2 ;
70 if (isreal(x)) y = inv( imag(expm(i*x)) );
71 else y = inv( expm(i*x)-expm(-i*x) )*2i;
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 );
95 error ("thfm doesn't support function %s - try to use funm instead.", M);