Using Chebyshev points for Lagrange interpolation


Lagrange interpolation is a useful polynomial interpolation formula that lets us fit to an arbitrary function. While we might expect that increasing the number of points would decrease the average error of our polynomial approximation, this intuition fails because of Runge's phenomenon. If we pick our points carefully, though, we can overcome this phenomenon and still achieve a great approximation.

Specifically, in this post, we're going to approximate the function 1 / (1 + x^6) over [0, 5] by using two different sets of points: a set of equidistant points across the interval vs. Chebyshev nodes. We'll compare the average error (MSE) and the maximum error over the interval.

Helper functions

We'll add a function that returns the Lagrange interpolating polynomial for a set of points:

from operator import mul
from functools import reduce

def lagrange(pts, f):
    """ Return a Lagrange interpolating polynomial for f
    at the points pts. """

    def basis(others):
        # w_{i}(x) = (x - x_{1}) * ... * (x - x_{i-1}) * (x - x_{i+1}) * ...
        return lambda x : reduce(mul, [x - o for o in others])

    # Get each basis polynomial and coefficients
    bases = [basis(pts[:i] + pts[i + 1:]) for i in range(len(pts))]
    # a_{i} = f(x_{i}) / w_{i}(x_{i})
    weights = [f(pt) / bases[i](pt) for i, pt in enumerate(pts)]

    return lambda x : sum([weights[i] * bases[i](x) for i in range(len(pts))])

This is a basic implementation of Lagrange interpolation - nothing special here.

Next, we'll add a couple of helper functions to compute the average and maximum error:

def average_error(p, f):
    # Get the MSE of p over 1000 points on the interval [0, 5]
    N = 1000
    error = 0
    for pt in [5 * i / N for i in range(N)]:
        error += (p(pt) - f(pt))**2
    return error / N

def max_error(p, f):
    error = float("-inf")
    for pt in [5 * i / 1000 for i in range(1000)]:
        curr_error = abs(p(pt) - f(pt))
        if curr_error > error:
            error = curr_error
    return curr_error

Computing for equidistant vs. Chebyshev points

Finally, we'll add our function to approximate:

def g(x):
    return 1 / (1 + x**6)

And some functions to compute and evaluate the interpolating polynomial for N points, whether equidistant or Chebyshev:

def evaluate(pts, f, tag):
    print("Interpolating for {0}".format(tag))
    poly = lagrange(pts, f)
    print("Average error: {0}".format(average_error(poly, f)))
    print("Max error: {0}".format(max_error(poly, f)))

def equidistant(n):
    # Interpolate for n points from [0, 5]
    pts = [5 * i / n for i in range(n)]
    evaluate(pts, g, "{0} equidistant points".format(n))

import math
def chebyshev(n):
    # Interpolate for n points from [0, 5]
    pts = [(5 + 5 * math.cos((2 * i + 1) * math.pi / (2 * n + 2))) / 2 for i in range(n)]
    evaluate(pts, g, "{0} Chebyshev points".format(n))

n_values = [3, 5, 10, 20, 30]
for n in n_values:


Here's what the results show:

Mean Squared Error:

Number of pointsEquidistantChebyshev

Maximum Error:

Number of pointsEquidistantChebyshev

I initially wanted to graph these results, but as you can see, the error for the Chebyshev nodes is so low that it would hardly be visible.

All in all, clearly it's preferable to use Chebyshev nodes. This is especially clear as the number of points increases. I hope this was helpful!

Tags: chebyshev lagrange interpolation python

See also: Adding GitHub to Your Hosts File

Back to all posts

Neel Somani

About the Author

Neel Somani is a software engineer in the Bay Area. In addition to computer science, he's interested in philosophy and entrepreneurship. You can follow him on LinkedIn and Twitter.

comments powered by Disqus