Chebyshev spectral method¶
AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.
These notebooks are all available on Github.
This notebook illustrates:
- Polynomial interpolation using monomial and Chebyshev polynomial bases
- The dangers of interpolating at equally spaced points using the Runge function
- How to do Chebyshev interpolation
- Approximating the first derivative of a function using the Chebyshev points.
%matplotlib inline
from pylab import *
Interpolation using the monomial basis¶
First we illustrate polynomial interpolation using the monomial basis $$ p(x) = c_0 + c_1 x + c_2 x^2 + \cdots c_{N-1}x^{N-1} + c_N x^N $$ Note that this can be rewritten as $$ p(x) = ((\cdots (c_N x + c_{N-1})x + c_{N-2})x + \cdots + c_1)x + c_0 $$ A nice way to evaluate $p(x)$ that is sometimes called Horner's method.
In the function below the Vandermonde matrix $B$ is set up, with $k$th column given by the values $x_j^k$ where $x_j$ are the desired interpolation points. The system $$ Bc = u $$ is then solved for the polynomial coefficients, where $u$ is the vector of function values at each point.
The resulting polynomials is then evaluated on a fine grid of points (using Horner's method) for plotting purposes and to estimate the error.
Any set of $N+1$ distinct points $x$ could be used, here it's set up to use either equally spaced or Chebyshev points.
def poly_interp_monomials(u_fcn, x):
"""
Peform polynomial interpolation of the function u_fcn at points x using the monomial basis.
"""
# Vandermonde matrix using monomial basis:
N = len(x) - 1
B = empty((N+1,N+1), dtype=float) # initialize storage for B
for k in range(N+1):
B[:,k] = x**k
print('The condition number of the Vandermonde matrix B is %g' % cond(B))
# Solve for polynomial coefficients:
u = u_fcn(x)
c = solve(B,u)
return c
def poly_plot_monomials(u_fcn, x, c):
"""
Plot the true function and the interpolating function
"""
xfine = linspace(-1,1,1000)
ufine = u_fcn(xfine)
figure(figsize=(10,4))
plot(xfine, ufine, 'b', label='Original function')
plot(x, u_fcn(x), 'bo', label="Interpolation points")
# Evaluate p(xfine) using Horner's method:
p = c[N]
for k in range(N-1,-1,-1):
p = p*xfine + c[k]
plot(xfine, p, 'r', label='Polynomial interpolant')
plot(x, zeros(x.shape), 'k-+')
legend()
grid(True)
title('Using %i interpolation points' % N)
print('Max absolute error = %g' % abs(p-ufine).max())
Test on the Runge function¶
This illustrate the danger of polynomial interpolation on a uniform grid.
u_fcn = lambda x: 1 / (1 + 16*x**2)
N = 12
x = linspace(-1,1,N+1)
c = poly_interp_monomials(u_fcn, x)
poly_plot_monomials(u_fcn, x, c)
If you increase $N=12$ to a larger value you will see that the approximation gets much worse near the boundaries.
Using Chebyshev points¶
If we instead use the Chebyshev interpolation points we get much better results:
N = 12
x = cos(linspace(0, pi,N+1))
c = poly_interp_monomials(u_fcn, x)
poly_plot_monomials(u_fcn, x, c)
N = 24
x = cos(linspace(0, pi,N+1))
c = poly_interp_monomials(u_fcn, x)
poly_plot_monomials(u_fcn, x, c)
Although this works fine for relatively small $N$, note that the condition number of the Vandermonde matrix is getting very large. That's because the monomial basis is poorly conditioned.
It is much better to use the basis of Chebyshev polynomials.
Evaluating Chebyshev polynomials¶
We can use the 3-term recurrence relation to evaluate the Chebyshev polynomial $T_k(x)$ for any order $k$ at a set of points x (assumed to be a numpy array) as in the function below.
The recurrence is: $$ T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x) $$ with starting values: $$ T_0(x) = 1, \quad T_1(x) = x. $$
def T(n,x):
"""
Evaluate Chebyshev polynomial T_n(x) using 3-term recurrence relation.
Assumes x is a numpy array of values.
Does not work as written if x is just a single value!
"""
# Initialize Tnm = 1 and Tn = x and then iterate
Tnm = ones(x.shape) # will be used for T_{n-1} in loop
if n==0:
return Tnm # special case T_0 = 1 forall x
Tn = x # will be used for T_n in loop
for k in range(2,n+1):
Tnp = 2*x*Tn - Tnm
Tnm = Tn.copy()
Tn = Tnp.copy()
return Tn
xfine = linspace(-1,1,1000)
k = 7
plot(xfine, T(k,xfine), 'b')
grid(True)
title('Chebyshev polynomial T_%s' % k);
Evaluating a polynomial written in the Chebyshev basis¶
Once we have computed coefficients $c$ by solving an interpolation problem in the Cheyshev basis, we need to be able to evaluate $$ p(x) = c_0T_0(x) + \cdots c_NT_N(x). $$ This routine performs this sum using the 3-term recurrence:
def T_sum(c,x):
"""
Evaluate sum c_n * T_n(x) (from n=0 to N) using 3-term recurrence relation
"""
N = len(c) - 1
Tnm = ones(x.shape) # will be used for T_{n-1} in loop
Tsum = c[0]*Tnm
if N==0:
return Tsum # special case T_0 = 1 forall x
Tn = x # will be used for T_n in loop
Tsum += c[1]*Tn
for n in range(2,N+1):
Tnp = 2*x*Tn - Tnm
Tsum += c[n]*Tnp
Tnm = Tn.copy()
Tn = Tnp.copy()
return Tsum
Another way to evaluate the polynomial at an arbitrary set of points such as xfine is to set up a matrix $B$ with $j$ th column $T_j(x)$ and the multiply $Bc$ where $c$ are the coefficients. Note this matrix can be rectangular $m \times (N+1)$ where $m$ is the number of points where it should be evaluated.
This idea will be useful for evaluating higher order derivatives, so we use this instead of the above:
def T_sum(c,x):
N = len(c) - 1
B = empty((len(x),N+1))
for k in range(N+1):
B[:,k] = T(k,x)
Tsum = dot(B,c)
return Tsum
Chebyshev interpolation¶
The function below sets up a Vandermonde matrix using the Chebyshev basis, solves for the coefficients c, and then evaluates the resulting polynomial on a fine grid of points for plotting purposes and to estimate the error:
def poly_interp_chebyshev(u_fcn, x):
# Vandermonde matrix using monomial basis:
N = len(x) - 1
B = empty((N+1,N+1), dtype=float) # initialize storage for B
for k in range(N+1):
B[:,k] = T(k,x)
print('The condition number of the Vandermonde matrix B is %g' % cond(B))
# Solve for polynomial coefficients:
u = u_fcn(x)
c = solve(B,u)
return c
def poly_plot_chebyshev(u_fcn, x, c):
xfine = linspace(-1,1,1000)
ufine = u_fcn(xfine)
figure(figsize=(10,4))
plot(xfine, ufine, 'b', label='Original function')
plot(x, u_fcn(x), 'bo', label="Interpolation points")
# Evaluate p(xfine) using T_sum:
p = T_sum(c,xfine)
plot(xfine, p, 'r', label='Polynomial interpolant')
plot(x, zeros(x.shape), 'k-+')
legend()
grid(True)
title('Using %i interpolation points' % N)
print('Max absolute error = %g' % abs(p-ufine).max())
Test this on the Runge function. Note that this should give the same polynomial approximation as obtained using the monomial basis, but with a much better conditioned Vandermonde matrix.
u_fcn = lambda x: 1 / (1 + 16*x**2)
N = 24
x = cos(linspace(0, pi,N+1))
c = poly_interp_chebyshev(u_fcn, x)
poly_plot_chebyshev(u_fcn, x, c)
Evaluating the derivative¶
First we define a function that evaluates the derivative $T_n'(x)$. It uses a 3-term recurrence derived from that of the $T_n$ by differentiating: $$ T'_{n+1}(x) = 2T_n(x) + 2xT_n'(x) - T_{n-1}'(x) $$ with starting values: $$ T_0'(x) = 0, \quad T_1'(x) = 1. $$
def Tx(n,x):
"""
Evaluate first derivative of Chebyshev polynomial T_n(x) using 3-term recurrence relation
Assumes x is a numpy array of points.
"""
Tnm = ones(x.shape) # will be used for T_{n-1} in loop
Txnm = zeros(x.shape)
if n==0:
return Txnm # special case T_0 = 1 forall x
Tn = x # will be used for T_n in loop
Txn = ones(x.shape)
for k in range(2,n+1):
Tnp = 2*x*Tn - Tnm
Txnp = 2*Tn + 2*x*Txn - Txnm
Tnm = Tn.copy()
Tn = Tnp.copy()
Txnm = Txn.copy()
Txn = Txnp.copy()
return Txn
Smooth test problem¶
Try it out on $u(x) = \sin(x)$:
u_fcn = lambda x: sin(2*x)
ux_fcn = lambda x: 2*cos(2*x)
N = 4
x = cos(linspace(0, pi, N+1))
c = poly_interp_chebyshev(u_fcn, x)
B = empty((len(xfine),N+1))
# evaluate sum c_n * Tx_n(x) at xfine points:
for k in range(N+1):
B[:,k] = Tx(k,xfine)
pfine = dot(B,c)
# true u'(x) at xfine points:
ux_fine = ux_fcn(xfine)
figure(figsize=(10,6))
plot(xfine, ux_fine, 'b', label="u'(x)")
plot(xfine, pfine, 'r', label="p'(x)")
plot(x, ux_fcn(x), 'bo')
plot(x, zeros(x.shape), 'k-+')
grid(True)
legend()
title('Derivative of u and polynomial interpolant')
print('Max absolute error in derivative = %g' % abs(pfine-ux_fine).max());
With $N=4$ there is a noticable error but increasing to $N=5$ shows much better agreement and $N=18$ gives nearly full machine precision.
Another test:¶
Here's a more oscillatory function where a few more points are needed for good accuracy.
u_fcn = lambda x: (x-0.5)*sin(10*x)
ux_fcn = lambda x: sin(10*x) + 10*(x-0.5)*cos(10*x)
N = 40
x = cos(linspace(0, pi, N+1))
c = poly_interp_chebyshev(u_fcn, x)
B = empty((len(xfine),N+1))
# evaluate sum c_n * Tx_n(x) at xfine points:
for k in range(N+1):
B[:,k] = Tx(k,xfine)
pfine = dot(B,c)
# true u'(x) at xfine points:
ux_fine = ux_fcn(xfine)
figure(figsize=(10,6))
plot(xfine, ux_fine, 'b', label="u'(x)")
plot(xfine, pfine, 'r', label="p'(x)")
plot(x, ux_fcn(x), 'bo')
plot(x, zeros(x.shape), 'k-+')
grid(True)
legend()
title('Derivative of u and polynomial interpolant')
print('Max absolute error in derivative = %g' % abs(pfine-ux_fine).max());
Increasing $N$ to around 40 gives error around $10^{-12}$.