Optimum small number for numerical differentiation 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?
 A: 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.)
A: 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$. 
A: 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.
