deriv.h

Go to the documentation of this file.

 
 
#ifndef DERIVATION_THEORETICA_H
#define DERIVATION_THEORETICA_H
 
#include "../core/function.h"
#include "../polynomial/polynomial.h"
 
 
namespace theoretica {
 
 
    template<typename Field = real>
    inline polynomial<Field> deriv(const polynomial<Field>& p) {
 
        if (p.size() == 0)
            return polynomial<Field>();
 
        if (p.size() == 1)
            return polynomial<Field>({static_cast<Field>(0.0)});
 
        polynomial<Field> Dp;
        Dp.resize(p.size() - 1);
 
        for (unsigned int i = 1; i < p.size(); ++i)
            Dp[i - 1] = p[i] * i;
 
        return Dp;
    }
 
 
    inline real deriv(const polynomial<real>& p, real x) {
 
        real dp = 0.0;
 
        for (unsigned int i = 0; i + 1 < p.size(); ++i) {
 
            const unsigned int pos = p.size() - i - 1;
            dp = pos * p[pos] + x * dp;
        }
 
        return dp;
    }
 
 
    template <
        typename RealFunction = std::function<real(real)>,
        enable_real_func<RealFunction> = true
    >
    inline real deriv_central(RealFunction f, real x, real h = CALCULUS_DERIV_STEP) {
        return (f(x + h) - f(x - h)) / (2.0 * h);
    }
 
 
    template <
        typename RealFunction = std::function<real(real)>,
        enable_real_func<RealFunction> = true
    >
    inline real deriv_forward(RealFunction f, real x, real h = CALCULUS_DERIV_STEP) {
        return (f(x + h) - f(x)) / h;
    }
 
 
    template <
        typename RealFunction = std::function<real(real)>,
        enable_real_func<RealFunction> = true
    >
    inline real deriv_backward(RealFunction f, real x, real h = CALCULUS_DERIV_STEP) {
        return (f(x) - f(x - h)) / h;
    }
 
 
    template <
        typename RealFunction = std::function<real(real)>,
        enable_real_func<RealFunction> = true
    >
    inline real deriv_ridders2(RealFunction f, real x, real h = CALCULUS_DERIV_STEP) {
        return (4.0 * deriv_central(f, x, h / 2.0) - deriv_central(f, x, h)) / 3.0;
    }
 
 
    template <
        typename RealFunction = std::function<real(real)>,
        enable_real_func<RealFunction> = true
    >
    inline real deriv_ridders(RealFunction f, real x, real h = 0.01, unsigned int degree = 3) {
 
        const unsigned int MAX_RIDDERS_DEGREE = 16;
        real A[MAX_RIDDERS_DEGREE][MAX_RIDDERS_DEGREE];
 
        if (degree == 0 || degree > MAX_RIDDERS_DEGREE) {
            TH_MATH_ERROR("deriv_ridders", degree, MathError::InvalidArgument);
            return nan();
        }
 
        for (unsigned int i = 0; i < degree; ++i) {
            
            for (unsigned int n = 0; n <= i; ++n) {
                
                unsigned int m = i - n;
                
                if(n == 0) {
                    A[n][m] = deriv_central(f, x, h / (1 << m));
                } else {
                    real coeff = square(1 << n);
                    A[n][m] = (coeff * A[n - 1][m + 1] - A[n - 1][m]) / (coeff - 1);
                }
            }
 
        }
 
        return A[degree - 1][0];
    }
 
 
    template <
        typename RealFunction = std::function<real(real)>,
        enable_real_func<RealFunction> = true
    >
    inline real deriv(RealFunction f, real x, real h = CALCULUS_DERIV_STEP) {
        return deriv_ridders2(f, x, h);
    }
 
 
    template <
        typename RealFunction = std::function<real(real)>,
        enable_real_func<RealFunction> = true
    >
    inline real deriv2(RealFunction f, real x, real h = CALCULUS_DERIV_STEP) {
        return (f(x + h) - (2 * f(x)) + f(x - h)) / (h * h);
    }
 
}
 
 
#endif