4th order accurate difference formula less accurate than 2nd order formula?

I understand that there are different difference expressions for calculating numerical derivatives. For example, this one is a central difference formula supposed to be 2nd order accurate, i.e. error decreases as O(h^2)

$\frac {dy}{dx} = \frac {y(x+h)-y(x-h)}{2} $

This formula is supposed to be 4th order accurate, i.e. error decreases as O(h^4)

$\frac {dy}{dx} = \frac {y(x-2h)-8y(x-h)+8y(x+h)-y(x+2h)}{12} $

Is the second formula always better than the first formula or not ? I have heard someone say that if increment h is not small enough 2nd formula will be more inaccurate than first one. If this is the case, is there a way to figure out the limit of h below which 2nd formula will give more accurate values of derivatives compared to the 1st formula ? Thanks.

Solution 1:

Your formula seems of. Set $$ Df(x,h)=\frac{f(x+h)-f(x-h)}{2h} $$ then $$ Df(x,h)=f'(x)+\frac{f'''(x)}{6}h^2+O(h^4) $$ To eliminate the quadratic term, you have to compute $$ 4Df(x,h)-f(x,2h)=3f'(x)+O(h^4) $$ which indeed results in $$ \iff f'(x)=\frac{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}{12h}+O(h^4) $$ This is also called Richardson extrapolation (more commonly known in the context of Romberg integration). Naming the expression $D_2f(x,h)$, the next interpolation step would be $$ 16D_2f(x,h)-D_2f(x,2h)=15f'(x)+O(h^6), $$ but one can also obtain that error order by combining the central differences for $h,2h,3h$.

Now to the errors: The function evaluation will always introduce some noise of the relative magnitude $\mu$, where $\mu$ is the machine constant of the floating point type. Thus your error is a combination of the approximation error and that noise, or about $$\frac{M_0\mu}h+M_3h^2$$ for the simple central difference which is smallest when both terms are about equal, which in most cases gives $h\sim\sqrt[3]\mu$ if the magnitudes $M_0$ of the function values and $M_3$ of the third derivative (and below $M_5$ of the fifth derivative) are about equal. For double this gives $h\sim 10^{-5}$ reducing the (relative) error to about $10^{-10}$ as a rather hard lower bound.

For the Richardson extrapolation formula, one gets in the same way the error about $$\frac{M_0\mu}h+M_5h^4,$$ which is minimal at about $h\sim\sqrt[5]\mu$, for double this is $h\sim 10^{-3}$, giving about $10^{-12}$ as the best relative error.

The following graph demonstrates the errors of the one-sided and central difference quotients and the first Richardson extrapolation.

To answer the question, at $h=10^{-5}$ where the simple central difference gives best results with error $10^{-10}$, the extrapolation formula will be dominated by the noise term $\mu/h$, which is also of size $10^{-10}$. The actual constants involved in the error terms are different, so that the second error might be slightly larger than the first.