%pylab inline
import numpy as np
import pylab as plt
Before we look at the spectral collocation method, let's see just why it's important not to use an equispaced grid. We'll try interpolating Runge's function:
$$f(x)=\frac{1}{1+25x^2}$$
on equispaced and Chebyshev grids. We can do this very easily using the NumPy functions polyfit() and poly1d(). Of course, under the hood polyfit() is just solving a certain linear system (a Vandermonde matrix, in fact) associated with the interpolation problem.
The script below plots $f(x)$ and these two interpolants. It also plots the Chebyshev point locations. Try changing $m$ and see what happens. Is the Chebyshev interpolant more accurate at every point?
m=10 # This is the number of (interior) points to be used.
# Runge's function:
f = lambda x: 1./(1+25*x**2)
# The grids:
x_eq = np.linspace(-1,1,m+2)
x_cheb= np.cos(np.pi*np.arange(m+2)/(m+1.))
# The polynomial interpolants:
p1=np.poly1d(np.polyfit(x_eq,f(x_eq),m+1))
p2=np.poly1d(np.polyfit(x_cheb,f(x_cheb),m+1))
# The rest is just plotting:
x_fine=np.linspace(-1,1,1000)
plt.clf()
plt.plot(x_fine,f(x_fine),x_fine,p1(x_fine),x_fine,p2(x_fine))
plt.legend(['f(x)','Equispaced interpolant','Chebyshev interpolant'],loc='best')
plt.plot(x_cheb,f(x_cheb),'ok')
plt.title('N='+str(m+2)+' points')
plt.axis([-1,1,-0.5,1.5])
Now we'll look at implementation of a Chebyshev spectral collocation method. The function below computes a matrix that gives a spectral approximation to the first derivative over the interval $[a,b]$. This is adapted from MATLAB code in Trefethen's text on spectral methods (see Chapter 6 of that text for an explanation). The function also returns the Chebyshev grid points over the interval.
def cheb(m,a,b):
x = np.cos(np.pi*np.arange(m+2)/(m+1)) # Chebyshev points
c = np.ones(m+2); c[0]=2; c[-1]=2;
c = c*np.power(-1,np.arange(m+2))
X = np.tile(x,[m+2,1]).T
dX = X-X.T
D = (np.outer(c,1./c))/(dX+np.eye(m+2))
D = D - np.diag(np.sum(D,axis=1))
D = -D*2./(b-a) # Rescale from [-1,1] to [a,b]
x=a + 0.5*(b-a)*(1.+np.cos(np.pi*(1.-np.arange(m+2)/(m+1.)))) #Shifted Chebyshev points
return D,x
D,x=cheb(1,-1,1)
print D
print x
The Chebyshev differentiation matrix $D$ approximates the first derivative. Let's try using it to approximately compute the first derivative of Runge's function:
m=100
D,x_cheb=cheb(m,-1,1)
F=f(x_cheb)
dF=np.dot(D,F)
plt.clf()
plt.plot(x_cheb,F,x_cheb,dF)
Now let's use $D$ to solve the BVP:
$$\begin{align*} u''(x) = \exp(4x) \quad \quad -1\le x \le1 \\ u(-1)=u(1)=0.\end{align*}$$
To approximate the second derivative, we just take $D^2$. To implement the Dirichlet boundary conditions, we change the first and last rows of the resulting system of equations:
def solvecheb(m,f,alpha,beta):
D,x=cheb(m,-1.,1.);
D2=np.dot(D,D);
D2=D2[1:-1,:];
A=np.zeros([m+2,m+2]); A[1:-1,:]=D2; A[0,0]=1; A[-1,-1]=1;
F=np.zeros(m+2)
F[0]=alpha; F[m+1]=beta
F[1:-1]=f(x[1:-1])
return np.linalg.solve(A,F),x
f = lambda x: np.exp(4.*x)
U,x=solvecheb(15,f,0,0)
plt.clf()
plt.plot(x,U)