src/NFfunction/muParser/muParser.cpp
//---------------------------------------------------------------------------
/** \brief Numerically differentiate with regard to a variable. 
      \param [in] a_Var Pointer to the differentiation variable.
      \param [in] a_fPos Position at which the differentiation should take place.
      \param [in] a_fEpsilon Epsilon used for the numerical differentiation.

    Numerical differentiation uses a 5 point operator yielding a 4th order 
    formula. The default value for epsilon is 0.00074 which is
    numerical_limits<double>::epsilon() ^ (1/5) as suggested in the muparser
    forum:

    http://sourceforge.net/forum/forum.php?thread_id=1994611&forum_id=462843
  */
\brief Numerically differentiate with regard to a variable.
\param [in] a_Var Pointer to the differentiation variable.
\param [in] a_fPos Position at which the differentiation should take place.
\param [in] a_fEpsilon Epsilon used for the numerical differentiation.

Numerical differentiation uses a 5 point operator yielding a 4th order
formula. The default value for epsilon is 0.00074 which is
numerical_limits::epsilon() ^ (1/5) as suggested in the muparser
forum.


\ brief Numerically differentiate with regard to a variable . \ param [ in ] a_Var Pointer to the differentiation variable . \ param [ in ] a_fPos Position at which the differentiation should take place . \ param [ in ] a_fEpsilon Epsilon used for the numerical differentiation . Numerical differentiation uses a 5 point operator yielding a 4th order formula . The default value for epsilon is 0 . 00074 which is numerical_limits :: epsilon () ^ ( 1 / 5 ) as suggested in the muparser forum .
value_type Parser::Diff(value_type *a_Var, 
                          value_type  a_fPos, 
                          value_type  a_fEpsilon) const
  {
    value_type fRes(0), 
               fBuf(*a_Var), 
               f[4] = {0,0,0,0};
    *a_Var = a_fPos+2 * a_fEpsilon;  f[0] = Eval();
    *a_Var = a_fPos+1 * a_fEpsilon;  f[1] = Eval();
    *a_Var = a_fPos-1 * a_fEpsilon;  f[2] = Eval();
    *a_Var = a_fPos-2 * a_fEpsilon;  f[3] = Eval();
    *a_Var = fBuf; 
    fRes = (-f[0] + 8*f[1] - 8*f[2] + f[3]) / (12*a_fEpsilon);
    return fRes;
  }
value_type Parser :: Diff ( value_type * a_Var , value_type a_fPos , value_type a_fEpsilon ) const { value_type fRes ( 0 ) , fBuf ( * a_Var ) , f [ 4 ] = { 0 , 0 , 0 , 0 } ; * a_Var = a_fPos + 2 * a_fEpsilon ; f [ 0 ] = Eval ( ) ; * a_Var = a_fPos + 1 * a_fEpsilon ; f [ 1 ] = Eval ( ) ; * a_Var = a_fPos - 1 * a_fEpsilon ; f [ 2 ] = Eval ( ) ; * a_Var = a_fPos - 2 * a_fEpsilon ; f [ 3 ] = Eval ( ) ; * a_Var = fBuf ; fRes = ( - f [ 0 ] + 8 * f [ 1 ] - 8 * f [ 2 ] + f [ 3 ] ) / ( 12 * a_fEpsilon ) ; return fRes ; }
\brief Numerically differentiate with regard to a variable.
\ brief Numerically differentiate with regard to a variable .
// restore variable
