Unit 3: Numerical Differentiation, Richardson Extrapolation & Advanced Quadrature

1. Finite Difference Approximations of Derivatives

When an analytical expression for $f(x)$ is unavailable or excessively complicated, derivatives must be computed from discrete evaluations. By truncating Taylor series expansions, we construct finite difference stencils.

Let $h > 0$ denote the discretization grid spacing.

1. Forward and Backward Difference Approximations ($\mathcal{O}(h)$):

Expanding $f(x + h)$ in Taylor series:

Solving for $f'(x)$:

Similarly, expanding $f(x - h)$:

2. Three-Point Centered Difference Formula ($\mathcal{O}(h^2)$):

Subtracting the backward Taylor expansion from the forward Taylor expansion:

Solving for $f'(x)$ yields the Centered Difference Formula:

The linear error term cancels identically due to symmetry, elevating the accuracy from $\mathcal{O}(h)$ to $\mathcal{O}(h^2)$.

3. Second Derivative Centered Difference Stencil ($\mathcal{O}(h^2)$):

Adding the two Taylor expansions:

2. The Fundamental Truncation-Roundoff Dilemma

In pure analysis, the derivative is obtained by taking $\lim_{h \to 0}$. In digital computing, as $h \to 0$, finite precision arithmetic destroys the solution!

Suppose each function evaluation is contaminated by roundoff error bounded by $\epsilon = \epsilon_{\text{mach}} |f(x)|$:

The computationally computed centered difference derivative is:

The total computational error $E_{\text{total}}(h) = |f'(x) - \tilde{D}_h^0 f(x)|$ satisfies:

where $M = \max |f'''(\xi)|$.

• As $h$ decreases, truncation error decays as $h^2$.
• But as $h$ decreases, round-off error explodes as $\frac{1}{h}$ due to division by near-zero!

Optimal Step Size ($h^*$):

Minimizing $E(h) = \frac{M h^2}{6} + \frac{\epsilon}{h}$ by setting $\frac{dE}{dh} = \frac{M h}{3} - \frac{\epsilon}{h^2} = 0$:

For double-precision floating point ($\epsilon \approx 10^{-16}$), the optimal step size for centered differences is $h^ \approx 10^{-16/3} \approx 10^{-5.3} \approx 10^{-5}$, and the achievable accuracy is bounded by $E(h^) \approx \mathcal{O}(\epsilon^{2/3}) \approx 10^{-11}$. One cannot achieve full 16-digit machine precision with simple finite differences!

1. The Principle of Richardson Error Cancellation

Richardson Extrapolation is a general numerical acceleration technique that systematically eliminates the dominant term in a truncation error expansion.

Suppose an approximation $N_1(h)$ to an unknown exact quantity $M$ has an asymptotic expansion with even powers of $h$:

where the coefficients $K_j$ are independent of $h$. Now evaluate the same approximation with half the step size, $h/2$:

Multiply the second equation by 4 and subtract the first equation:

Dividing by 3:

Define the new approximation $N_2(h) = N_1(h/2) + \frac{N_1(h/2) - N_1(h)}{3}$. Its error is $\mathcal{O}(h^4)$, having completely eliminated the $\mathcal{O}(h^2)$ term without evaluating higher derivatives!

2. General Extrapolation Tableau

Repeating this process recursively generates a triangular tableau:

• Column 1: $\mathcal{O}(h^2)$ approximations ($N_1(h), N_1(h/2), N_1(h/4), \dots$).
• Column 2: $\mathcal{O}(h^4)$ approximations.
• Column 3: $\mathcal{O}(h^6)$ approximations.
• Column $m$: $\mathcal{O}(h^{2m})$ approximations.

1. General Framework of Interpolatory Quadrature

Numerical integration (quadrature) replaces a definite integral $I = \int_a^b f(x) dx$ with a finite weighted sum:

In Newton-Cotes formulas, the nodes $x_i$ are equispaced: $x_i = a + i h$ with $h = \frac{b - a}{n}$. The weights $w_i$ are obtained by integrating the cardinal Lagrange basis polynomials:

2. The Trapezoidal Rule ($n = 1$)

Connecting $(a, f(a))$ and $(b, f(b))$ with a linear interpolant:

Composite Trapezoidal Rule:

Dividing $[a, b]$ into $m$ subintervals of length $h = \frac{b - a}{m}$:

Theorem 3.1 (Trapezoidal Truncation Error):

If $f \in C^2[a, b]$, the global truncation error is:

3. Simpson's 1/3 Rule ($n = 2$)

Fitting a parabola through $(x_0, f_0), (x_1, f_1), (x_2, f_2)$ with $h = \frac{b - a}{2}$:

Composite Simpson's 1/3 Rule ($m$ even):

Dividing $[a, b]$ into an even number $m$ of subintervals:

Theorem 3.2 (Simpson's Error Formula & Precision Miracle):

If $f \in C^4[a, b]$, the global error is:

Significance: Because Simpson's rule uses a symmetric quadratic interpolant, the cubic error term $\int_{-h}^h x^3 dx = 0$ vanishes by symmetry. Consequently, Simpson's rule integrates polynomials of degree up to 3 exactly without any error!

4. Simpson's 3/8 Rule ($n = 3$) and Weddle's Rule ($n = 6$)

• Simpson's 3/8 Rule: Fits a cubic through 4 equispaced points:

Composite error over $[a, b]$: $E_{3/8} = -\frac{(b - a) h^4}{80} f^{(4)}(\xi)$.

• Weddle's Rule ($n = 6$): Six-panel formula with optimal integer weights:

Global error: $E_W = -\frac{(b - a) h^6}{1400} f^{(6)}(\xi) = \mathcal{O}(h^6)$.

1. Romberg Quadrature

Romberg integration applies Richardson extrapolation systematically to the composite Trapezoidal rule:

The column entries $R_{k, j}$ have error order $\mathcal{O}(h^{2j})$. The diagonal entries $R_{k, k}$ converge with extraordinary rapidity.

2. Gauss-Legendre Quadrature Theory

Newton-Cotes formulas fix nodes $x_i$ to be equispaced and determine $n + 1$ weights, achieving degree of precision $n$ (or $n+1$ for even $n$). Gaussian Quadrature frees both the $n$ nodes $x_1, \dots, x_n$ and the $n$ weights $w_1, \dots, w_n$ (a total of $2n$ degrees of freedom), achieving maximal algebraic degree of precision.

Theorem 3.3 (Maximal Precision of Gaussian Quadrature):

Let $\{P_n(x)\}$ be the family of orthogonal Legendre polynomials on $[-1, 1]$ satisfying:

Let $x_1, x_2, \dots, x_n$ be the $n$ distinct roots of $P_n(x) = 0$ in $(-1, 1)$, and let:

Then the quadrature formula:

is strictly exact for all polynomials of degree up to $2n - 1$.

Proof: Let $p(x) \in \mathbb{P}_{2n-1}$. By polynomial division, divide $p(x)$ by $P_n(x)$:

where $\deg(q) \le n - 1$ and $\deg(r) \le n - 1$. Integrate both sides over $[-1, 1]$:

Since $\deg(q) \le n - 1$ and $P_n$ is orthogonal to all polynomials of degree $< n$, the first integral vanishes identically: $\int_{-1}^1 q(x) P_n(x) dx = 0$. Thus: $\int_{-1}^1 p(x) dx = \int_{-1}^1 r(x) dx$. Now evaluate the quadrature rule on $p(x)$:

Because the nodes $x_i$ are roots of $P_n(x)$, $P_n(x_i) = 0$. Thus:

Since $\deg(r) \le n - 1$, the interpolatory quadrature rule with $n$ nodes integrates $r(x)$ exactly:

Hence $\int_{-1}^1 p(x) dx = \sum_{i=1}^n w_i p(x_i)$, establishing precision $2n - 1$. $\blacksquare$

Transformation to Arbitrary Interval $[a, b]$:

To evaluate $\int_a^b f(t) dt$, map $t \in [a, b]$ to $x \in [-1, 1]$ via substitution: