Question

http://en.wikipedia.org/wiki/Numerical_differentiation#Practical_considerations mentions the formula

$h=\sqrt \epsilon * x$

where $\epsilon$ is the machine epsilon (approx. 2.2e-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

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.

Answer 3

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.)