This library can be used to compute derivatives of any order to machine precision.
If ε is not 0 and ε2 = 0 then for any sufficiently differentiable function f, f(x + ε) = f(x) + f'(x) ε + f''(x) ε2/2 + · · · = f(x) + f'(x) ε. The coefficient of ε is the derivative of f.
For example, if f(x) = x2 then (x + ε)2 = x2 + 2 x ε + ε2 = x2 + 2 x ε, so the derivative of x2 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)