Module 7: Numerical Differentiation

1. Introduction

Numerical differentiation is the process of estimating the derivative of a mathematical function using values of the function. In many real-world physics and engineering applications, an analytical expression for a function may not be available. Instead, we might only have discrete data points from an experiment or simulation, such as position measurements over time, from which we need to derive velocity and acceleration. Alternatively, the function might be known but analytically intractable to differentiate. Numerical differentiation provides robust methods for approximating these derivatives and understanding the intrinsic trade-offs between mathematical truncation errors and computational round-off errors.

2. Theory

The foundation of most numerical differentiation formulas is the Taylor series expansion. For a sufficiently smooth function \(f(x)\), its value near a point \(x\) can be expanded as:

\[f(x+h) = f(x) + hf'(x) + \frac{h^2}{2!}f''(x) + \frac{h^3}{3!}f'''(x) + \mathcal{O}(h^4)\]

where \(h\) is a small step size.

Finite Difference Approximations

Forward Difference

By rearranging the Taylor series and solving for \(f'(x)\), we get the forward difference approximation:

\[f'(x) = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f''(x) + \dots\] \[f'(x) \approx \frac{f(x+h) - f(x)}{h}\]

The leading error term is proportional to \(h\), so we say this method is \(\mathcal{O}(h)\) accurate.

Backward Difference

Similarly, using the expansion for \(f(x-h)\): \[f(x-h) = f(x) - hf'(x) + \frac{h^2}{2!}f''(x) - \frac{h^3}{3!}f'''(x) + \mathcal{O}(h^4)\] \[f'(x) \approx \frac{f(x) - f(x-h)}{h}\] This is also an \(\mathcal{O}(h)\) approximation.

Central Difference

If we subtract the \(f(x-h)\) expansion from the \(f(x+h)\) expansion, the even terms cancel out: \[f(x+h) - f(x-h) = 2hf'(x) + \frac{2h^3}{3!}f'''(x) + \dots\] \[f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}\] The leading error term is proportional to \(h^2\), making the central difference an \(\mathcal{O}(h^2)\) method, which is significantly more accurate than forward or backward differences for small \(h\).

Higher-Order Derivatives

To find the second derivative, we can add the Taylor expansions for \(f(x+h)\) and \(f(x-h)\): \[f(x+h) + f(x-h) = 2f(x) + h^2 f''(x) + \mathcal{O}(h^4)\] \[f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}\] This is an \(\mathcal{O}(h^2)\) approximation for the second derivative.

Richardson Extrapolation

Richardson extrapolation is a technique to improve the accuracy of a numerical approximation. If we have an \(\mathcal{O}(h^2)\) central difference estimate \(D(h)\) for \(f'(x)\): \[D(h) = f'(x) + c_1 h^2 + c_2 h^4 + \dots\] Evaluating with a step size of \(h/2\): \[D(h/2) = f'(x) + c_1 (h/2)^2 + c_2 (h/2)^4 + \dots = f'(x) + \frac{c_1}{4} h^2 + \frac{c_2}{16} h^4 + \dots\] We can eliminate the \(h^2\) error term by computing: \[D_{\text{ext}} = \frac{4D(h/2) - D(h)}{3}\] This new estimate is \(\mathcal{O}(h^4)\) accurate.

Optimal Step Size: Truncation vs. Round-off Error

When computing numerically, decreasing \(h\) reduces the mathematical truncation error (the terms ignored in the Taylor series). However, as \(h \to 0\), \(f(x+h)\) and \(f(x)\) become nearly equal, leading to subtractive cancellation (a form of round-off error in floating-point arithmetic). Dividing by a very small \(h\) amplifies this round-off error. The total error is approximately \(E(h) = \frac{\epsilon}{h} + c h^n\), where \(\epsilon\) is machine precision. There exists an optimal step size \(h_{\text{opt}}\) that minimizes this total error.

3. Geometric/Visual Explanation

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return x**2 + 1

x0 = 1.0
h = 1.5

x_vals = np.linspace(-1, 3.5, 400)
y_vals = f(x_vals)

plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_vals, 'k-', label='f(x) = x^2 + 1')

# Points
plt.plot(x0, f(x0), 'ro')
plt.text(x0-0.2, f(x0)+0.5, '(x, f(x))')
plt.plot(x0+h, f(x0+h), 'bo')
plt.text(x0+h-0.2, f(x0+h)+0.5, '(x+h, f(x+h))')
plt.plot(x0-h, f(x0-h), 'go')
plt.text(x0-h-0.2, f(x0-h)+0.5, '(x-h, f(x-h))')

# True tangent
def tangent(x): return 2*x0*(x - x0) + f(x0)
plt.plot(x_vals, tangent(x_vals), 'r--', label='True Tangent')

# Forward diff secant
m_fwd = (f(x0+h) - f(x0)) / h
def sec_fwd(x): return m_fwd*(x - x0) + f(x0)
plt.plot(x_vals, sec_fwd(x_vals), 'b-', alpha=0.6, label='Forward Diff Secant')

# Backward diff secant
m_bwd = (f(x0) - f(x0-h)) / h
def sec_bwd(x): return m_bwd*(x - x0) + f(x0)
plt.plot(x_vals, sec_bwd(x_vals), 'g-', alpha=0.6, label='Backward Diff Secant')

# Central diff secant
m_cen = (f(x0+h) - f(x0-h)) / (2*h)
def sec_cen(x): return m_cen*(x - x0) + f(x0)
plt.plot(x_vals, sec_cen(x_vals), 'm-.', alpha=0.8, label='Central Diff Secant')

plt.ylim(0, 10)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Geometric Interpretation of Finite Differences')
plt.legend()
plt.grid(True)
plt.show()
Figure 1: Geometric representation of Finite Differences. 
import numpy as np
import matplotlib.pyplot as plt

def f2(x): return np.sin(x)
def df2(x): return np.cos(x)

x0 = 1.0
hs = np.logspace(-5, -1, 50)
err_fwd = []
err_cen = []

exact = df2(x0)
for h_val in hs:
    fwd = (f2(x0+h_val) - f2(x0)) / h_val
    cen = (f2(x0+h_val) - f2(x0-h_val)) / (2*h_val)
    err_fwd.append(abs(fwd - exact))
    err_cen.append(abs(cen - exact))

plt.figure(figsize=(8, 5))
plt.loglog(hs, err_fwd, label='Forward Diff Error (O(h))', marker='.')
plt.loglog(hs, err_cen, label='Central Diff Error (O(h^2))', marker='.')
plt.loglog(hs, hs, 'k--', label='Slope = 1')
plt.loglog(hs, hs**2, 'r--', label='Slope = 2')

plt.xlabel('Step size h')
plt.ylabel('Absolute Error')
plt.title('Error Scaling of Numerical Differentiation')
plt.legend()
plt.grid(True, which="both", ls="--")
plt.show()
Figure 2: Log-log plot of truncation errors for Forward vs Central Differences.

4. Python Implementation

import numpy as np

def forward_difference(f, x, h=1e-5):
    """Computes the forward difference approximation of f'(x)."""
    return (f(x + h) - f(x)) / h

def backward_difference(f, x, h=1e-5):
    """Computes the backward difference approximation of f'(x)."""
    return (f(x) - f(x - h)) / h

def central_difference(f, x, h=1e-5):
    """Computes the central difference approximation of f'(x)."""
    return (f(x + h) - f(x - h)) / (2 * h)

def second_derivative(f, x, h=1e-5):
    """Computes the central difference approximation of f''(x)."""
    return (f(x + h) - 2 * f(x) + f(x - h)) / (h**2)

def richardson_extrapolation(f, x, h=1e-1):
    """Improves the central difference estimate using Richardson Extrapolation."""
    D_h = central_difference(f, x, h)
    D_h_half = central_difference(f, x, h / 2)
    return (4 * D_h_half - D_h) / 3

5. Solved Examples

🟢 Easy: Finite Differences for \(f(x) = x^3\)

Problem: Compute the forward, backward, and central difference approximations of \(f'(1)\) for \(f(x) = x^3\) using \(h=0.1\). Compare with the exact value.

Solution: The exact derivative is \(f'(x) = 3x^2\), so \(f'(1) = 3(1)^2 = 3\). Given \(h=0.1\): \(f(1) = 1^3 = 1\) \(f(1.1) = 1.1^3 = 1.331\) \(f(0.9) = 0.9^3 = 0.729\)

  • Forward Difference: \(f'(1) \approx \frac{f(1.1) - f(1)}{0.1} = \frac{1.331 - 1}{0.1} = 3.31\) Error: \(|3 - 3.31| = 0.31\)

  • Backward Difference: \(f'(1) \approx \frac{f(1) - f(0.9)}{0.1} = \frac{1 - 0.729}{0.1} = 2.71\) Error: \(|3 - 2.71| = 0.29\)

  • Central Difference: \(f'(1) \approx \frac{f(1.1) - f(0.9)}{0.2} = \frac{1.331 - 0.729}{0.2} = 3.01\) Error: \(|3 - 3.01| = 0.01\) (Much more accurate!)

🟡 Medium: Richardson Extrapolation for \(\cos(x)\)

Problem: Use Richardson Extrapolation to improve a central difference estimate of the derivative of \(f(x) = \cos(x)\) at \(x = \pi/4\), starting with \(h = 0.2\).

Solution: Exact derivative: \(f'(x) = -\sin(x)\), so \(f'(\pi/4) = -\sin(\pi/4) \approx -0.70710678\).

First, calculate central difference with \(h = 0.2\): \(D(0.2) = \frac{\cos(\pi/4 + 0.2) - \cos(\pi/4 - 0.2)}{2(0.2)}\) \(D(0.2) \approx \frac{0.5516847 - 0.8335028}{0.4} \approx -0.704545\)

Next, calculate central difference with \(h = 0.1\): \(D(0.1) = \frac{\cos(\pi/4 + 0.1) - \cos(\pi/4 - 0.1)}{2(0.1)}\) \(D(0.1) \approx \frac{0.6333192 - 0.7741549}{0.2} \approx -0.7041785\)

Actually, \(D(0.1)\) is closer to exact, the extrapolation yields: \(D_{\text{ext}} = \frac{4D(0.1) - D(0.2)}{3} \approx \frac{4(-0.7041785) - (-0.704545)}{3} \approx -0.704056...\) \(D_{\text{ext}}\) will yield an \(\mathcal{O}(h^4)\) approximation, matching the true value to more decimal places.

🔴 Hard: The Optimal Step Size for \(e^x\)

Problem: Compute the central difference derivative of \(f(x) = e^x\) at \(x=1\) for \(h\) ranging from \(10^{-1}\) down to \(10^{-16}\). Plot the error and find the optimal step size \(h_{\text{opt}}\).

Solution:

import numpy as np
import matplotlib.pyplot as plt

def f3(x): return np.exp(x)

x0 = 1.0
exact = np.exp(x0)
hs = np.logspace(-16, -1, 100)
errors = []

for h_val in hs:
    approx = (f3(x0+h_val) - f3(x0-h_val)) / (2*h_val)
    errors.append(abs(approx - exact))

plt.figure(figsize=(8, 5))
plt.loglog(hs, errors, marker='.', color='purple')
plt.xlabel('Step size h')
plt.ylabel('Absolute Error')
plt.title('Round-off vs Truncation Error (U-Curve)')
plt.grid(True)
plt.axvline(x=1e-5, color='r', linestyle='--', label='Approx Optimal h')
plt.legend()
plt.show()
Figure 3: Trade-off between truncation error and round-off error.

The plot clearly demonstrates the “U-curve”. As \(h\) decreases from \(10^{-1}\), the error drops because the mathematical truncation error \(\mathcal{O}(h^2)\) decreases. However, around \(h = 10^{-5}\), the error hits a minimum. For \(h < 10^{-5}\), the error rapidly increases due to the round-off error dominating the subtractive cancellation in floating-point arithmetic.

6. Practice Problems

  1. Find the forward, backward, and central difference approximations for the derivative of \(f(x) = \ln(x)\) at \(x=2\) with \(h=0.1\).
  2. Derive the \(\mathcal{O}(h^2)\) forward difference formula for the first derivative using Taylor series expansions for \(f(x+h)\) and \(f(x+2h)\).
  3. The central difference formula for the second derivative is given in the text. Determine its leading error term order by substituting the Taylor series for \(f(x+h)\) and \(f(x-h)\).
  4. Apply Richardson Extrapolation twice (to eliminate \(h^2\) and \(h^4\) terms) starting with central differences for \(f(x) = \sin(x)\) at \(x=\pi/3\) using \(h=0.4, 0.2, 0.1\).
  5. A student attempts to calculate the derivative of a very noisy experimental dataset using the central difference method. What problems will they encounter, and why might they need to use a curve-fitting method (like least squares) before differentiating?

7. Summary

Key Takeaways
  • Forward / Backward Difference: First-order approximations \(\mathcal{O}(h)\). \[f'(x) \approx \frac{f(x+h) - f(x)}{h}\]
  • Central Difference: Second-order approximation \(\mathcal{O}(h^2)\). \[f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}\]
  • Second Derivative (Central): Second-order approximation. \[f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}\]
  • Richardson Extrapolation: A technique to boost the order of accuracy of an approximation method (e.g., from \(\mathcal{O}(h^2)\) to \(\mathcal{O}(h^4)\)).
  • Optimal Step Size (\(h\)): We must balance mathematical truncation error (requires small \(h\)) and computational round-off error (requires sufficiently large \(h\)). For central difference, optimal \(h \approx \sqrt[3]{\epsilon}\), where \(\epsilon\) is machine precision.

← Back to Course Index