keithalewis/epsilon

epsilon

This library can be used to compute derivatives of any order to machine precision.

If ε is not 0 and ε² = 0 then for any sufficiently differentiable function f, f(x + ε) = f(x) + f'(x) ε + f''(x) ε²/2 + · · · = f(x) + f'(x) ε. The coefficient of ε is the derivative of f.

For example, if f(x) = x² then (x + ε)² = x² + 2 x ε + ε² = x² + 2 x ε, so the derivative of x² is 2 x. No need to take limits of difference quotients.

Of course ε cannot be a real number, but the 2 x 2 matrix ε = [0 1; 0 0] satisfies these conditions. This is a special case of a Toeplitz matrix. The class epsilon implements the arithmetic datatype required to compute derivatives to any order.

In order to use this library you must define functions that take generic arguments. For example:

template<class X>
X f(X x)
{
    return x*x + 2*x + 3;
}

Now you can use epsilon to compute derivatives.

double x = 1.23;
auto y = f(epsilon<3>(x, 1));
assert (y[0] == f(x));
assert (y[1] == 2*x + 2); // f'(x)
assert (y[2] == 2);       // f''(x)

epsilon's Issues

Isnt this just forward mode automatic differentiation implemented with operator overloading on taylor series.