[Math] Derivation of fourth-order accurate formula for the second derivative

calculusfinite differencesnumerical methods

I am trying to derive / prove the fourth order accurate formula for the second derivative:

$f''(x) = \frac{-f(x + 2h) + 16f(x + h) – 30f(x) + 16f(x – h) – f(x -2h)}{12h^2}$.

I know that in order to do this I need to take some linear combination for the Taylor expansions of $f(x + 2h)$, $f(x + h)$, $f(x – h)$, $f(x -2h)$. For example, when deriving the the centered-difference formula for the first derivative, the Taylor expansion of $f(x + h)$ minus $f(x-h)$ can be computed to give the desired result of $f'(x)$, in that case.

In what way would I have to combine these Taylor expansions above to obtain the required result?

Best Answer

Exactly as Gammatester says, Taylor expand the terms upto order $4$ and verify. \begin{eqnarray*} -f(x+2h) &=& -f(x) &-& 2h f'(x) &-& 2h^2 f''(x) &-& \frac{4}{3} h^3 f'''(x) &-& \frac{2}{3} h^4 f''''(x) &+& O(h^5) \\ 16f(x+h) &=& 16 f(x)&+& 16h f'(x) &+& 8h^2 f''(x) &+& \frac{8}{3} h^3 f'''(x) &+& \frac{2}{3} h^4 f''''(x) &+& O(h^5) \\ -30f(x) &=& -30f(x) & & & & & & & & & & \\ 16f(x-h) &=& 16 f(x)&-& 16h f'(x) &+& 8h^2 f''(x) &-& \frac{8}{3} h^3 f'''(x) &+& \frac{2}{3} h^4 f''''(x) &+& O(h^5) \\ -f(x-2h) &=& -f(x) &+& 2h f'(x) &-& 2h^2 f''(x) &+& \frac{4}{3} h^3 f'''(x) &-& \frac{2}{3} h^4 f''''(x) &+& O(h^5) \\ \end{eqnarray*}

Related Solutions

Numerical Methods – Finite Difference Approximations Past Fourth Derivative

For even orders the finite-difference derivative approximation has a simple form.

Consider following finite-difference operator $\Delta$ $$ \Delta f(x) = \frac{f(x+h/2) - f(x-h/2)}{h}. $$ It's a second order first derivative operator. You can apply it several times $$ \Delta^2 f(x) = \Delta \Delta f(x) = \Delta \left(\frac{f(x+h/2) - f(x-h/2)}{h}\right) = \frac{f(x+h) - f(x)}{h^2} - \frac{f(x) - f(x-h)}{h^2}= \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} \sim f''(x). $$ Note that $\Delta$ can be written using shift operator $T_a = \exp\left(a \frac{d}{dx}\right)$ $$ \Delta = \frac{T_{h/2} - T_{-h/2}}{h} = \frac{e^{\frac{h}{2}\frac{d}{dx}} - e^{-\frac{h}{2}\frac{d}{dx}}}{h}. $$ Thus powers of $\Delta$ can be easily expressed using binomial theorem (note that $T_a$ and $T_b$ commute) $$ \Delta^{2k} = \frac{1}{h^{2k}}\sum_{m=0}^{2k} \begin{pmatrix}2k\\m\end{pmatrix} T_{h/2}^{2k-m} (-1)^m T_{-h/2}^{m} = \frac{1}{h^{2k}}\sum_{m=0}^{2k} \begin{pmatrix}2k\\m\end{pmatrix} (-1)^m T_{(2k-m)h/2} T_{-mh/2} = \\ =\frac{1}{h^{2k}}\sum_{m=0}^{2k} \begin{pmatrix}2k\\m\end{pmatrix} (-1)^m T_{(k-m)h} = \frac{1}{h^{2k}} \sum_{m=-k}^k \begin{pmatrix}2k\\k+m\end{pmatrix} (-1)^{m+k} T_{mh}. $$ Thus $$ \Delta^{2k} f(x) = \frac{1}{h^{2k}} \sum_{m=-k}^{k} \begin{pmatrix}2k\\k+m\end{pmatrix} (-1)^{m+k} f(x + mh). $$ This will be a second order formula, since it represents $2k$ times applying second order formula to $f(x)$. Actually, the truncation error is $$ \Delta^{2k} f(x) - f^{(2k)}(x) = \frac{kh^2}{12} f^{(2k+2)}(x) + O(h^4). $$

But from practical point of view, the formulas for high order derivatives (maybe more than 4) are almost never used. Let $f(x)$ be evaluated on a machine with limited number of binary digits. Let's assume double precision, which has $53$ bits of mantissa. Thus every $f(x)$ value would have also small error $f(x) \pm \epsilon |f(x)|$ with $\varepsilon =2^{-53} \approx 10^{-16}$. Applying $\Delta^{2k}$ to this function will yield $$ \Delta^{2k} (f(x) + \delta f(x)) = \Delta^{2k} f(x) + \Delta^{2k} \delta f(x) = f^{(2k)}(x) + \epsilon_{\text{trunc}} + \Delta^{2k} \delta f(x). $$ For the $\delta f$ it is known that $|\delta f(x)| \leq \varepsilon |f(x)|$, but the sign can be arbitrary. Thus the best we can write for $\Delta^{2k} \delta f(x)$ is $$ \left|\Delta^{2k} \delta f(x)\right| \leq \frac{1}{h^{2k}} \sum_{m=-k}^{k} \begin{pmatrix}2k\\k+m\end{pmatrix} |(-1)^{m+k}| \cdot |\delta f(x + mh)| \leq \frac{\epsilon \max |f(x)|}{h^{2k}} \sum_{m=-k}^{k} \begin{pmatrix}2k\\k+m\end{pmatrix} = \\ = 2^{2k}\frac{\varepsilon \max |f(x)|}{h^{2k}}. $$ Note that the bound is tight, it's not hard to construct such $\delta f(x)$. Since $\epsilon_{\text{trunc}} \leq \frac{kh^2}{12}\max |f^{(2k+2)|}(x)$ the total error consists of two terms $$ \epsilon \leq \frac{kh^2}{12} M_{2k+2} + 2^{2k}\frac{\varepsilon M_0}{h^{2k}}. $$ Here for simplicity $M_k = \max |f^(k)(x)|$. From this estimate one can find the optimal $h_\star$ and minimal achievable $\epsilon_\star$: $$ 0 = \frac{d\epsilon}{dh} = \frac{kh}{6}M_{2k+2} - 2^{2k+1}k \frac{\varepsilon M_0}{h^{2k+1}}\\ h_\star^{2k+2} = \frac{3\varepsilon 2^{2k+2}M_0}{M_{2k+2}}, \qquad h_\star = 2 \sqrt[2k+2]\frac{3\varepsilon M_0}{M_{2k+2}}\\ \epsilon_\star = \frac{kh_\star ^2}{12} M_{2k+2} \left( 1 + 2^{2k+2}\frac{\varepsilon M_0}{h_\star^{2k+2}}\frac{3}{k M_{2k+2}} \right) = \frac{kh_\star ^2}{12} M_{2k+2} \frac{k+1}{k} = \frac{(k+1)h_\star ^2}{12} M_{2k+2} $$ Here's a table for different $k$ and $M_0 = M_{2k+2} = 1$ ($f(x) = \sin x$ for example), $\varepsilon = 2^{-53}$: $$ \begin{array}{ccc} k & h_\star & \epsilon_\star\\ \hline 1 & 2.70186\times 10^{-4} & 1.21667\times 10^{-8} \\ 2 & 5.26565\times 10^{-3} & 6.93176\times 10^{-6} \\ 3 & 2.32459\times 10^{-2} & 1.80124\times 10^{-4} \\ 4 & 5.66609\times 10^{-2} & 1.33769\times 10^{-3} \\ 5 & 1.02622\times 10^{-1} & 5.26565\times 10^{-3} \\ 6 & 1.56854\times 10^{-1} & 1.43519\times 10^{-2} \\ 7 & 2.1562\times 10^{-1} & 3.09945\times 10^{-2} \\ 8 & 2.76166\times 10^{-1} & 5.72009\times 10^{-2} \\ 9 & 3.36633\times 10^{-1} & 9.44348\times 10^{-2} \\ 10 & 3.9583\times 10^{-1} & 1.43625\times 10^{-1} \\ \end{array} $$ Here's also results from numerical experiments with $k=2$ (fourth order) and $k=6$ (12th order)

Note that 12th order formula hardly gave two significant digits of the derivative, while 4th order gave 6. Optimal values agree with theoretical values pretty well. To get at least 5 right digits for the 12th derivative you need double-double (128 bytes) floating point arithmetic.

Note that if you need such finite difference for some ODE, for example $$ y^{(12)}(x) + a y^{(6)}(x) + b y'(x) + c y(x) = f(x) $$ you can always rewrite it as a system of the first order with 12 unknows $$ y_{11}'(x) + a y_6(x) + b y_1(x) + c y_0(x) = f(x)\\ y_{10}'(x) = y_{11}(x)\\ \vdots\\ y_{0}'(x) = y_1(x). $$

Fourth order accurate approx

No, your expression evaluates to $12h^2\color{red}{f''(x_0)+O(h^4)}$, not $12h^2$, which doesn't give you the error estimate $+O(h^4)$ (only $+O(h^2)$) after you divide by the denominator $12h^2$.

To get the explicit $O(h^4)$ term means you need to expand to order (at least) 6. To save space, I'll note that since the numerator $$S:=-f(x_0+2h)+16 f(x_0+h)-30 f(x_0)+16 f(x_0-h)-f(x_0-2h)$$ is even in $h$ only the even-order derivatives appear. Writing $[h^n]F$ for the coefficient of $h^n$ in the Taylor series expansion of a (sufficiently smooth) function $F(h)$, we have \begin{align} [h^n](f(x_0+mh))&=[h^n]\left(\sum_{j=0}^n f^{(j)}(x_0)\frac{(mh)^j}{j!}+o(h^n)\right)\\ &=\frac{f^{(n)}(x_0)}{n!}m^n \end{align} so $$ \boxed{\color{blue}{[h^n]S=\frac{f^{(n)}(x_0)}{n!}\bigg[-(2)^n+16(1)^n-30(0)^n+16(-1)^n-(-2)^n\bigg]}} $$ For even $n>0$, $$ -(2)^n+16(1)^n-30(0)^n+16(-1)^n-(-2)^n =-2(2^n-16) $$ So we can write down the coefficients immediately \begin{align} [1]S&=0\\ [h^2]S&=\frac{-2(2^2-16)}{2!}f^{(2)}(x_0)=12f''(x_0)\\ [h^4]S&=0\\ [h^6]S&=\frac{-2(2^6-16)}{6!}f^{(6)}(x_0)=-\frac{2}{15}f^{(6)}(x_0). \end{align} Hence $$ S=12f''(x_0)h^2-\frac{2}{15}h^6f^{(6)}(x_0)+o(h^6), $$ which rearranges to $$ f''(x_0)=\frac{S}{12h^2}+\left(\frac{1}{90}f^{(6)}(x_0)+o(1)\right)h^4. $$ I.e., the leading order error term is $\frac{1}{90}f^{(6)}(x_0)h^4$.

Best Answer

Related Solutions

Numerical Methods – Finite Difference Approximations Past Fourth Derivative

Fourth order accurate approx

Related Question