The following function (from L.N. Trefethen, Spectral Methods in MATLAB, with slight modifications) solves the 2nd order wave equation in 2 dimensions () 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, , and periodic boundary conditions at the ends, .

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.