A Chebyshev spectral method for advection using Chebfun¶

Developed by Randy LeVeque for a course on Approximation Theory and Spectral Methods at the University of Washington.

See http://faculty.washington.edu/rjl/classes/am590a2013/codes.html for more IPython Notebook examples.

This is a simple example using Forward Euler time-stepping.

The equation solved is

$u_t + u_x=0$ on the interval

$[-1,1]$ with

$u(x,0) = \exp(-50 (x+0.5)^2)$ and boundary condition

$u(0,t) = 0$ .

In [17]:

import pymatbridge
ip = get_ipython()
pymatbridge.load_ipython_extension(ip)

In [18]:

%%matlab
path(path,'/Users/rjl/chebfun/chebfun_v4.2.2889/chebfun')

Set up the matrix and investigate eigenvalues:¶

Set up the Chebyshev points and differentiation matrix. The first point in $x$ and the first row and column of the matrix $D$ are deleted in order to implement the boundary condition. $D$ is also negated since we are solving $u_t = -u_x$ .

In [19]:

%%matlab
N = 100;
x = chebpts(N);
x = x(2:N);
Dop = chebop(@(u) diff(u));
D = Dop(N);
D = -D(2:N,2:N);

In [20]:

%%matlab
[V,Lam] = eig(D);
Lam = diag(Lam);
Lmax = max(abs(Lam))

Lmax =
     8.693384277076647e+02

Choose a time step dt $= \Delta t$ , and then plot the eigenvalues of $\Delta t D$ . Also plot the stability region for forward Euler and zoom in near the origin. If the time step is small enough, the points should lie inside the circle.

In [21]:

%%matlab
nmax = 2000;  % number of time steps to reach t=1
dt = 1./nmax;
plot(real(Lam*dt), imag(Lam*dt),'r.')
eitheta = exp(1i*linspace(0,2*pi,1000));
hold on
plot(eitheta-1,'k')
axis 'equal'
xylim = dt*Lmax + 0.1;
axis([-xylim,0.1,-xylim,xylim])

Pseudospectra¶

Since the matrix D is highly non-normal, we should also check that the $\lambda \Delta t$ lies in the stability region for $\lambda$ in the pseudospectrum of $D$ for small $\epsilon$ ...

In [22]:

%%matlab
x1 = linspace(-xylim,0.1,100);
y1 = linspace(-xylim,xylim,200);
[xx,yy] = meshgrid(x1,y1);
zz = xx + 1i*yy;
I = eye(size(D,1));
sigmin = zeros(size(zz));
for j=1:length(x1)
    for i=1:length(y1)
        sigmin(i,j) = min(svd(zz(i,j)*I - dt*D));
        end
    end

plot(real(Lam*dt), imag(Lam*dt),'r.')
hold on
plot(eitheta-1,'k')
axis 'equal'
axis([-xylim,0.1,-xylim,xylim])
contour(xx,yy,sigmin,[1e-2,1e-3,1e-4],'b')
contour(xx,yy,sigmin,[1e-5,1e-6,1e-7],'k')
contour(xx,yy,sigmin,[1e-8,1e-9,1e-10],'r')

Initial conditions:¶

In [23]:

%%matlab
u = exp(-50*(x+0.5).^2);
plot(x,u,'b');

Time-stepping loop:¶

In [24]:

%%matlab
for n=1:nmax
    u = u + dt*D*u;
    end

plot(x,u,'r-o')
hold on
utrue = exp(-50*(x-0.5).^2);
plot(x,utrue,'b')

Compute the maximum error at the gridpoints.

In [25]:

%%matlab
norm(u-utrue,inf)

ans =
   0.025976954585357

Alternatively we can interpolate with a polynomial and compute the max-norm of the error in this function:

In [26]:

%%matlab
un = [0; u];
pn = chebfun(u)
phat = chebfun('exp(-50*(x-0.5).^2)')
max_error = norm(pn-phat)

pn = 
   chebfun column (1 smooth piece)
       interval       length   endpoint values   
[      -1,       1]       99 -5.6e-08  2.9e-06   
vertical scale =   1 
phat = 
   chebfun column (1 smooth piece)
       interval       length   endpoint values   
[      -1,       1]       86  1.4e-49  3.7e-06   
vertical scale =   1 
max_error =
   0.029392125102472

Try changing N or nmax in the cells above and re-executing everything to see the effect of changing the resolution or time step.

Leapfrog¶

We do not expect Leapfrog to work well for this problem since the eigenvalues are not pure imaginary. Even with a small dt the solution soon blows up.

In [27]:

%%matlab
un = exp(-50*(x+0.5).^2);
plot(x,un,'b');

dt = 1e-5;
nmax = 3000;
unm1 = exp(-50*(x+0.5-dt).^2);

for n=1:nmax
    unp1 = unm1 + 2*dt*D*un;
    unm1 = un;
    un = unp1;
    end
    
hold on
plot(x,un,'r-o')

Note that there are unstable eigenvalues for which the eigenvector looks like this...

In [28]:

%%matlab
k = 10;
plot(x,real(V(:,k)))
Lam(k)

ans =
     -1.471839354591068e+02 - 1.608047492945914e+01i

In [28]: