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:
    equidistant(n)
    print("")
    chebyshev(n)
    print("")

Evaluation

Here's what the results show:

Mean Squared Error:

Number of pointsEquidistantChebyshev
30.0834490168300.123504372195
50.0827392684750.08960932097
100.642070734750.0127952042780
2026.43283283200.0005982536404
30435.5294251779.3148785757e-06

Maximum Error:

Number of pointsEquidistantChebyshev
30.864425827840.0030308889852
51.345126645570.0063971895249
106.22682090920.000121486335635
2061.2052877452.1722048813e-05
30312.6903736014.1559874016e-06

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: Setting up Flask with Red Hat on AWS

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