here is a parametric version of fsolve

[Craig Earls]

attached is the source code for a drop in replacement for fsolve. It add
the following capability:

out=fsolve("foo",x,a,b,c,...)

a,b,c etc are parameters that will be passed directly to foo, while
fsolve manipulates the vector x to solve the equation.

To use this you must either recompile Octave with this in place of
src/fsolve.cc
 
or

mkoctfile fsolve.cc

if you have dynamic loading enabled already. Make certain to save your
old version, and you will need to symbolically link fsolve_options.oct
to it in order to set the tolerance

If you would like to try this without replacing your fsolve, simply do a
search and replace for fsolve and reaplce it with some other name such
as fsolve2, it will then have the same syntax with the exception of the
name change. (change the file name to the same thing)

-------begin fsolve.cc
/*

Copyright (C) 1996 John W. Eaton

Modified 3 Mar 97 by Craig P. Earls, all rights released to John W.
Eaton

This file is 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.

*/

#include <octave/config.h>

#include <string>

#include <iostream.h>

#include <octave/NLEqn.h>

#include <octave/defun-dld.h>
#include <octave/error.h>
#include <octave/gripes.h>
#include <octave/help.h>
#include <octave/pager.h>
#include <octave/pt-fvc.h>
#include <octave/oct-obj.h>
#include <octave/utils.h>
#include <octave/variables.h>

// Global pointer for user defined function required by hybrd1.
static tree_fvc *fsolve_fcn;

static NLEqn_options fsolve_opts;
static octave_value_list extras;
static int numextras;

int
hybrd_info_to_fsolve_info (int info)
{
  switch (info)
    {
    case -1:
      info = -2;
      break;

    case 0:
      info = -1;
      break;

    case 1:
      break;

    case 2:
      info = 4;
      break;

    case 3:
    case 4:
    case 5:
      info = 3;
      break;

    default:
      panic_impossible ();
      break;
    }

  return info;
}

ColumnVector
fsolve_user_function (const ColumnVector& x)
{
  ColumnVector retval;

  int n = x.capacity ();

  octave_value_list args;
  
//modified 3 Mar 97 CPE
  
  args.resize (1+numextras);

  if (n > 1)
    {
      Matrix m (n, 1);
      for (int i = 0; i < n; i++)
        m (i, 0) = x (i);
      octave_value vars (m);
      args(0) = vars;
    }
  else
    {
      double d = x (0);
      octave_value vars (d);
      args(0) = vars;
    }
  
      //modified 3mar 97 CPE
  if(numextras){
      for(int i=0;i<numextras;i++){
          args(i+1)=extras(i);
      }
  }
  
      
  if (fsolve_fcn)
    {
      octave_value_list tmp = fsolve_fcn->eval (0, 1, args);
      if (tmp.length () > 0 && tmp(0).is_defined ())
        {
          retval = tmp(0).vector_value ();

          if (error_state || retval.length () <= 0)
            gripe_user_supplied_eval ("fsolve");
        }
      else
        gripe_user_supplied_eval ("fsolve");
    }

  return retval;
}

DEFUN_DLD (fsolve, args, nargout,
  "Solve nonlinear equations using Minpack.  Usage:\n\
\n\
  [X, INFO] = fsolve (F, X0, X1, X2,...)\n\
\n\
Where the first argument is the name of the  function to call to\n\
compute the vector of function values.  It must have the form\n\
\n\
  y = f (x0,x1,x2,x3,...)
\n\
where y and x0 are vectors. Additional arguments are passed to \n\
the function as parameters")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin < 2 || nargout > 3)
    {
      print_usage ("fsolve");
      return retval;
    }

  fsolve_fcn = is_valid_function (args(0), "fsolve", 1);
  if (! fsolve_fcn)
    return retval;

  ColumnVector x = args(1).vector_value ();

      //add code to copy extra parameters to pass through
  if (nargin>2){
      extras.resize(nargin-2);// extras will pass the extra arguments
through
      numextras=nargin-2;
      
      for(int i=2;i<nargin;i++){
          extras(i-2)=args(i).vector_value();
      }
  }
  else{
      extras.resize(0);
      numextras=0;
  }
  
      
  if (error_state)
    {
      error ("fsolve: expecting vector as second argument");
      return retval;
    }

  if (nargout > 2)
    warning ("fsolve: can't compute path output yet");

  NLFunc foo_fcn (fsolve_user_function);
  NLEqn foo (x, foo_fcn);
  foo.set_options (fsolve_opts);

  int info;
  ColumnVector soln = foo.solve (info);

  info = hybrd_info_to_fsolve_info (info);

  retval.resize (nargout ? nargout : 1);
  retval(0) = soln, 1;

  if (nargout > 1)
    retval(1) = (double) info;

  return retval;
}

typedef void (NLEqn_options::*d_set_opt_mf) (double);
typedef double (NLEqn_options::*d_get_opt_mf) (void);

#define MAX_TOKENS 1

struct NLEQN_OPTIONS
{
  const char *keyword;
  const char *kw_tok[MAX_TOKENS + 1];
  int min_len[MAX_TOKENS + 1];
  int min_toks_to_match;
  d_set_opt_mf d_set_fcn;
  d_get_opt_mf d_get_fcn;
};

static NLEQN_OPTIONS fsolve_option_table [] =
{
  { "tolerance",
    { "tolerance", 0, },
    { 1, 0, }, 1,
    NLEqn_options::set_tolerance,
    NLEqn_options::tolerance, },

  { 0,
    { 0, 0, },
    { 0, 0, }, 0,
    0, 0, },
};

static void
print_fsolve_option_list (ostream& os)
{
  print_usage ("fsolve_options", 1);

  os << "\n"
     << "Options for fsolve include:\n\n"
     << "  keyword                                  value\n"
     << "  -------                                  -----\n\n";

  NLEQN_OPTIONS *list = fsolve_option_table;

  const char *keyword;
  while ((keyword = list->keyword) != 0)
    {
      os.form ("  %-40s ", keyword);

      double val = (fsolve_opts.*list->d_get_fcn) ();
      if (val < 0.0)
        os << "computed automatically";
      else
        os << val;

      os << "\n";
      list++;
    }

  os << "\n";
}

static void
set_fsolve_option (const string& keyword, double val)
{
  NLEQN_OPTIONS *list = fsolve_option_table;

  while (list->keyword != 0)
    {
      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
                                list->min_toks_to_match, MAX_TOKENS))
        {
          (fsolve_opts.*list->d_set_fcn) (val);

          return;
        }
      list++;
    }

  warning ("fsolve_options: no match for `%s'", keyword.c_str ());
}

static octave_value_list
show_fsolve_option (const string& keyword)
{
  octave_value_list retval;

  NLEQN_OPTIONS *list = fsolve_option_table;

  while (list->keyword != 0)
    {
      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
                                list->min_toks_to_match, MAX_TOKENS))
        {
          return (fsolve_opts.*list->d_get_fcn) ();
        }
      list++;
    }

  warning ("fsolve_options: no match for `%s'", keyword.c_str ());

  return retval;
}

DEFUN_DLD (fsolve_options, args, ,
  "fsolve_options (KEYWORD, VALUE)\n\
\n\
Set or show options for fsolve.  Keywords may be abbreviated\n\
to the shortest match.")
{
  octave_value_list retval;

  int nargin = args.length ();

  if (nargin == 0)
    {
      print_fsolve_option_list (octave_stdout);
      return retval;
    }
  else if (nargin == 1 || nargin == 2)
    {
      string keyword = args(0).string_value ();

      if (! error_state)
        {
          if (nargin == 1)
            return show_fsolve_option (keyword);
          else
            {
              double val = args(1).double_value ();

              if (! error_state)
                {
                  set_fsolve_option (keyword, val);
                  return retval;
                }
            }
        }
    }

  print_usage ("fsolve_options");

  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
-- 
-----------------------------------------------------------------
Craig P Earls                                 address@hidden
LT US Navy, MIT Ocean Engineering             address@hidden
-----------------------------------------------------------------