Re: numerical differentiation?

[Gerard Jungman <jungman at lanl dot gov>]

I'm not sure which methods you are talking about. There are
different notions, depending on usage. Do you mean one
or all of the following?

  /* differentiate a function at a point */
  gsl_diff_function(const gsl_function * f, double x0, double *
dydx_x0);

  /* differentiate equally spaced data in a table */
  gsl_diff_table_unif(
    const double * y_data,
    size_t length,
    double delta_x,
    double * dydx_data
    );

  /* differentiate a table */
  gsl_diff_table(
    const double * y_data,
    const double * x_data,
    size_t length,
    double * dydx_data
    );

  /* evaluate a derivative poinwise, using given data */
  gsl_diff_point(
    const double * y_data,
    const double * x_data,
    size_t length,
    double x0,
    double * dydx_x0
    );


Making a general purpose interface for this sort of functionality
will require a little thought. The following issues arise:

  o Both constant and non-constant spacing should be supported
    for the table methods.

  o 'gsl_diff_table' and 'gsl_diff_table_unif' are not complete
interfaces.
    They leave open the question of what happens at the endpoints,
    for instance whether or not the data is assumed periodic (which will
    in itself be confusing when the spacing is not constant).

  o 'gsl_diff_function' might have trouble with the domain of the
function.
    You don't want the needed evaluations of the function to wander
    outside the domain [ex. f(x) = log(x), x in [0, infinity), should
    never try to evaluate f(-1.0e-56) or something]. Brian: is there any
    facility hovering around gsl_function for declaring a domain?

  o People who use 'gsl_diff_point' repeatedly on the same data,
    but with different x0, may want an optimization based on
    some kind of cached state, depending on the complexity
    of the evaluation. We have to decide whether or not
    that is worth supporting. There is already something like that
    in the chebyshev evaluation functions, where you get much
    of it for free because you have to cache the chebyshev fit
    anyway. Of course, in the case of chebshev, you only get the
    one evaluation method, which is only obliguely related to
    difference quotients. My feeling is that it is not worth
    supporting any sort of cached state for differencing, since
    the implementation complexity will probably be low. But that
    really depends on the set of methods to be encapsulated;
    if somebody comes along with a whiz-bang but expensive method,
    then we might feel bad about not having an interface flexible
    enough to support it properly.

  o As Brian points out, you really want to have some error estimate.
    I'm not sure what the theory is for this. For 'gsl_diff_function',
    presumably you do some kind of extrapolation (Richardson?), and
produce an error
    estimate from that. But I don't know what you do for fixed data;
that
    requires some assumption about the class of functions. Also,
    error estimates will complicate the interface for something
    like 'gsl_diff_table'; but maybe it could be something simple like
      gsl_diff_table(
        const double * y_data,
        const double * x_data,
        size_t length,
        double * dydx_data,
        double * dydx_error
        );
    where dydx_error==0 indicates no error estimate will be provided.


-- 
G. Jungman


Brian Gough wrote:
> 
> Yes, that would be very useful, especially if there is an estimate of
> the error in the derivative.
> 
> Don't feel constrained by omissions in the TODO file -- it is not an
> exhaustive list.  As a general rule anything that's covered by
> existing proprietary libraries is desirable for GSL.
> 
> Dave Morrison writes:
>  > Hi all,
>  >
>  > Is there any interest in adding routines for numerical differentiation
>  > to GSL?