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 ...
    .5*(-1).^(2:Nx)./sin((pi*dx/A)*(1:Nx-1)/2).^2]);

% 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
    end
    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
end

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
multiplication.

See also

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

Advertisements
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:

WordPress.com Logo

You are commenting using your WordPress.com 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 )

Google+ photo

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

Connecting to %s