Part V: Interpolation & IntegrationChapter 5 of 6

Interpolation & Numerical Integration

Reconstructing functions from discrete data, computing definite integrals to high precision, and the Fast Fourier Transform: the most important algorithm of the 20th century.

1. Lagrange Interpolation

Given \(n+1\) data points \((x_0, y_0), \ldots, (x_n, y_n)\), the unique polynomial of degree \(\leq n\) passing through all points is:

\(P(x) = \sum_{i=0}^{n} y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}\)

Each Lagrange basis polynomial \(L_i(x)\) equals 1 at \(x_i\) and 0 at all other nodes.

Runge's Phenomenon

High-degree polynomial interpolation on equally-spaced points can oscillate wildly near the edges. For \(f(x) = 1/(1+25x^2)\) on \([-1,1]\), the interpolation error grows as the degree increases! Fix: use Chebyshev nodes\(x_k = \cos(k\pi/n)\) which cluster near the endpoints.

Lagrange Interpolation and Runge's Phenomenon

Python

script.py40 lines

import numpy as np

def lagrange_interp(x_nodes, y_nodes, x_eval):
    """Evaluate Lagrange interpolating polynomial."""
    n = len(x_nodes)
    result = np.zeros_like(x_eval, dtype=float)
    for i in range(n):
        Li = np.ones_like(x_eval, dtype=float)
        for j in range(n):
            if j != i:
                Li *= (x_eval - x_nodes[j]) / (x_nodes[i] - x_nodes[j])
        result += y_nodes[i] * Li
    return result

# Runge's function
f = lambda x: 1.0 / (1 + 25 * x**2)
x_fine = np.linspace(-1, 1, 200)

print("Runge's Phenomenon: f(x) = 1/(1+25x^2)")
print(f"{'n':>4}  {'Max err (uniform)':>18}  {'Max err (Chebyshev)':>20}")
print("-" * 46)

for n in [5, 8, 10, 15, 20]:
    # Uniform nodes
    x_unif = np.linspace(-1, 1, n + 1)
    y_unif = f(x_unif)
    p_unif = lagrange_interp(x_unif, y_unif, x_fine)
    err_unif = np.max(np.abs(p_unif - f(x_fine)))

# Chebyshev nodes
    x_cheb = np.cos(np.pi * np.arange(n + 1) / n)
    y_cheb = f(x_cheb)
    p_cheb = lagrange_interp(x_cheb, y_cheb, x_fine)
    err_cheb = np.max(np.abs(p_cheb - f(x_fine)))

print(f"{n:4d}  {err_unif:18.6e}  {err_cheb:20.6e}")

print("\nUniform nodes: error GROWS with degree (Runge's phenomenon)")
print("Chebyshev nodes: error DECREASES with degree")
print("Lesson: node placement matters more than polynomial degree!")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

2. Newton's Divided Differences

A more computationally efficient form using divided differences. Unlike Lagrange, adding a new data point requires only one new coefficient (no recomputation):

\(P(x) = f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)(x-x_1) + \cdots\)

Newton's Divided Differences

Python

script.py38 lines

import numpy as np

def divided_diff(x, y):
    """Compute Newton's divided difference coefficients."""
    n = len(x)
    coeff = y.copy().astype(float)
    for j in range(1, n):
        for i in range(n - 1, j - 1, -1):
            coeff[i] = (coeff[i] - coeff[i - 1]) / (x[i] - x[i - j])
    return coeff

def newton_eval(coeff, x_nodes, x_eval):
    """Evaluate Newton interpolating polynomial."""
    n = len(coeff)
    result = coeff[-1] * np.ones_like(x_eval, dtype=float)
    for i in range(n - 2, -1, -1):
        result = result * (x_eval - x_nodes[i]) + coeff[i]
    return result

# Interpolate sin(x) with increasing number of points
f = lambda x: np.sin(x)
x_test = np.linspace(0, np.pi, 100)

print("Newton Interpolation of sin(x) on [0, pi]")
print(f"{'n':>4}  {'Max error':>14}  {'Coefficients (first 5)':>30}")
print("-" * 52)

for n in [3, 5, 7, 9, 11]:
    x_nodes = np.linspace(0, np.pi, n)
    y_nodes = f(x_nodes)
    coeff = divided_diff(x_nodes, y_nodes)
    p_eval = newton_eval(coeff, x_nodes, x_test)
    err = np.max(np.abs(p_eval - f(x_test)))
    coeff_str = ", ".join(f"{c:.4f}" for c in coeff[:5])
    print(f"{n:4d}  {err:14.6e}  {coeff_str}")

print("\nNewton form is ideal for incremental data:")
print("adding a point only needs one new coefficient.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

3. Cubic Splines

Instead of one high-degree polynomial, use piecewise cubics joined smoothly (\(C^2\)continuity). No Runge phenomenon, local control, and the interpolant minimizes the total bending energy \(\int |u''(x)|^2 dx\).

Cubic Spline Interpolation

Python

script.py53 lines

import numpy as np
from scipy.interpolate import CubicSpline

# Compare polynomial vs spline for Runge's function
f = lambda x: 1.0 / (1 + 25 * x**2)
x_fine = np.linspace(-1, 1, 500)

print("Cubic Spline vs Polynomial: Runge's function")
print(f"{'n':>4}  {'Poly error':>14}  {'Spline error':>14}  {'Winner':>8}")
print("-" * 46)

for n in [5, 10, 15, 20, 30, 50]:
    x_nodes = np.linspace(-1, 1, n + 1)
    y_nodes = f(x_nodes)

# Cubic spline
    cs = CubicSpline(x_nodes, y_nodes)
    spline_eval = cs(x_fine)
    err_spline = np.max(np.abs(spline_eval - f(x_fine)))

# Polynomial (using numpy)
    try:
        p = np.polyfit(x_nodes, y_nodes, n)
        poly_eval = np.polyval(p, x_fine)
        err_poly = np.max(np.abs(poly_eval - f(x_fine)))
        if err_poly > 1e6:
            err_poly_str = "DIVERGED"
            winner = "Spline"
        else:
            err_poly_str = f"{err_poly:.6e}"
            winner = "Spline" if err_spline < err_poly else "Poly"
    except:
        err_poly_str = "FAILED"
        winner = "Spline"

print(f"{n:4d}  {err_poly_str:>14}  {err_spline:14.6e}  {winner:>8}")

print("\nSplines always converge; polynomials can diverge!")

# Show spline smoothness
print("\nSpline derivative continuity check (n=10):")
x_nodes = np.linspace(-1, 1, 11)
cs = CubicSpline(x_nodes, f(x_nodes))
for i in range(1, len(x_nodes) - 1):
    xi = x_nodes[i]
    eps = 1e-8
    d1_left = cs(xi - eps, 1)  # first derivative from left
    d1_right = cs(xi + eps, 1)
    d2_left = cs(xi - eps, 2)  # second derivative from left
    d2_right = cs(xi + eps, 2)
    if i <= 3:
        print(f"  Node x={xi:.1f}: |d1_jump| = {abs(d1_left-d1_right):.2e}, |d2_jump| = {abs(d2_left-d2_right):.2e}")
print("  ... (C^2 continuous at all interior knots)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

4. Gaussian Quadrature

The optimal choice of nodes and weights for numerical integration. An \(n\)-point Gauss-Legendre rule integrates polynomials of degree \(\leq 2n-1\) exactly!

\(\int_{-1}^{1} f(x) \, dx \approx \sum_{i=1}^{n} w_i f(x_i)\)

Nodes \(x_i\) are roots of Legendre polynomials. Error: \(O(h^{2n})\) for \(n\)-point rule.

Gaussian Quadrature: Optimal Integration

Python

script.py44 lines

import numpy as np
from scipy.integrate import quad, fixed_quad

# Compare quadrature rules
def trapezoid(f, a, b, n):
    x = np.linspace(a, b, n + 1)
    h = (b - a) / n
    return h * (0.5*f(x[0]) + np.sum(f(x[1:-1])) + 0.5*f(x[-1]))

def simpson(f, a, b, n):
    if n % 2 != 0: n += 1
    x = np.linspace(a, b, n + 1)
    h = (b - a) / n
    return h/3 * (f(x[0]) + 4*np.sum(f(x[1::2])) + 2*np.sum(f(x[2:-1:2])) + f(x[-1]))

# Test: integral of exp(-x^2) from 0 to 1
f = lambda x: np.exp(-x**2)
exact, _ = quad(f, 0, 1)

print("Integration of exp(-x^2) from 0 to 1")
print(f"Exact (scipy.quad): {exact:.15f}\n")

print(f"{'Method':>20}  {'n points':>10}  {'Result':>16}  {'|Error|':>12}")
print("-" * 62)

for n in [4, 8, 16, 32]:
    # Trapezoid: O(h^2)
    I_trap = trapezoid(f, 0, 1, n)
    print(f"{'Trapezoid':>20}  {n:10d}  {I_trap:16.12f}  {abs(I_trap-exact):12.2e}")

print()
for n in [4, 8, 16, 32]:
    # Simpson: O(h^4)
    I_simp = simpson(f, 0, 1, n)
    print(f"{'Simpson':>20}  {n:10d}  {I_simp:16.12f}  {abs(I_simp-exact):12.2e}")

print()
for n in [2, 4, 6, 8]:
    # Gauss-Legendre: O(h^{2n})
    I_gauss, _ = fixed_quad(f, 0, 1, n=n)
    print(f"{'Gauss-Legendre':>20}  {n:10d}  {I_gauss:16.12f}  {abs(I_gauss-exact):12.2e}")

print("\nGauss-Legendre with 8 points beats Simpson with 32 points!")
print("Optimal node placement is the key insight.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. Monte Carlo Integration

For high-dimensional integrals, deterministic methods suffer the curse of dimensionality(cost grows as \(N^d\)). Monte Carlo converges as \(O(1/\sqrt{N})\) regardless of dimension.

\(\int_\Omega f(x) \, dx \approx \frac{V(\Omega)}{N} \sum_{i=1}^{N} f(x_i), \quad x_i \sim \text{Uniform}(\Omega)\)

Error: \(O(\sigma / \sqrt{N})\) by the Central Limit Theorem. Dimension-independent!

Monte Carlo Integration: High Dimensions

Python

script.py43 lines

import numpy as np

def monte_carlo_integrate(f, bounds, N):
    """Monte Carlo integration over a hyperrectangle."""
    d = len(bounds)
    # Sample uniformly
    rng = np.random.RandomState(42)
    samples = np.zeros((N, d))
    volume = 1.0
    for i, (lo, hi) in enumerate(bounds):
        samples[:, i] = rng.uniform(lo, hi, N)
        volume *= (hi - lo)

values = np.array([f(x) for x in samples])
    estimate = volume * np.mean(values)
    std_err = volume * np.std(values) / np.sqrt(N)
    return estimate, std_err

# 1D test: integral of sin(x) from 0 to pi = 2
print("=== 1D: integral of sin(x) from 0 to pi ===")
f1d = lambda x: np.sin(x[0])
for N in [100, 1000, 10000, 100000]:
    I, err = monte_carlo_integrate(f1d, [(0, np.pi)], N)
    print(f"  N = {N:>7}: I = {I:.6f} +/- {err:.6f} (exact = 2.0, error = {abs(I-2):.4f})")

# High-dimensional: volume of unit ball in d dimensions
print("\n=== Volume of unit ball in d dimensions ===")
print(f"{'d':>4}  {'MC estimate':>12}  {'Exact':>12}  {'Rel error':>10}")
print("-" * 42)
for d in [2, 4, 6, 8, 10]:
    def indicator(x):
        return 1.0 if np.sum(x**2) <= 1.0 else 0.0

bounds = [(-1, 1)] * d
    I, err = monte_carlo_integrate(indicator, bounds, N=100000)
    # Exact: pi^(d/2) / Gamma(d/2 + 1)
    from math import gamma, pi
    exact = pi**(d/2) / gamma(d/2 + 1)
    print(f"{d:4d}  {I:12.4f}  {exact:12.4f}  {abs(I-exact)/exact:10.4f}")

print("\nMonte Carlo: error ~ 1/sqrt(N) regardless of dimension")
print("Deterministic (trapezoid) in 10D with 10 points/dim: 10^10 evaluations!")
print("Monte Carlo: 100,000 evaluations gives reasonable accuracy.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. The Fast Fourier Transform (FFT)

The Discrete Fourier Transform converts a signal from time/space domain to frequency domain. The naive DFT costs \(O(N^2)\); the FFT (Cooley-Tukey, 1965) achieves \(O(N \log N)\), making spectral analysis practical for large datasets.

\(\hat{f}_k = \sum_{n=0}^{N-1} f_n \, e^{-2\pi i k n / N}\)

DFT: \(O(N^2)\) multiplications. FFT: \(O(N \log N)\) via divide-and-conquer. For \(N = 10^6\): 10¹² vs \(2 \times 10^7\) -- a 50,000x speedup!

FFT: Signal Analysis and Filtering

Python

script.py48 lines

import numpy as np

# Create a signal: sum of sinusoids + noise
N = 1024
dt = 0.001
t = np.arange(N) * dt

# Signal: 50 Hz + 120 Hz + noise
f_signal = 0.7*np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)
noise = 0.5 * np.random.RandomState(42).randn(N)
signal = f_signal + noise

# FFT
spectrum = np.fft.fft(signal)
freqs = np.fft.fftfreq(N, dt)

# Power spectrum (positive frequencies only)
pos_mask = freqs > 0
power = 2.0/N * np.abs(spectrum[pos_mask])
pos_freqs = freqs[pos_mask]

# Find dominant frequencies
peak_indices = np.argsort(power)[-5:][::-1]

print("FFT Signal Analysis")
print(f"Signal: 0.7*sin(2pi*50*t) + sin(2pi*120*t) + noise")
print(f"Samples: N = {N}, dt = {dt}s, T = {N*dt}s")
print(f"Frequency resolution: {1/(N*dt):.1f} Hz")
print(f"Nyquist frequency: {1/(2*dt):.0f} Hz\n")

print("Top 5 frequency components:")
print(f"{'Freq (Hz)':>10}  {'Amplitude':>10}  {'Expected':>10}")
print("-" * 34)
for idx in peak_indices:
    expected = ""
    if abs(pos_freqs[idx] - 50) < 2:
        expected = "~0.70"
    elif abs(pos_freqs[idx] - 120) < 2:
        expected = "~1.00"
    print(f"{pos_freqs[idx]:10.1f}  {power[idx]:10.4f}  {expected:>10}")

# FFT vs DFT timing comparison
print(f"\n--- Complexity Comparison ---")
for n in [2**8, 2**10, 2**12, 2**14, 2**16]:
    dft_ops = n * n
    fft_ops = n * np.log2(n)
    speedup = dft_ops / fft_ops
    print(f"N = {n:>6}: DFT = {dft_ops:>12.0f}, FFT = {fft_ops:>10.0f}, speedup = {speedup:>7.0f}x")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Method Selection Guide

Task	Best Method	When to Use
Smooth interpolation	Cubic splines	Default choice, avoids Runge phenomenon
Low-D integration	Gaussian quadrature	Smooth integrands, d = 1-3
High-D integration	Monte Carlo	d > 4, irregular domains
Spectral analysis	FFT	\(O(N \log N)\), periodic data
Adaptive integration	scipy.integrate.quad	Black-box, handles singularities