derivative.hpp

Go to the documentation of this file.

// ====================================================================
// This file is part of FlexibleSUSY.
//
// FlexibleSUSY 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 3 of the License,
// or (at your option) any later version.
//
// FlexibleSUSY 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 FlexibleSUSY.  If not, see
// <http://www.gnu.org/licenses/>.
// ====================================================================
 
#ifndef DERIVATIVE_H
#define DERIVATIVE_H
 
#include <cmath>
#include <limits>
#include <vector>
 
namespace flexiblesusy {
 
template <class F, class A>
auto derivative_forward_fx(const F& f, A x, decltype(f(x)) fx, A eps = std::numeric_limits<A>::epsilon())
   -> decltype(f(x))
{
   const A h = std::fabs(x) < eps ? eps : std::sqrt(eps) * x;
   volatile const A xph = x + h; // avoid away optimization
   const A dx = xph - x;
   return (f(xph) - fx) / dx;
}
 
template <class F, class A>
auto derivative_backward_fx(const F& f, A x, decltype(f(x)) fx, A eps = std::numeric_limits<A>::epsilon())
   -> decltype(f(x))
{
   const A h = std::fabs(x) < eps ? eps : std::sqrt(eps) * x;
   volatile const A xph = x - h; // avoid away optimization
   const A dx = x - xph;
   return (fx - f(xph)) / dx;
}
 
template <int Order, int sign, class F, class A>
auto derivative_one_sided(const F& f, A x, A eps = std::numeric_limits<A>::epsilon()) -> decltype(f(x))
{
   static_assert(Order <= 7, "1st forward derivative with order > 7 not implemented");
   static_assert(sign == -1 || sign == +1, "sign must be either +1 or -1 for one-sided derivative");
 
   using return_type = decltype(f(x));
 
   // coefficients from Math. Comp. 51 (1988), 699-706, Table 3
   // DOI: http://dx.doi.org/10.1090/S0025-5718-1988-0935077-0
   static const std::vector<std::vector<double> > coeffs = {
      {-1., 1.},
      {-3./2., 2., -1./2.},
      {-11./6., 3., -3./2., 1./3.},
      {-25./12., 4., -3., 4./3., -1./4.},
      {-137./60., 5., -5., 10./3., -5./4., 1./5.},
      {-49./20., 6., -15./2., 20./3., -15./4., 6./5., -1./6.},
      {-363./140., 7., -21./2., 35./3., -35./4., 21./5., -7./6., 1./7.},
      {-761./280., 8., -14., 56./3., -35./2., 56./5., -14./3., 8./7., -1./8.}
   };
 
   const A h = std::fabs(x) < eps ? eps : std::sqrt(eps) * x;
   return_type result = 0;
 
   for (int i = 0; i < Order + 2; i++) {
      const double coeff = coeffs[Order][i];
      result += sign * coeff * f(x + sign*i*h);
   }
 
   return result / h;
}
 
template <int Order, class F, class A>
auto derivative_forward(const F& f, A x, A eps = std::numeric_limits<A>::epsilon()) -> decltype(f(x))
{
   return derivative_one_sided<Order, -1>(f, x, eps);
}
 
template <int Order, class F, class A>
auto derivative_backward(const F& f, A x, A eps = std::numeric_limits<A>::epsilon()) -> decltype(f(x))
{
   return derivative_one_sided<Order, +1>(f, x, eps);
}
 
template <int Order, class F, class A>
auto derivative_central(const F& f, A x, A eps = std::numeric_limits<A>::epsilon())
   -> decltype(f(x))
{
   static_assert(Order <= 3, "1st central derivative with order > 3 not implemented");
 
   using return_type = decltype(f(x));
 
   // coefficients from Math. Comp. 51 (1988), 699-706, Table 1
   // DOI: http://dx.doi.org/10.1090/S0025-5718-1988-0935077-0
   static const std::vector<std::vector<double> > coeffs = {
      {0.5},
      {2./3., -1./12.},
      {3./4., -3./20., 1./60.},
      {4./5., -1./5., 4./105., -1./280.}
   };
 
   const A h = std::fabs(x) < eps ? eps : std::sqrt(eps) * x;
   return_type result = 0;
 
   for (int i = 0; i < Order + 1; i++) {
      const double coeff = coeffs[Order][i];
      const A step = (i + 1) * h;
      result += coeff * (f(x + step) - f(x - step));
   }
 
   return result / h;
}
 
} // namespace flexiblesusy
 
#endif