Pseudospectral Approximation of Functions
My graduate research at BYU relied on good boundary value problem (BVP) solvers. So I had the chance to review a lot of what was already implemented in Python, and build several of my own implementations. One approach that I've really enjoyed learning about is the pseudospectral method. This blog post will focus on approximating functions with Chebychev polynomials. I plan to eventually discuss a pseudospectral BVP solver.
Any smooth function defined on $[-1,1]$ can be expanded as a infinite series of Chebychev polynomials $f(x) = \sum_{n=0}^{\infty} a_n T_n(x),$ where the Chebychev polynomials are defined by a recurrence relation: $T_0 = 1$, $T_1 = x$, with $T_{n+1} = 2xT_n(x) - T_{n-1}(x)$. In numerical calculations we always truncate the sum at some degree $N$: $f \approx \sum_{n=0}^{N} a_n T_n(x)$.
Another way to approximate $f$ is to use polynomial interpolation. Polynomial interpolation has problems if an evenly spaced set of grid points is used. A much better set of grid points is the Gauss-Lobatto grid, which consists of the Chebychev extremum plus the endpoints $\pm 1$.
Let's approximate the function $f(x) = x \cos(5x^2 + 2x +1)$. We begin by plotting the function:
import numpy as np
from scipy.interpolate import BarycentricInterpolator
import matplotlib.pyplot as plt
def f(x):
return x*np.cos(5*x**2 + 2*x + 1)
X = np.linspace(-1,1,200)
plt.plot(X,f(X),'-k')
plt.show()
The Chebychev extrema plus endpoints grid can be calculated quickly. We use SciPy's BarycentricInterpolator to plot the interpolating polynomials over grids of increasing size:
N = 5
x1 = np.cos((np.pi/N)*np.linspace(0,N,N+1))
interp = BarycentricInterpolator(x1,f(x1))
u1 = interp.__call__(X)
N = 10
x2 = np.cos((np.pi/N)*np.linspace(0,N,N+1))
interp = BarycentricInterpolator(x2,f(x2))
u2 = interp.__call__(X)
N = 15
x3 = np.cos((np.pi/N)*np.linspace(0,N,N+1))
interp = BarycentricInterpolator(x3,f(x3))
u3 = interp.__call__(X)
plt.plot(X,f(X),'-k',label="$f$")
plt.plot(x1,f(x1),'*g')
plt.plot(X,u1,'-g',label="$p_5$")
plt.plot(x2,f(x2),'*b')
plt.plot(X,u2,'-b',label="$p_{10}$")
plt.legend()
plt.show()
The errors in these approximating polynomials get small remarkably fast:
plt.plot(X,(u1-f(X)),'-g',label="$E_5$")
plt.plot(X,(u2-f(X)),'-b',label="$E_{10}$")
plt.plot(X,(u3-f(X)),'-r',label="$E_{15}$")
plt.ylabel("Error")
plt.xlabel("$x$")
plt.legend()
plt.show()
Math Details: Is Interpolation Safe?
Uniformly spaced interpolation points do cause lots of problems! Uniformly spaced nodes make polynomial interpolation ill-conditioned for any interpolation algorithm. Rounding errors can then grow exponentially. However, Chebychev points provide a good set of interpolation points because of how they are distributed throughout the interval $[-1,1]$. Lloyd Trefethen has an excellent book that discusses interpolation over Chebychev points, Approximation Theory and Approximation Practice.
We used SciPy's BarycentricInterpolator, which is based on a paper by Berrut and Trefethen (2004), Barycentric Lagrange Interpolation. Barycentric interpolation is a stable algorithm with several interesting properties. It does not compute the coefficients of the interpolating polynomial. Given interpolating data $(x_i,y_i)_i$, the interpolator is initialized in $\mathcal{O}(n^2)$ time and evaluated in $\mathcal{O}(n)$ time. However, the initialization does not depend on the function data $y_i$, and so it is possible to update $y_i$ quickly.
The Chebychev interpolator in general has exponential convergence (faster than $n^k$ for every $k$). John Boyd has an excellent summary of these results, and discussions of related numerical methods, in his book Chebychev and Fourier Spectral Methods.