# Is there a general formula for estimating the step size h in numerical differentiation formulas?

Using three-point central-difference formula

$$f^{\prime}(x_0)\approx \frac{f(x_0+h)-f(x_0-h)}{2h}$$

and for $f(x)=\exp(x)$ at $x_0=0$ we have

$$\begin{array}{c, l, r} h & f^{\prime}(0) & error \\ \hline 10^{-01} & 1.0017 & 1.6675\times 10^{-03} \\ 10^{-02} & 1 & 1.6667\times 10^{-05} \\ 10^{-03} & 1 & 1.6667\times 10^{-07} \\ 10^{-04} & 1 & 1.6669\times 10^{-09} \\ 10^{-05} & 1 & 1.2102\times 10^{-11} \\ 10^{-06} & 1 & -2.6755\times 10^{-11} \\ 10^{-07} & 1 & -5.2636\times 10^{-10} \\ 10^{-08} & 1 & -6.0775\times 10^{-09} \\ 10^{-09} & 1 & 2.7229\times 10^{-08} \\ 10^{-10} & 1 & 8.2740\times 10^{-08} \\ 10^{-11} & 1 & 8.2740\times 10^{-08} \\ 10^{-12} & 1 & 3.3389\times 10^{-05} \\ 10^{-13} & 9.9976\times 10^{-01} & -2.4417\times 10^{-04} \\ 10^{-14} & 9.9920\times 10^{-01} & -7.9928\times 10^{-04} \\ 10^{-15} & 1.0547 & 5.4712\times 10^{-02} \\ 10^{-16} & 5.5511\times 10^{-01} & -4.4489\times 10^{-01} \\ \end{array}$$

From $10^{-1}$ down to $10^{-5}$ the results are evident (because the rate of convergence of the three-point central-difference formula is $O(h^2)$). As you see because of the round-off error, the error deteriorate rapidly as $h$ decrease. My question is: Is there a general formula for estimating the step size $h$ in numerical differentiation formulas to get the best result?

• Usually you know how much error you can tolerate and you use the error formula to compute the value of h that will suffice – bobbym May 30 '14 at 17:02
• It depends rather how you are computing - your formula is $\frac {\sinh h}h=1+\frac {h^2}{3!}+\frac {h^4}{5!}+\dots$ which I can compute pretty accurately by hand without the kind of problem your method has. So the error is an artefact of the computation rather than the particular problem. – Mark Bennet May 30 '14 at 17:04
• Note that the error is expected to be close to $\frac {h^2}6$ so your method is drifting off already at $10^{-5}$ – Mark Bennet May 30 '14 at 17:28
• @bobbym Any tips on how to set up that inequality? – Vaughan Hilts Mar 22 '15 at 1:33

Well, there is a method to estimate the optimal $h$ when a roundoff error is present. I'll show that for the central difference formula at hand, but you can repeat the same process for any similar method: Assume that for all $i \in \Bbb N$ holds : $|f^{(i)}| < M_i$. The function evaluations we have contain a roundoff error, denote it $\delta(x)$. Let us assume that $|\delta(x)| < \epsilon$. Denote the measured function value at $x_0$ as $\bar{f}(x_0)$ and the actual value as $f(x_0)$.

Then, we estimate the error

$$error(h) = \left|f'(x_0) - \frac{\bar{f}(x_0+h) - \bar{f}(x_0-h)}{2h}\right| = \left|f'(x_0) - \frac{f(x_0+h) - f(x_0-h)}{2h} + \frac{\delta(x_0+h) - \delta(x_0-h)}{2h}\right| \le \left|f'(x_0) - \frac{f(x_0+h) - f(x_0-h)}{2h}\right| + \left|\frac{\delta(x_0+h) - \delta(x_0-h)}{2h}\right| \le \frac{h^2M_3}{12} + \frac{2\epsilon}{2h} = \frac{h^2M_3}{12} + \frac{\epsilon}{h} := g(h)$$

Now, you can differentiate $g$ with respect to $h$ and find the $h$ that minimizes it.

In this case we get that $h_{opt} = \sqrt[3]{\frac{6\epsilon}{M_3}} \approx 10^{-6}$, assuming $\epsilon = 10^{-16}$ and a reasonable number for $M_3$.

• How can I determine $M_3$? – Dante Jun 3 '14 at 15:32
• There are no definitive rules. In your case, $f = e^x$, which means that $f^{(3)} = e^x$, and you are interested in a derivative around $0$. So, for example, your calculations never take outside the interval[-1,1] (A very overkill estimate), so you can say that $|f^{(i)}| < e$ for all $i$. Also, the roundoff error is less additive, than relative. What I want to say, that in the calculations above, it would be more proper to say that $\left| \frac{\bar{f} - f}{f}\right| \le \epsilon$ – Aahz Jun 3 '14 at 17:26
• (Sorry, a bug with the site. This comment is a firect continuation of the previous one. If someone can merge them, that would be nice.) But that only will add an $M_0$ factor in front of the $\epsilon$ in the final statement - that does not affect much. In the end $g(h) = \frac{h^2M_3}{12} + \frac{\epsilon M_0}{h}$, and $h_{opt} = \sqrt[3]{\frac{6 \epsilon M_0}{M_3}}$ But note the calculations in the original post still hold, the difference is only the estimate of the error, and even then the end result only differs my a constant. – Aahz Jun 3 '14 at 17:37

Rule of thumb: take the relative step $\frac h{x_0}$ to be the square root of the machine ulp. $$h=x_0\sqrt{ulp}$$ The rationale is that the truncation and roundoff errors are then of the same order.

When $x_0=0$, no rule :(

• If I am reading that wiki correctly, that would imply ulp for current double precision is 2^-53? – dashnick Jan 20 '17 at 2:52
• @dashnick: exactly, as there are 53 bits in the representation of the number. This is a little more than $10^{-16}$, giving hope for derivatives with $8$ accurate digits. – Yves Daoust Jan 20 '17 at 7:35