Octave interface to gsl_deriv for finite differences with adaptive step size

This is a multi-part message in MIME format.
--------------000509070608050809000108
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 7bit

G'day all,

Please find attached an Octave interface to the GSL numerical
differentiation function gsl_deriv(). This interface provides an
adaptive step size forward difference approximation to the gradient of a
vector-valued function. Compile it using

mkoctfile -v gsl_deriv.cc -lgsl

For an relevant article on the topic, see
http://www.nezumi.demon.co.uk/consult/deriv.htm

Cheers,
Paul.



--------------000509070608050809000108
Content-Type: text/x-c++src;
 name="gsl_deriv.cc"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="gsl_deriv.cc"

#include <octave/oct.h>
#include <octave/parse.h>
#include <gsl/gsl_deriv.h>

//$Id: gsl_deriv.cc,v 1.3 2006/07/25 13:40:02 pkoufalas Exp $

std::string fstr="";
octave_value_list* fargsinptr;
double* x_=NULL;
double fun (double xj, void * jay) {
  int* j = (int *)jay;
  x_[*j]=xj;
  octave_value_list fargsout;
  fargsout=feval(fstr,*fargsinptr,1);
#ifdef DEBUG
  octave_value((*fargsinptr)(0)).print(std::cout);
  fargsout(0).print(std::cout);
#endif
  return fargsout(0).scalar_value();
}

DEFUN_DLD (gsl_deriv, args, ,
           "[deriv,abserr] = gsl_deriv (fun,x,h0,method,varargin)\n\
Calculate finite differences using central, forward or backward methods\n\
with adaptive step size using GSL routine gsl_deriv.\n\
 fun, name of function, string\n\
 x, point at which to evaluate derivative, column vector\n\
 h0, initial step size, scalar\n\
 method, 0 for forward,\n\
         1 for backward,\n\
         2 for central.\n\
 varargin, additional args to pass to fun.\n") {

    octave_value_list retval; 

    if (args.length() < 4) { 
        print_usage("gsl_deriv");
        return retval;
    }
    if (error_state) {
      print_usage("gsl_deriv: invalid parameters specified.");
      return retval;
    }

    octave_value_list fargsin;
    fargsinptr = &fargsin;
    const Matrix x ( args(1).matrix_value() );
    fargsin.append(x);

    x_ = (double *)x.fortran_vec();

    int i;
    int nargs = args.length()-4;
    for (i=0;i<nargs;i++)
      fargsin.append(args(i+4));

    gsl_function F;
    fstr = args(0).string_value();

    F.function = &fun;

    double h = args(2).scalar_value();
    int method = args(3).int_value();
    int ndim = x.rows();
    double result[ndim]; 
    double abserr[ndim];
    double x__[ndim];
    for (i=0;i<ndim;i++) // save x!
      x__[i]=x_[i]; 

    switch (method) {
    case 0:
      for (i=0;i<ndim;i++) {
	F.params = &i;
	gsl_deriv_forward(&F,x(i,0),h,&(result[i]),&(abserr[i]));
      }
      break;
    case 1:
    case 2:
    default:
      error("sorry, only forward differences in this version.");
    }

    Matrix d (ndim,1);
    Matrix e (ndim,1);
    Matrix f (ndim,1);
    for (i=0;i<ndim;i++) {
      d(i,0)=result[i];
      e(i,0)=abserr[i];
      f(i,0)=x_[i]-x__[i];
      x_[i]=x__[i]; // restore x! i.e. undo changes to x through pointer
      // trickery in using const Matrix x, x_ = (double *)x.fortran_vec()...
    }

    retval.append(d);
    retval.append(e);
    retval.append(f);
    return retval;

}

--------------000509070608050809000108
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

_______________________________________________
Octave-sources mailing list
[email protected]
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

--------------000509070608050809000108--