For higher-order derivatives, we proceed similarly to first-order derivatives by differentiating the interpolation polynomial as many times as needed. A major difficulty arises in error analysis since differentiating the error term multiple times is tedious. A Taylor series approach is often preferred.

Using the previously considered degree 2 polynomial:

p₂(x) = f(x₀) + f[x₀,x₁](x − x₀) + f[x₀,x₁,x₂](x − x₀)(x − x₁)
p″₂(x) = 2 f[x₀,x₁,x₂] = (f(x₂) − 2f(x₁) + f(x₀)) / h²

This approximates the second derivative f″(x) in [x₀, x₂]. Its order depends on the point of approximation.

Case 1 – Approximation at x₀:

f″(x₀) ≃ p″₂(x₀) = (f(x₀ + 2h) − 2f(x₀ + h) + f(x₀)) / h²

Using the Taylor expansion, we get: f″(x₀) + O. This is a first-order forward difference.

Case 2 – Approximation at x₁:

f″(x₁) ≃ p″₂(x₁) = (f(x₁ + h) − 2f(x₁) + f(x₁ − h)) / h²

This is a centered difference. Taylor expansion shows the error is O(h²), giving a second-order approximation.

Case 3 – Approximation at x₂:

f″(x₂) ≃ p″₂(x₂)

Similar reasoning to case 1 shows this is a first-order backward difference.

Remark: The same degree 2 polynomial gives first-order differences at x₀ and x₂, and a second-order difference at x₁ due to the symmetry of centered differences.

Finite Difference Formulas for Second Derivatives

f″(x) = (f(x−2h) − 2f(x−h) + f(x)) / h² + O (backward)
f″(x) = (f(x+2h) − 2f(x+h) + f(x)) / h² + O (forward)
f″(x) = (f(x+h) − 2f(x) + f(x−h)) / h² + O(h²) (centered)
f″(x) = (−f(x+2h) + 16f(x+h) − 30f(x) + 16f(x−h) − f(x−2h)) / 12h² + O(h⁴) (centered)

f⁽⁴⁾(x) = (f(x+2h) − 4f(x+h) + 6f(x) − 4f(x−h) + f(x−2h)) / h⁴ + O(h²) (centered)

Differentiation is numerically unstable. Finite difference formulas depend on the step h. Reducing h does not always improve accuracy. In first derivative centered differences:

f′(x₀) ≃ (f(x₀ + h) − f(x₀ − h)) / 2h

When h → 0, subtraction of nearly equal numbers causes loss of significant digits, combined with division by a small number, illustrating the limitations of numerical differentiation.