From octave-sources-request at bevo dot che dot wisc dot edu  Wed Nov  1 18:36:30 2000
Subject: funm
From: Paul Kienzle <pkienzle at kienzle dot powernet dot co dot uk>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Wed, 1 Nov 2000 10:19:45 +0000

Here's an implementation of matlab's funm function.

Paul Kienzle
pkienzle at kienzle dot powernet dot co dot uk

## funm:  Matrix equivalent of function 'name'
##
## Usage:    B = funm(A, name)
##  where    A = square non-singular matrix, provisionally
##               real-valued
##           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
##               element-wise to that vector.
##
## Example:  To compute sqrtm(A), you could use funm(A, 'sqrt')

## 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 eigen-analysis of the network to decouple
## the systems, apply the scalar functions to the eigenvalues,
## and then recombine the systems into a network.

## First version Oct. 2000 by P.R. Nienhuis, 106130 dot 1515 at compuserve dot com
## Additions by P. Kienzle, .........
##
## This file is not part of Octave.
##
## Octave 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.
##
## Octave 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,
## write to the Free Software Foundation, 59 Temple Place -
## Suite 330, Boston, MA 02111-1307, USA.

function B = funm(A, name, ...)

  if (nargin < 2)
    usage ("B = funm (A, 'f' [, args])");
  endif

  dfi = do_fortran_indexing;
  unwind_protect
    do_fortran_indexing = 1;
    n = size (A, 1);
    idx = 1 : n+1 : n*n;
    [V, D] = eig (A);
    D (idx) = feval (name, D (idx), all_va_args);
    B = V * D * inv (V);
  unwind_protect_cleanup
    do_fortran_indexing = dfi;
  end_unwind_protect

endfunction



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------