optimal way to approximate second derivative

Question

Suppose there is a function $f: \mathbb R\to \mathbb R$ and that we only know $f(0),f(h),f'(h),f(2h)$ for some $h>0$. and we can't know the value of $f$ with $100$% accuracy at any other point.

What is the optimal way of approximating $f''(0)$ with the given data?

I'd say that $f''(0)=\frac{f'(h)-f'(0)}{h}+O(h)$ and $f'(0)=\frac{f(h)-f(0)}{h}+O(h)$, therefor we get

$$f''(0)=\frac{f'(h)-\frac{f(h)-f(0)}{h}}{h}+O(h)$$

But that can't be the optimal way since we know $f(2h)$ and i didn't use it at all.

Could someone shed some light?

Answer 1

Let's say that $$y_k=\{f(0),f(h),f'(h),f(2h)\}$$ are given and $$x_k=\{f(0),f'(0),f''(0),f'''(0)\}$$ are unknown. Taylor's theorem gives a way of writing each $y_k$ as a linear combination of $x_j$'s, dropping all the terms starting with $f^{(4)}(0)$. For example: $$f(2h)=f(0)+f'(0)(2h)+\frac12f''(0)(2h)^2+\frac16f'''(0)(2h)^3. $$

Conversely, let $y=Ax$ be the linear equations relating $y$'s and $x$'s. These are four linear equations in four unknowns. To find the solution $x_2=f''(0)$ in terms of $y$'s just solve the system of linear equations for $x$. To find how accurate the approximation is, compute the Taylor expansion of $y$'s to a higher order and substitute into the expression for $x_2$. The answer will always have the form $$ f''(0) = \alpha_1 y_1+\cdots+\alpha_4 y_4 + \Theta(h^m), $$ where $m$ is the order of the approximation.

This approximation is "best" in the sense that it is the highest-order approximation possible (in this case $x_2 = f''(0) + O(h^2)$) with the given data. This is also the way all the standard finite-difference formulas are derived.

Answer 2

I believe the answer is something along the lines of: $$f''(x) \approx \frac{f(x+2h)-2f(x+h)+f(x)}{h^2}$$ For $x=0$, this becomes: $$f''(0) \approx \frac{f(2h)-2f(h)+f(0)}{h^2}$$

Source: Wikipedia finite Difference.

Answer 3

Because you only have four pieces of information - your finite difference scheme is essentially fitting a cubic. Additionally, because most of your information is centered around $x=h$ I would build the interpolating polynomial around that point and fit it with the information at hand: $$f(h+\delta) \approx a_0 + a_1\delta + a_2\delta^2 + a_3\delta^3 +O(\delta^4)$$ then implementing the information that you are given: $$f(0) = a_0 - a_1h + a_2h^2 - a_3h^3\\ f(h) = a_0\\ f(2h) = a_0 + a_1h + a_2h^2 + a_3h^3\\ f'(h) = a_1\\ $$

which gives a system of equations: $$ \begin{pmatrix} 1 & -h & h^2 & -h^3 \\ 1 & 0 & 0 & 0 \\ 1 & h & h^2 & h^3 \\ 0 & 1 & 0 & 0 \\ \end{pmatrix}\begin{pmatrix}a_0\\a_1\\a_2\\a_3\end{pmatrix}=\begin{pmatrix}f(0)\\f(h)\\f(2h)\\f'(h)\end{pmatrix} $$

after solving for the $a_k$ (there are really only two unknowns here)

I get $$a_2 = \frac{f(0)-2f(h)+f(2h)}{2h^2}$$ and $$a_3 = \frac{f(2h)-2f'(h)h - f(0)}{2h^3}$$

then you can simply compute $f''(0) = 2a_2 - 6a_3 h$.

This solution will minimize the energy of the error.