Neumann Boundary Conditions, Decoded

The following function (from L.N. Trefethen, Spectral Methods in MATLAB, with slight modifications) solves the 2nd order wave equation in 2 dimensions (u_{tt} = u_{xx}+u_{yy}) using spectral methods, Fourier for x and Chebyshev for y direction. On its rectangular domain, the equation is subject to Neumann boundary conditions along the sides, u_{y}(x,\pm1,t) = 0, and periodic boundary conditions at the ends, u(+3,y,t) = u(-3,y,t).

function wavetank
%WAVETANK 2D "wave tank" with Neumann BCs for |y| = 1

% x variable in [-A,A], Fourier:
A = 3; Nx = 50; dx = 2*A/Nx; x = -A + dx*(1:Nx)';
D2x = (pi/A)^2*toeplitz([-1/(3*(dx/A)^2) - 1/6 ...

% y variable in [-1,1], Chebyshev:
Ny = 15; [Dy,y] = cheb(Ny); D2y = Dy^2;
BC = -Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny);

% Grid and initial data:
[xx,yy] = meshgrid(x,y);
vv = exp(-8*((xx + 1.5).^2 + yy.^2));
dt = 5/(Nx + Ny^2);
vvold = exp(-8*((xx + dt + 1.5).^2 + yy.^2));

% Time-stepping by leap frog formula:
clf, plotgap = round(2/dt); dt = 2/plotgap;
for n = 0:2*plotgap
    t = n*dt;
    if rem(n + .5,plotgap)<1
        subplot(3,1,n/plotgap + 1), mesh(xx,yy,vv), view(-10,60)
        axis([-A A -1 1 -0.15 1]), colormap(1e-6*[1 1 1]);
        text(-2.5,1,.5,['t = ' num2str(t)],'fontsize',18),
        set(gca,'ztick',[]), grid off, drawnow
    vvnew = 2*vv - vvold + dt^2*(vv*D2x + D2y*vv);
    vvold = vv; vv = vvnew;
    vv([1 Ny+1],:) = BC*vv(2:Ny,:);       % Neumann BCs for |y| = 1

function [D,x] = cheb(N)
%CHEB  compute D = differentiation matrix, x = Chebyshev grid

if N == 0, D = 0; x = 1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X - X';
D  = (c*(1./c)')./(dX + (eye(N+1)));      % off-diagonal entries
D  = D - diag(sum(D,2));                  % diagonal entries

The way Neumann boundary conditions are implemented here (highlighted lines) is elegant, and compressed. With hints from Prof. Trefethen in a private communication, here is a “decompression” of the implementation.

The goal of the statement

vv([1 Ny+1],:) = BC*vv(2:Ny,:)

is to set the data on the top and bottom of the tank to the right values, i.e., values that will make the normal derivative equal to zero. The aim is to end up with values on the top and bottom such that when Dy is applied to do the differentiation, i.e., we compute

Dy * vv(:,:)

the result is zero in rows 1 and Ny+1. In other words, we want

Dy([1 Ny+1],:) * vv(:,:)

to be a matrix consisting of two rows of zeros. Another way to say that is that

Dy([1 Ny+1],[1 Ny+1]) * vv([1 Ny+1],:)

should be the negative of

Dy([1 Ny+1],2:Ny) * vv(2:Ny,:)

In other words,

vv([1 Ny+1],:) = -Dy([1 Ny+1],[1 Ny+1]) \  ( Dy([1 Ny+1,2:Ny]) * vv(2:Ny,:) )

And there in effect we have the matrix BC, i.e.,

BC = -Dy([1 Ny+1],[1 Ny+1])\Dy([1 Ny+1],2:Ny)
vv([1 Ny+1],:) = BC*vv(2:Ny,:)

The key to “decoding” these lines of code is looking at the matrix-matrix multiplication as a linear combination of the rows of the second matrix with the weights coming from the rows of the first, much the same way the application of elimination matrices is interpreted in Gaussian elimination and one of the various interpretations of matrix-matrix or matrix-vector

See also

Lecture 1 of Numerical Linear Algebra (Trefethen and Bau) discusses some of various interpretations of matrix-matrix and matrix-vector multiplications.

This entry was posted in High Order and Spectral Methods and tagged , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s