This is a Numerical Analysis textbook question on Numerical Differentiation (like https://en.wikipedia.org/wiki/Numerical_differentiation)
(fyi, I asked a similar but different question on the same textbook section earlier. This is the last one, promise 🙂 )
Official Textbook Question
Derive a method for approximating $f'''(x_0)$ whose error term is of order $h^2$ by expanding the function $f$ in a fourth Taylor polynomial about $x_0$ and evaluating at $x_0 \pm h$ and $x_0 \pm 2h$
My Question
I believe I can do this problem, but I end up with an error term that has a $h^4$ term while the question asks for a solution with an error term of order $h^2$. Is this right, am I doing something wrong?
My Work
First basic four term Taylor expansion and remainder term:
\begin{align*}
f(x) &= f(x_0) + \frac{f'(x_0)}{1!} (x – x_0) + \frac{f''(x_0)}{2!} (x – x_0)^2 + \frac{f'''(x_0)}{3!} (x – x_0)^3 + \frac{f''''(\xi(x))}{4!} (x – x_0)^4 \\
f(x) &= f(x_0) + f'(x_0) (x – x_0) + \frac{1}{2} f''(x_0) (x – x_0)^2 + \frac{1}{6} f'''(x_0) (x – x_0)^3 + \frac{1}{24} f''''(\xi(x)) (x – x_0)^4 \\
\end{align*}
Apply for $x_0 \pm h, x_0 \pm 2h$:
\begin{align*}
f(x_0 – 2h) &= f(x_0) – 2 f'(x_0) h + 2 f''(x_0) h^2 – \frac{4}{3} f'''(x_0) h^3 + \frac{2}{3} f''''(\xi(x_0 – 2h)) h^4 \\
f(x_0 – h) &= f(x_0) – f'(x_0) h + \frac{1}{2} f''(x_0) h^2 – \frac{1}{6} f'''(x_0) h^3 + \frac{1}{24} f''''(\xi(x_0 – h)) h^4 \\
f(x_0 + h) &= f(x_0) + f'(x_0) h + \frac{1}{2} f''(x_0) h^2 + \frac{1}{6} f'''(x_0) h^3 + \frac{1}{24} f''''(\xi(x_0 + h)) h^4 \\
f(x_0 + 2h) &= f(x_0) + 2 f'(x_0) h + 2 f''(x_0) h^2 + \frac{4}{3} f'''(x_0) h^3 + \frac{2}{3} f''''(\xi(x_0 + 2h)) h^4 \\
\end{align*}
Dropping the remainder terms and converting to matrix form:
\begin{align*}
\begin{bmatrix}
1 & -2 & 2 & -\frac{4}{3} \\
1 & -1 & \frac{1}{2} & -\frac{1}{6} \\
1 & 1 & \frac{1}{2} & \frac{1}{6} \\
1 & 2 & 2 & \frac{4}{3} \\
\end{bmatrix}
\begin{bmatrix}
f(x_0) \\
h \cdot f'(x_0) \\
h^2 \cdot f''(x_0) \\
h^3 \cdot f'''(x_0) \\
\end{bmatrix}
&\approx
\begin{bmatrix}
f(x_0 – 2h) \\
f(x_0 – h) \\
f(x_0 + h) \\
f(x_0 + 2h) \\
\end{bmatrix} \\
\end{align*}
Inverting the matrix:
Wolfram Alpha
inverse {{1, -2, 2, -4/3}, {1, -1, 1/2, -1/6}, {1, 1, 1/2, 1/6}, {1, 2, 2, 4/3}}
\begin{align*}
\begin{bmatrix}
-\frac{1}{6} & \frac{2}{3} & \frac{2}{3} & -\frac{1}{6} \\
\frac{1}{12} & -\frac{2}{3} & \frac{2}{3} & -\frac{1}{12} \\
\frac{1}{3} & -\frac{1}{3} & -\frac{1}{3} & \frac{1}{3} \\
-\frac{1}{2} & 1 & -1 & \frac{1}{2} \\
\end{bmatrix}
\begin{bmatrix}
f(x_0 – 2h) \\
f(x_0 – h) \\
f(x_0 + h) \\
f(x_0 + 2h) \\
\end{bmatrix}
&\approx
\begin{bmatrix}
f(x_0) \\
h \cdot f'(x_0) \\
h^2 \cdot f''(x_0) \\
h^3 \cdot f'''(x_0) \\
\end{bmatrix}
\end{align*}
This gives us:
\begin{align*}
f'''(x_0) &\approx \frac{1}{h^3} \left[ -\frac{1}{2} f(x_0 – 2h) + f(x_0 – h) – f(x_0 + h) + \frac{1}{2} f(x_0 + 2h) \right] \\
\end{align*}
This gives an error term of:
\begin{align*}
h \left[ -\frac{1}{3} f''''(\xi(x_0 – 2h))
+ \frac{1}{24} f''''(\xi(x_0 – h))
– \frac{1}{24} f''''(\xi(x_0 + h))
+ \frac{1}{3} f''''(\xi(x_0 + 2h)) \right] \\
\end{align*}
Is this wrong because my error term has $h$ instead of $h^2$?
UPDATE: Fixed issue where I was off by a factor of $h^3$
Best Answer
You are missing something in establishing your system. In extracting the rational constants, the powers of $h$ are still bound to the derivatives, thus the system has to be \begin{align*} \begin{bmatrix} 1 & -2 & 2 & -\frac{4}{3} \\ 1 & -1 & \frac{1}{2} & -\frac{1}{6} \\ 1 & 1 & \frac{1}{2} & \frac{1}{6} \\ 1 & 2 & 2 & \frac{4}{3} \\ \end{bmatrix} \begin{bmatrix} f(x_0) \\ f'(x_0)h \\ f''(x_0)h^2 \\ f'''(x_0)h^3 \\ \end{bmatrix} &\approx \begin{bmatrix} f(x_0 - 2h) \\ f(x_0 - h) \\ f(x_0 + h) \\ f(x_0 + 2h) \\ \end{bmatrix} \\ \end{align*} with the solution formula \begin{align*} \begin{bmatrix} -\frac{1}{6} & \frac{2}{3} & \frac{2}{3} & -\frac{1}{6} \\ \frac{1}{12} & -\frac{2}{3} & \frac{2}{3} & -\frac{1}{12} \\ \frac{1}{3} & -\frac{1}{3} & -\frac{1}{3} & \frac{1}{3} \\ -\frac{1}{2} & 1 & -1 & \frac{1}{2} \\ \end{bmatrix} \begin{bmatrix} f(x_0 - 2h) \\ f(x_0 - h) \\ f(x_0 + h) \\ f(x_0 + 2h) \\ \end{bmatrix} &\approx \begin{bmatrix} f(x_0) \\ f'(x_0)h \\ f''(x_0)h^2 \\ f'''(x_0)h^3 \\ \end{bmatrix} \\ \end{align*} You should also recognize that because of the symmetry, the coefficients of the 4th derivative cancel so that the error term for $f'''(x_0)h^3$ is $O(h^5)$. As you have to divide by $h^3$, the error term for $f'''(x_0)$ is $O(h^2)$.
You can greatly reduce the computations by recognizing the symmetry of the situation and start with the Taylor expansions of $$ f(x_0+kh)-f(x_0-kh)=2f'(x_0)kh+\frac13f'''(x_0)k^3h^3+\frac1{60}f^{(5)}(x_0)k^5h^5+O(h^7) $$ so that you get for the coefficients $$2A_2+A_1=0\\8A_2+A_1=3$$