Interpolation¶
Given a set of $n+1$ data points $(x_0, y_0), (x_1, y_1), \ldots, (x_n, y_n)$, we want to find a polynomial $P(x)$ of degree at most $n$ such that $P(x_i) = y_i$ for $i = 0, 1, \ldots, n$. This is called interpolation.
Chebyshev Interpolation¶
Usually we use the Lagrange interpolation formula to find the interpolating polynomial. However, the Lagrange interpolation formula is numerically unstable when the data points are closely spaced. In this case, we can use the Chebyshev interpolation formula to find the interpolating polynomial.
Chebyshev Points¶
The Chebyshev points are defined as $$ x_i = \cos\left(\frac{2i+1}{2(n+1)}\pi\right), \quad i = 0, 1, \ldots, n. $$ The Chebyshev points are the roots of the Chebyshev polynomial of the first kind $T_{n+1}(x) = \cos((n+1)\arccos(x))$.
Chebyshev Interpolation Formula¶
The Chebyshev interpolation formula is given by $$ P(x) = \sum_{i=0}^n y_i T_i(x), $$ where $T_i(x)$ is the Chebyshev polynomial of the first kind of degree $i$.
Example¶
Let's consider the following data points: $$ (0,1), (1,2), (2,3), (3,4), (4,5). $$ We want to find the interpolating polynomial using the Chebyshev interpolation formula.
Implementation¶
We first define the Chebyshev polynomial of the first kind and the Chebyshev interpolation formula. Then we find the interpolating polynomial using the Chebyshev interpolation formula.
from numpy import cos, arccos, pi
def T(n, x):
return cos(n*arccos(x))
def chebyshev_interpolation(x, y):
n = len(x) - 1
p = 0
for i in range(n+1):
p += y[i]*T(i, x)
return p
x = [0, 1, 2, 3, 4]
y = [1, 2, 3, 4, 5]
p = chebyshev_interpolation(x, y)
print(p)
# Suppose now we want to use chebyshev interpolation for a custom function in the interval [-1, 1]
# We will use the Chebyshev nodes and the Chebyshev polynomials of the first kind
# The Chebyshev nodes are defined as x_i = cos((2i-1)*pi/(2n)), i = 1, 2, ..., n
# The Chebyshev polynomials of the first kind are defined as T_i(x) = cos(i*arccos(x)), i = 0, 1, 2, ...
# The Chebyshev interpolation of a function f(x) in the interval [-1, 1] is given by
# p(x) = sum_{i=0}^{n} a_i*T_i(x), where a_i = (2/n)*sum_{j=1}^{n} f(x_j)*T_i(x_j)/(1 - delta_{i0})
# delta_{ij} is the Kronecker delta
# We will use the numpy library to perform the calculations
import numpy as np
import matplotlib.pyplot as plt
# Define the custom function
def func(x):
return np.tanh(np.exp(x))+np.sin(x)**2-4*x*np.cos(np.exp(4*x**2))
y = np.polynomial.chebyshev.chebfit(np.linspace(-1, 1, 1000), func(np.linspace(-1, 1, 1000)), 14) # 14 is the degree of the polynomial approximation
p = np.polynomial.chebyshev.Chebyshev(y) # p is the Chebyshev interpolation
x = np.linspace(-1, 1, 1000)
plt.figure(figsize=(12,8))
plt.title('Chebyshev Interpolation for a Custom Function')
plt.scatter(x, func(x), color = 'salmon', label='Custom Function')
# Plot the Chebyshev interpolation with different color
plt.plot(x, p(x), color = 'teal', label='Chebyshev Interpolation')
plt.legend()
plt.show()
# Suppose now we use the Lagrange interpolation for the same custom function just for comparison
# We will use the numpy library to perform the calculations
import scipy.interpolate as spi # SciPy
y = func(np.linspace(-1, 1, 10))
p2 = spi.lagrange(np.linspace(-1, 1, 10), y)
plt.figure(figsize=(12,8))
plt.title('Lagrange Interpolation and Chebyshev Interpolation for a Custom Function')
plt.scatter(x, func(x), color = 'salmon', label='Custom Function')
plt.plot(x, p(x), color='teal', label='Chebyshev interpolation')
plt.plot(x, p2(x), color='navy', label='Lagrange interpolation')
plt.legend()
plt.show()
As we can see, Chebyshev Interpolation and Lagrange Interpolation gives different results. Since the data points created from the custom function are not easy to interpret, and the ill-conditioned nature of the Vandermonde matrix, the Chebyshev Interpolation is more stable and accurate than Lagrange Interpolation.
However, one may wonder why we can use some polynomials to interpolate the data points (which is the case of interpolation) and not others. One may refer to the Weierstrass Approximation Theorem for more information.
The Theorem is stated as follows:
Weierstrass Approximation Theorem¶
Let $f(x)$ be a continuous function on $[a, b]$. Then for any $\epsilon > 0$, there exists a polynomial $P(x)$ such that $|f(x) - P(x)| < \epsilon$ for all $x \in [a, b]$.
This theorem guarantees the existence of a polynomial that can approximate any continuous function on a closed interval. The Chebyshev Interpolation is one of the methods to find such a polynomial.
Here comes another question: Is it true that we use higher-degree polynomials to approximate the function better? Let us think of the following example:
import numpy as np
import matplotlib.pyplot as plt
# Define the custom function
def func(x):
return 1/(1 + 25*x**2)
# Chebyshev Interpolation for different number of nodes
plt.figure(figsize=(12,8),dpi=100)
plt.title('Chebyshev Interpolation for a Custom Function')
x = np.linspace(-1, 1, 1000)
plt.plot(x, func(x), color = 'salmon', label='Custom Function')
for n in [5, 10, 15]:
y = func(np.linspace(-1, 1, n))
p2 = np.polynomial.chebyshev.chebfit(np.linspace(-1, 1, n), func(np.linspace(-1, 1, n)), n-1)
plt.plot(x, np.polynomial.chebyshev.Chebyshev(p2)(x), label='Chebyshev Interpolation for n = {}'.format(n), linewidth=2, linestyle='--', alpha=0.7)
plt.legend()
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.show()
Therefore, the answer to the question is quite obvious. The higher-degree polynomial does not always give a better approximation. Also, values at the endpoints are even far from the actual values. This is called Runge's phenomenon.
Runge's Phenomenon¶
Runge's phenomenon is the oscillation of the interpolating polynomial at the endpoints when using higher-degree polynomials. This phenomenon is due to the fact that the higher-degree polynomial oscillates more than the lower-degree polynomial. Therefore, the higher-degree polynomial does not always give a better approximation.
# Here we use Chebyshev Nodes and Chebyshev Interpolation to see whether non-uniformly distributed nodes can improve the approximation
# Define the custom function
def func(x):
return 1/(1 + 25*x**2)
# Chebyshev Nodes
def cheb_nodes(n):
return np.cos((2*np.arange(1, n+1)-1)*np.pi/(2*n))
# Chebyshev Interpolation for different number of nodes
plt.figure(figsize=(12,8),dpi=100)
plt.title('Chebyshev Interpolation for a Custom Function')
x = np.linspace(-1, 1, 1000)
plt.plot(x, func(x), color = 'salmon', label='Custom Function')
for n in range(5, 31, 5):
y = func(cheb_nodes(n))
p2 = np.polynomial.chebyshev.chebfit(cheb_nodes(n), y, n-1)
plt.plot(x, np.polynomial.chebyshev.Chebyshev(p2)(x), label='Chebyshev Interpolation for n = {}'.format(n), linewidth=2, linestyle='--', alpha=0.7)
plt.legend()
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.show()
The result now seems better than the previous one. The difference we make here is to use Chebyshev nodes instead of equispaced nodes. The Chebyshev nodes are defined as $x_i = \cos\left(\frac{2i+1}{2(n+1)}\pi\right)$, $i = 0, 1, \ldots, n$. The Chebyshev nodes are the roots of the Chebyshev polynomial of the first kind $T_{n+1}(x) = \cos((n+1)\arccos(x))$. Therefore, we may conclude that polynomial interpolation in equispaced points is significantly ill-conditioned (One may know that it is exponentially ill-conditioned, but to discuss it analytically is out of the scope of this notebook). Therefore, we may use Chebyshev nodes to avoid this ill-conditioned nature of the interpolation.
Conclusion¶
In this notebook, we have discussed the Chebyshev Interpolation and its implementation. We have also discussed the Weierstrass Approximation Theorem and Runge's Phenomenon. Chebyshev Interpolation is more stable and accurate than Lagrange Interpolation when the data points are closely spaced. However, the higher-degree polynomial does not always give a better approximation.
References¶
- Lloyd N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2013. (I find this book interesting to read; though, it uses MATLAB instead of Python)