Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.
Close
From: <prnienhuis@us...>  20110729 09:26:10

Revision: 8422 http://octave.svn.sourceforge.net/octave/?rev=8422&view=rev Author: prnienhuis Date: 20110729 09:26:03 +0000 (Fri, 29 Jul 2011) Log Message:  Calls thfm.m for trig & hyperbolic functions w/o invoking diagonalization. Cleaned up texinfo. Can cope with function handles for functions in thfm Modified Paths:  trunk/octaveforge/main/linearalgebra/inst/funm.m Modified: trunk/octaveforge/main/linearalgebra/inst/funm.m ===================================================================  trunk/octaveforge/main/linearalgebra/inst/funm.m 20110729 09:24:04 UTC (rev 8421) +++ trunk/octaveforge/main/linearalgebra/inst/funm.m 20110729 09:26:03 UTC (rev 8422) @@ 1,62 +1,98 @@ ## Copyright (C) 2000 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, 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 this program; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>;.  ## funm: Matrix equivalent of function 'name' ## ## Usage: B = funm(A, name) ## where A = square nonsingular matrix, provisionally ## realvalued ## B = square result matrix ## name = string, name of function to apply to A. ## args = any arguments to pass to function 'name' ## The function must accept a vector and apply ## elementwise to that vector. ## ## Example: To compute sqrtm(A), you could use funm(A, 'sqrt') ## ## Note that you should not use funm for 'sqrt', 'log' or 'exp'; instead ## use sqrtm, logm and expm which are more robust. Similarly, ## trigonometric and hyperbolic functions (cos, sin, tan, cot, sec, csc, ## cosh, sinh, tanh, coth, sech, csch) are better handled by thfm(A, ## name), which defines them in terms of the more robust expm.  ## NOTE: the following comments are withheld until they can be verified ## ## If you have a network of coupled systems, where for the individual ## systems a solution exists in terms of scalar variables, in many ## cases the network might be solved using the same form of the ## solution but with substituting the matrix equivalent of the function ## applied to the scalar variables. ## The approach is to do an eigenanalysis of the network to decouple ## the systems, apply the scalar functions to the eigenvalues, ## and then recombine the systems into a network.  ## Author: P.R. Nienhuis <106130.1515@...> ## Additions by P. Kienzle, ......... ## 20010301 Paul Kienzle ## * generate error for repeated eigenvalues  function B = funm(A, name)   if (nargin != 2  !ischar(name)  ischar(A))  usage ("B = funm (A, 'f' [, args])");  endif   [V, D] = eig (A);  D = diag (feval (name, diag(D)));  B = V * D * inv (V);  endfunction +## Copyright (C) 2000, 2011 P.R. Nienhuis <prnienhuis@...> +## +## 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 <prnienhuis@...> (somewhere in 2000) +## Additions by P. Kienzle, ......... +## 20010301 Paul Kienzle +## * Many code improvements +## 20110327 Philip Nienhuis +## * Function handles +## * Texinfo header +## * Fallback to thfm for trig & hyperb funcs to avoid diagonalization +## 20110729 Philip Nienhuis +## * Layout, cleanup (tabs, alignment, ...) + +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' where A = square matrix and f = function name"); + 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.e15; + [V, D] = eig (A + tol * randn (size(A))); + D = diag (feval (name, diag(D))); + B = V * D / V; + endif + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 