Interpolation methods in Python

What you will learn

This page is written for anyone who wants to understand how to interpolate data points in Python, not just copy a formula.

Chebyshev nodes vs equispaced nodes

Node placement can matter as much as the interpolation formula. With high-degree polynomial interpolation, equally spaced points may create large oscillations near the interval edges. Chebyshev nodes cluster near the endpoints and often reduce that behavior.

x_nodes = np.linspace(-1, 1, n)

k = np.arange(1, n + 1)
x_nodes = np.cos((2*k - 1) * np.pi / (2*n))

Newton interpolation in Python

Newton interpolation writes the polynomial using coefficients from a divided-difference table. It is useful because adding one new data point does not require rebuilding the whole polynomial from scratch.

import numpy as np
from polynomial_interpolation import (
    evaluate_newton_polynomial,
    newton_coefficients,
)

x_nodes = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
y_nodes = 1 / (1 + 25 * x_nodes**2)

coefficients = newton_coefficients(x_nodes, y_nodes)
value = evaluate_newton_polynomial(coefficients, x_nodes, 0.25)

print(float(value))

Lagrange interpolation in Python

Lagrange interpolation builds the same polynomial as Newton when the same nodes are used. The formula is direct and easy to read, so it is common in numerical methods courses.

import numpy as np
from scipy.interpolate import lagrange

x_nodes = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
y_nodes = 1 / (1 + 25 * x_nodes**2)

polynomial = lagrange(x_nodes, y_nodes)
print(polynomial(0.25))

Inside the Python implementation

The repository keeps the Newton method as small, readable NumPy functions. Here are the core pieces, straight from polynomial_interpolation.py.

Building the divided-difference table

The table starts with the y values in the first column. Each next column divides the difference of the previous column by the distance between the matching x-nodes. The Newton coefficients are simply the top row of the finished table.

def divided_difference_table(x, y):
    n = len(x)
    table = np.zeros((n, n))
    table[:, 0] = y

    for order in range(1, n):
        for row in range(n - order):
            numerator = table[row + 1, order - 1] - table[row, order - 1]
            denominator = x[row + order] - x[row]
            table[row, order] = numerator / denominator

    return table

# Newton coefficients are the top row of the table:
coefficients = divided_difference_table(x, y)[0, :]

Evaluating with nested multiplication (Horner's scheme)

Instead of expanding the polynomial, the evaluator walks the coefficients backwards and multiplies as it goes. This is numerically stable, works on a single point or a whole NumPy array, and keeps the cost linear in the number of nodes.

def evaluate_newton_polynomial(coefficients, x_nodes, x_eval):
    result = np.full_like(x_eval, coefficients[-1], dtype=float)

    for i in range(len(coefficients) - 2, -1, -1):
        result = result * (x_eval - x_nodes[i]) + coefficients[i]

    return result

Comparing four methods on the same nodes

The benchmark builds the Newton polynomial from this repository and the Lagrange, barycentric and cubic-spline interpolants from SciPy, then measures the maximum and L2 error against the true function on a dense grid.

from scipy.interpolate import (
    BarycentricInterpolator,
    InterpolatedUnivariateSpline,
    lagrange,
)

# Newton (implemented in this repository)
coefficients = newton_coefficients(x_nodes, y_nodes)
y_newton = evaluate_newton_polynomial(coefficients, x_nodes, x_grid)

# Lagrange, barycentric and cubic spline (SciPy)
y_lagrange = lagrange(x_nodes, y_nodes)(x_grid)
y_bary = BarycentricInterpolator(x_nodes, y_nodes)(x_grid)
y_spline = InterpolatedUnivariateSpline(x_nodes, y_nodes, k=3)(x_grid)

# Error against the true function on a dense grid
max_error = np.max(np.abs(y_newton - y_true))
l2_error = np.sqrt(np.trapezoid((y_newton - y_true)**2, x_grid))

Interactive interpolation demo

Change the function, nodes, and evaluation point. The plot shows why Newton and Lagrange give the same polynomial for the same nodes, and why Chebyshev nodes often reduce oscillation.

Newton and Lagrange interpolation in Python

python -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt
python polynomial_interpolation.py