Improve the order of accuracy for second derivative using central difference approximation

Question

The formula $$D^{m}(h)=(4^mD^{m-1}(h/2)-D^{m-1}(h))/3$$ for $m=2,3,\ldots$improves the order of accuracy for first derivative using central diff approx. where $$D^{1}(h)=f'(x)=(f(x+h)-f(x-h))/(2*h)$$ is of the order $O(h^2)$. Using the above formula once the order of accuracy of 1st derivative will improve from $O(h^2)$ to $O(h^4)$. How can we improve the order of accuracy for 2nd derivative $$f''(x)=(f(x+h)-2f(x)+f(x-h))/(h^2)?$$

Answer 1

As a first step, consider first the Taylor expansions of $$f(x-2h), \: f(x-h), \: f(x+h), \text{ and } f(x+2h)$$ which are: \begin{align} f(x-2 h)&\approx\frac{2}{3} h^4 f^{(4)}(x)-\frac{4}{3} h^3 f^{(3)}(x)+2 h^2 f''(x)-2 h f'(x)+f(x)\\ f(x-h)&\approx\frac{1}{24} h^4 f^{(4)}(x)-\frac{1}{6} h^3 f^{(3)}(x)+\frac{1}{2} h^2 f''(x)-h f'(x)+f(x)\\ f(x+h)&\approx\frac{1}{24} h^4 f^{(4)}(x)+\frac{1}{6} h^3 f^{(3)}(x)+\frac{1}{2} h^2 f''(x)+h f'(x)+f(x)\\ f(x+2 h)&\approx\frac{2}{3} h^4 f^{(4)}(x)+\frac{4}{3} h^3 f^{(3)}(x)+2 h^2 f''(x)+2 h f'(x)+f(x) \end{align}

Those (approximate) equations may be solved for the derivatives to obtain:

\begin{align} f'(x)&\approx \frac{f(x-2 h)-8 f(x-h)+8 f(h+x)-f(2 h+x)}{12 h}\\ f''(x)&\approx \frac{-f(x-2 h)+16 f(x-h)+16 f(h+x)-f(2 h+x)-30 f(x)}{12 h^2}\\ f^{(3)}(x)&\approx \frac{-f(x-2 h)+2 f(x-h)-2 f(h+x)+f(2 h+x)}{2 h^3}\\ f^{(4)}(x)&\approx \frac{f(x-2 h)-4 f(x-h)-4 f(h+x)+f(2 h+x)+6 f(x)}{h^4} \end{align}

If you consider that the Taylor expansions are $O(h^5)$ and carry that through, you can even work out associated error estimates for the difference quotients.

Answer 2

This is called Richardson extrapolation.

The property used in both the first and second derivative is the invariance under time reversion due to the symmetry of the formula. This gives that the error is an even function in $h$, thus the same extrapolation formula applies for the second as for the first derivative, i.e., $$ D^{m+1}(h)=(4^mD^{m}(h/2)−D^{m}(h))/(4^m-1). $$ in notable difference to your formula.

---

The derivation of this formula follows simply from the assumption $D^m(h)-D^*=O(h^{2m})$ which gets quantified to the next order as $$ D^m(h)=D^*+C·h^{2m}+O(h^{2(m+1)}) $$ where $D^*$ is the exact value. The leading error term can now be eliminated as $$ 4^mD^m(h/2)-D^m(h)=(4^m-1)D^*+O(h^{2(m+1)}). $$