Interpolation is a common technique used in numerical analysis and engineering to estimate values between known data points. It is a powerful tool that can be used to create smooth curves and surfaces based on a limited amount of data. However, traditional interpolation methods, such as polynomial interpolation, can suffer from high error rates when working with large datasets or when trying to approximate highly complex functions. Chebyshev grids offer a solution to this problem by providing a more accurate and efficient way of interpolating data. In this post, we will explain the problem with equidistant grids, the idea behind Chebyshev grids and show you how to implement them in Python. We will also demonstrate how to use Chebyshev grids to interpolate a sample dataset and compare the results to traditional interpolation methods. By the end of this post, you will have a solid understanding of how to use Chebyshev grids for interpolation with extreme accuracy in Python.

The problem with equidistant grids

Suppose you have a function given on an equidistant grid and you want to find the interpolating function using polynomials. This problem is easy to solve using the Lagrange method (see this article. ).

But there is a problem. Indeed, a serious problem. Let’s look at an example. Let’s image we are sampling the function

on the domain x∈[−1,1] using an equidistant grid of N=15 points:

import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np


N = 15

def f(x):
    return np.exp(-(5*x)**2) + x

x = np.linspace(-1, 1, N)
xx = np.linspace(-1, 1, 500)
f_eq = f(x) 

plt.plot(x, f_eq, 'o')
plt.plot(xx, f(xx), label='f(x)')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')
plt.legend();

Now let’s find an interpolating function for the sampled values. We will use the Lagrange method described in this article . For convenience, here is a function that gives you the Lagrange interpolating function:

def lagrange_interpolant(x, f):
    N = len(x)
    a = np.ones(N)
    for i in range(N):
        for j in range(N):
            if i != j:
                a[i] *= 1/(x[i] - x[j])
    def p(x_):
        result = 0.
        for i in range(N):
            term_i = a[i] * f[i]
            for j in range(N):
                if i != j:
                    term_i *= x_ - x[j]
            result += term_i
        return result
    
    return p

So we give our equidistant grid and the sampled values to this function and we get our interpolating function:

N = 15

def f(x):
    return np.exp(-(5*x)**2) + x

x = np.linspace(-1, 1, N)
f_eq = f(x) 
f_inter = lagrange_interpolant(x, f_eq)

xx = np.linspace(-1, 1, 300)

plt.plot(x, f_eq, 'o')
plt.plot(xx, f_inter(xx))
plt.ylim(None, 2.2)

As you can see, in the middle of the grid, the interpolation is kind of okay, but at the boundaries, it is going insane. This is called the Runge phenomenon and it occurs whenever you try to interpolate a function by polynomials when you are using an equidistant grid. The more grid points you have, the worse the problem. But what is the origin of that problem? The problem is: equidistance. You can avoid the Runge phenomenon by using a suitably non-equidistant grid. One of those grids, which are working very well, is the Chebyshev grid, which we will introduce in the next section.

The solution: Chebyshev grids

We have just seen that equidistant grids can lead to very big errors when interpolating with polynomials. An non-equidistant grid like the Chebyshev grid, however, can solve that problem easily. It is easy to visualize. Suppose you want to sample the domain from -1 to 1. Draw a semicircle above the axis and mark equistant points on that semicircle . Then use the x-coordinates of those points as your computational grid. The situation looks like this:

The points thus defined are located at

As you can see, the points are denser at the boundaries ±1 and look rather equidistant in the middle. Here is a function that gives you a Chebyshev grid for an arbitrary interval:

def chebyshevspace(x_left, x_right, N):
    radius = (x_right - x_left) / 2.
    center = (x_right + x_left) / 2.    
    return center + radius * np.cos(np.pi - np.arange(N+1)*np.pi/N)

Ok, so let’s apply the Lagrange method of interpolation using a Chebyshev grid instead of an equidistant one:

x = chebyshevspace(-1, 1, 20)

f_cheb = f(x) 
f_inter = lagrange_interpolant(x, f_cheb)

xx = np.linspace(-1, 1, 300)

plt.plot(x, f_cheb, 'o')
plt.plot(xx, f_inter(xx))

That looks much better. The Runge phenomenon is gone. But how much better is it really? Let’s vary the number of grid points and compute the error of the interpolating function with the analytically known solution in each case. And we compare equidistant grid and Chebyshev grids. Here is a Python code that is doing the job:

Ns = np.arange(10, 70, 2)

errs_eq = []
errs_cheb = []
xx = np.linspace(-1, 1, 300)

def calc_error(x_k, f_k):
    """Calculate the error of interpolated solution"""
    f_inter = lagrange_interpolant(x_k, f_k)
    return np.max(np.abs(f_inter(xx) - f(xx)))

    
for N in Ns:
    
    # Equidistant polynomial interpolation:    
    x = np.linspace(-1, 1, N)
    
    errs_eq.append(
        calc_error(x, f(x))
    )
    
    # Chebyshev polynomial interpolation:
    x = chebyshevspace(-1, 1, N)
    
    errs_cheb.append(
        calc_error(x, f(x))
    )

plt.semilogy(Ns, errs_eq, label='equidistant', lw=2)
plt.semilogy(Ns, errs_cheb, label='Chebyshev', lw=2)
plt.grid()
plt.xlabel('N')
plt.ylabel('max. interpolation error')
plt.legend();

As you can see, the interpolation with equidistant points is doing a really bad job. At first, the error is increasing with more grid points, then it is decreasing slightly, and the it is diverging wildly. The interpolation with Chebyshev grids, on the other hand, is behaving sensibly. With finer grids, it quickly and monotonously gives smaller errors and we already reach machine precision with about 65 points. Amazing!

Conclusion

In conclusion, Chebyshev grids provide a powerful method for interpolation with extreme accuracy. By using the Chebyshev nodes as the interpolation points, we can achieve high precision even for functions with large variations. The implementation in Python is relatively simple, and the use of the NumPy and SciPy libraries makes it easy to perform calculations. We hope that this article has given you a good understanding of how to use Chebyshev grids for interpolation and has provided you with the tools to implement it in your own projects.