What is the optimum epsilon/dx value to use within the finite difference method?

Question

double MyClass::dx = ?????;

double MyClass::f(double x)
{
    return 3.0*x*x*x - 2.0*x*x + x - 5.0;
}

double MyClass::fp(double x)    // derivative of f(x), that is f'(x)
{
    return (f(x + dx) - f(x)) / dx;
}

When using finite difference method for derivation, it is critical to choose an optimum dx value. Mathematically, dx must be as small as possible. However, I'm not sure if it is a correct choice to choose it the smallest positive double precision number (i.e.; 2.2250738585072014 x 10⁻³⁰⁸).

Is there an optimal numeric interval or exact value to choose a dx in to make the calculation error as small as possible?

(I'm using 64-bit compiler. I will run my program on a Intel i5 processor.)

Answer 1

Choosing the smallest possible value is almost certainly wrong: if dx were that smallest number, then f(x + dx) would be exactly equal to f(x) due to rounding.

So you have a tradeoff: Choose dx too small, and you lose precision to rounding errors. Choose it too large, and your result will be imprecise due to changes in the derivative as x changes.

To judge the numeric errors, consider (f(x + dx) - f(x))/f(x)¹ mathematically. The numerator denotes the difference you want to compute, but the denominator denotes the magnitude of numbers you're dealing with. If that fraction is about 2^‒k, then you can expect approximately k bits of precision in your result.

If you know your function, you can compute what error you'd get from choosing dx too large. You can then balence things, so that the error incurred from this is about the same as the error incurred from rounding. But if you know the function, you might be better off by providing a function that directly computes the derivative, like in your example with the polygonal f.

The Wikipedia section that pogorskiy pointed out suggests a value of sqrt(ε)x, or approximately 1.5e-8 * x. Without any more detailed knowledge about the function, such a rule of thumb will provide a reasonable default. Also note that that same section suggests not dividing by dx, but instead by (x + dx) - x, as this takes rounding errors incurred by computing x + dx into account. But I guess that whole article is full of suggestions you might use.

---

¹ This formula really should divide by f(x), not by dx, even though a past editor thought differently. I'm attempting to compare the amount of significant bits remaining after the division, not the slope of the tangent.

Answer 2

Why not just use the Power Rule to derive the derivative, you'll get an exact answer:

f(x) = 3x^3 - 2x^2 + x - 5
f'(x) = 9x^2 - 4x + 1

Therefore:

f(x) = 3.0 * x * x * x - 2.0 * x * x + x - 5.0
fp(x) = 9.0 * x * x - 4.0 * x + 1.0