# 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 $$\frac{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 points Equidistant Chebyshev 3 0.083449016830 0.123504372195 5 0.082739268475 0.08960932097 10 0.64207073475 0.0127952042780 20 26.4328328320 0.0005982536404 30 435.529425177 9.3148785757e-06

Maximum Error:

 Number of points Equidistant Chebyshev 3 0.86442582784 0.0030308889852 5 1.34512664557 0.0063971895249 10 6.2268209092 0.000121486335635 20 61.205287745 2.1722048813e-05 30 312.690373601 4.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 