Interpolation and curve fitting

Engineers constantly work with functions known only at discrete points: properties in a steam table, a pump's head at tested flow rates, survey elevations along a road, sensor readings every second. Two tasks follow. Interpolation finds values between exact data points, with a curve that passes through every one of them. Curve fitting finds a simple formula that follows noisy data closely without passing through each point. Choosing the right one, and the right curve, is the subject of this chapter.

Linear interpolation

The simplest interpolant joins neighbouring points with straight lines. Between \((x_0, y_0)\) and \((x_1, y_1)\):

Every engineering student does this with property tables. Its error is about \(\frac{1}{8}h^2|f''|\) for spacing h, small when the function is nearly straight over the interval and larger where it curves.

Polynomial interpolation

Through n + 1 points with distinct x values there is exactly one polynomial of degree at most n. Two convenient ways to write it:

• Lagrange form: \(p(x) = \sum_i y_i L_i(x)\), where each basis polynomial \(L_i(x) = \prod_{j \ne i}\frac{x - x_j}{x_i - x_j}\) equals 1 at \(x_i\) and 0 at the other nodes. Elegant for theory.
• Newton form with divided differences: \(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\). Adding a point adds one term, without recomputing the others. Better for hand and computer work.

Divided differences are computed in a table: \(f[x_i, x_{i+1}] = \dfrac{f[x_{i+1}] - f[x_i]}{x_{i+1} - x_i}\), and each higher difference is the difference of the two below it divided by the span of x values involved.

The danger of high degree: Runge's phenomenon

It is tempting to think more points and a higher degree always help. They do not. The error of polynomial interpolation is

and with equally spaced points the product term becomes huge near the ends of the interval. For the smooth function \(f(x) = 1/(1 + 25x^2)\) on [−1, 1], a degree-10 polynomial through 11 equally spaced points has a maximum error of 1.92, far larger than the function itself, and the error grows as the degree increases.

Two cures: cluster the points toward the ends (Chebyshev nodes \(x_k = \cos\frac{(2k + 1)\pi}{2(n + 1)}\) reduce the degree-10 error to 0.11), or, far more common in engineering, use piecewise polynomials of low degree.

Cubic splines

A cubic spline uses a separate cubic between each pair of adjacent points, chosen so that the curve and its first and second derivatives are continuous at every interior point. The result is smooth to the eye and to the touch, like the flexible strips of wood (splines) that shipwrights and draftsmen once bent through pegs to draw hull lines.

The continuity conditions give a tridiagonal linear system for the second derivatives at the nodes (Lesson 3), solved in O(n) operations. Two end conditions complete it: natural (zero second derivative at the ends), clamped (specified end slopes) or not-a-knot (the default in many libraries).

Through the same 11 points as above, a cubic spline has maximum error 0.022 for the Runge function, almost 100 times better than the degree-10 polynomial. Splines are used for road and railway alignments, cam profiles, CAD curves and surfaces, robot trajectories and smooth interpolation of tabulated data.

Least-squares curve fitting

Measured data contain noise, so a curve through every point would chase the noise. Instead, choose a simple model with a few parameters and fit it by least squares: minimize the sum of squared residuals \(\sum (y_i - \hat{y}_i)^2\), as in STAT-235's straight-line regression.

For a polynomial model \(y = c_0 + c_1x + \cdots + c_mx^m\) with m much smaller than the number of points, the fit is a linear least-squares problem \(V\mathbf{c} \approx \mathbf{y}\), where V is the matrix whose columns are \(1, x, x^2, \ldots\). It can be solved through the normal equations \(V^TV\mathbf{c} = V^T\mathbf{y}\), but \(V^TV\) squares the condition number of V (Lesson 3), so libraries use the more stable QR factorization instead (np.polyfit, np.linalg.lstsq, MATLAB polyfit).

Overfitting

A degree-5 polynomial passes exactly through all six pump points (R² = 1). Inside the data range it differs from the quadratic by only a few centimetres, but beyond it the two separate fast: at 65 L/s the quadratic predicts 9.6 m of head, while the degree-5 curve predicts −1.6 m, which is physically impossible. Extra parameters fit the noise, not the physics. Choose the simplest model that leaves residuals without a pattern, prefer models with a physical basis, and never trust extrapolation of a high-degree fit.

Nonlinear models

Many physical laws are not polynomials: exponential decay, power laws, saturation curves. Two approaches:

• Linearize when possible: \(y = ae^{bx}\) becomes \(\ln y = \ln a + bx\); \(y = ax^b\) becomes \(\ln y = \ln a + b\ln x\). Fit a line to the transformed data. Quick, but it weights the errors differently from the original scale.
• Nonlinear least squares: minimize the squared residuals in the original form directly, iteratively (Gauss–Newton or Levenberg–Marquardt, scipy.optimize.curve_fit). It needs reasonable starting values and returns parameter uncertainties too.

import numpy as np
from scipy.interpolate import CubicSpline
from scipy.optimize import curve_fit

T = [100, 110, 120]; P = [101.418, 143.376, 198.665]      # steam table
print(f"linear at 105 C: {np.interp(105, T, P):.2f} kPa")

Q = np.array([0, 10, 20, 30, 40, 50.0])
H = np.array([42.0, 41.2, 38.9, 34.8, 29.6, 22.7])
c = np.polyfit(Q, H, 2)
print("quadratic coefficients:", ", ".join(f"{v:.5f}" for v in c), f" H(35) = {np.polyval(c, 35):.1f} m")

x = np.linspace(-1, 1, 11)
spline = CubicSpline(x, 1 / (1 + 25*x**2))
print(f"spline at 0.95: {spline(0.95):.4f}  (true {1/(1 + 25*0.95**2):.4f})")

t = np.arange(1, 8.0)
y = np.array([48.8, 91.5, 131.3, 153.0, 164.3, 187.1, 197.5])
(L0, k), cov = curve_fit(lambda t, L, k: L*(1 - np.exp(-k*t)), t, y, p0=[200, 0.1])
print(f"L0 = {L0:.0f} mg/L, k = {k:.3f} /day")

Output: linear at 105 C: 122.40 kPa; quadratic coefficients −0.00752, −0.01096, 42.03214, H(35) = 32.4 m; spline at 0.95: 0.0436 (true 0.0424); L0 = 240 mg/L, k = 0.247 /day.