centered_difference.h

Go to the documentation of this file.

// -*- C++ -*-

#if !defined(__numerical_centered_difference_h__)
#define __numerical_centered_difference_h__

#include "../defs.h"

#include "../../ads/array/FixedArray.h"

#include <functional>
#include <limits>

#include <cmath>

BEGIN_NAMESPACE_NUMERICAL

//-----------------------------------------------------------------------------
// @{

template <class Functor>
inline
void
derivative_centered_difference
( const Functor& f, 
  const typename Functor::argument_type x, 
  typename Functor::result_type& deriv, 
  const typename Functor::argument_type 
  delta = std::pow( std::numeric_limits<typename Functor::argument_type>::
                    epsilon(), 1.0 / 3.0 ) ) 
{
  const typename Functor::argument_type p = x + delta;
  const typename Functor::argument_type n = x - delta;
  deriv = ( f( p ) - f( n ) ) / ( p - n );
}

template <class Functor>
inline
typename Functor::result_type
derivative_centered_difference
( const Functor& f, 
  const typename Functor::argument_type x, 
  const typename Functor::argument_type 
  delta = std::pow( std::numeric_limits<typename Functor::argument_type>::
                    epsilon(), 1.0 / 3.0 ) ) 
{
  typename Functor::result_type deriv;
  derivative_centered_difference( f, x, deriv, delta );
  return deriv;
}

template <class Functor>
class DerivativeCenteredDifference
  : public std::unary_function<typename Functor::argument_type,
                               typename Functor::result_type>
{
  //
  // Types
  //

private:

  typedef std::unary_function<typename Functor::argument_type,
                              typename Functor::result_type> base_type;

public:

  typedef Functor function_type;
  typedef typename base_type::argument_type argument_type;
  typedef typename base_type::result_type result_type;

private:
  
  //
  // Member data.
  //

  const Functor& _f;
  
  argument_type _delta;

  //
  // Not implemented.
  //

  // Default constructor not implemented.
  DerivativeCenteredDifference();

  // Assignment operator not implemented.
  DerivativeCenteredDifference&
  operator=( const DerivativeCenteredDifference& x );

public:

  DerivativeCenteredDifference
  ( const function_type& f, 
    const argument_type delta = 
    std::pow( std::numeric_limits<argument_type>::epsilon(), 1.0 / 3.0 ) ) : 
    _f( f ),
    _delta( delta )
  {}

  DerivativeCenteredDifference( const DerivativeCenteredDifference& x ) : 
    _f( x._f ),
    _delta( x._delta )
  {}
  
  void
  operator()( const argument_type x, result_type& deriv ) const 
  {
    derivative_centered_difference( _f, x, deriv, _delta );
  }

  result_type
  operator()( const argument_type x ) const 
  {
    return derivative_centered_difference( _f, x, _delta );
  }

  argument_type
  delta() const
  {
    return _delta;
  }

  void
  set_delta( const argument_type delta )
  {
    _delta = delta;
  }
};

// @}


//-----------------------------------------------------------------------------
// @{

template <int N, class Functor, typename T>
inline
void
gradient_centered_difference
( const Functor& f, 
  typename Functor::argument_type x,
  ads::FixedArray<N,typename Functor::result_type>& gradient, 
  const T delta = std::pow( std::numeric_limits<T>::epsilon(), 1.0 / 3.0 ) ) 
{
  T temp, xp, xn;
  typename Functor::result_type fp, fn;
  for ( int n = 0; n != N; ++n ) {
    temp = x[n];
    // Positive offset.
    xp = x[n] += delta;
    fp = f( x );
    // Negative offset.
    x[n] = temp;
    xn = x[n] -= delta;
    fn = f( x );
    // Reset
    x[n] = temp;
    // Centered difference.
    gradient[n] = ( fp - fn ) / ( xp - xn );
  }
}


template <int N, class Functor, typename T>
inline
ads::FixedArray<N,typename Functor::result_type>
gradient_centered_difference
( const Functor& f, 
  typename Functor::argument_type x,
  const T delta = std::pow( std::numeric_limits<T>::epsilon(), 1.0 / 3.0 ) ) 
{
  ads::FixedArray<N,typename Functor::result_type> gradient;
  gradient_centered_difference( f, x, gradient, delta );
  return gradient;
}

template <int N, 
          class Functor, 
          typename T = typename Functor::argument_type::value_type>
class GradientCenteredDifference : 
  public std::unary_function
< typename Functor::argument_type,
  ads::FixedArray<N,typename Functor::result_type> >
{
private:
  typedef std::unary_function
  < typename Functor::argument_type,
    ads::FixedArray<N,typename Functor::result_type> >
  base_type;
  typedef T number_type;

private:
  
  //
  // Member data.
  //

  const Functor& _f;
  
  number_type _delta;

  //
  // Not implemented.
  //

  // Default constructor not implemented.
  GradientCenteredDifference();

  // Assignment operator not implemented.
  GradientCenteredDifference&
  operator=( const GradientCenteredDifference& x );

public:

  typedef Functor function_type;
  typedef typename base_type::argument_type argument_type;
  typedef typename base_type::result_type result_type;

  GradientCenteredDifference
  ( const function_type& f, 
    const number_type delta = 
    std::pow( std::numeric_limits<T>::epsilon(), 1.0 / 3.0 ) ) : 
    _f( f ),
    _delta( delta )
  {}

  GradientCenteredDifference( const GradientCenteredDifference& x ) : 
    _f( x._f ),
    _delta( x._delta )
  {}
  
  void
  operator()( const argument_type x, result_type& deriv ) const 
  {
    gradient_centered_difference( _f, x, deriv, _delta );
  }

  result_type
  operator()( const argument_type x ) const 
  {
    return gradient_centered_difference<N>( _f, x, _delta );
  }

  number_type
  delta() const
  {
    return _delta;
  }

  void
  set_delta( const number_type delta )
  {
    _delta = delta;
  }
};

// @}

END_NAMESPACE_NUMERICAL

#endif

---