Optimum small number for numerical differentiation

Question

The Wikipedia article on numerical differentiation mentions the formula

$$ h=\sqrt \epsilon \times x $$

where $\epsilon$ is the machine epsilon (approx. $2.2\times 10^{-16}$ for 64-bit IEEE 754 doubles), to calculate the optimum "small number" $h$ to be used in differentiation, such as $$ \frac{f(x+h)-f(x)}{h} $$ But what if $x$ is zero? Then $h$ will be zero too, and division by zero is certainly not a way to do numerical differentiation. Is the article wrong? Is it otherwise correct, except that near zero (how near?) some small enough constant (how small?) should be used?

Answer 1

If you click through to the reference given for the Wikipedia piece, you'll find an answer. The formula given there is to take $h$ to be roughly $\sqrt{\epsilon_f}x_c$, where $\epsilon_f$ isn't necessarily "machine epsilon," but more to the point, where $x_c$ isn't necessarily $x$.

Answer 2

What you want is a number $\epsilon$ that is small enough that it isn't zero, but its square is zero. Because then we have $f(x+\epsilon)=f(x)+\epsilon f'(x)+0$, for suitable $f$. No real number has this property, but if you're implementing a numerical method on a computer it's straightforward to implement a type containing numbers that do have this property. It'll probably give you better results than anything involving small real numbers.

(BTW There have also been papers published that propose using a complex $\epsilon$ but I think these are misguided.)

Answer 3

I use the technique given by J.C. Nash in the book Compact Numerical Methods for Computers.

On page 219, there is the formula

$h=\sqrt{\epsilon}\left(|x|+\sqrt{\epsilon}\right)$

where $\epsilon$ is of course machine epsilon. I have also seen the following formula used (apologies, I no longer recall where this was from):

$h=\sqrt{\epsilon}\max\left(|x|,\sqrt{\epsilon}\right)$

The recipe for second derivatives is similar except that one now uses the cube root of machine epsilon instead of its square root.