Main Content

Finite Difference Laplacian

This example shows how to compute and represent the finite difference Laplacian on an L-shaped domain.

Domain

The numgrid function numbers points within an L-shaped domain. The spy function is a useful tool for visualizing the pattern of nonzero elements in a matrix. Use these two functions to generate and display an L-shaped domain.

n = 32;
R = 'L';
G = numgrid(R,n);
spy(G)
title('A Finite Difference Grid')

Figure contains an axes object. The axes object with title A Finite Difference Grid, xlabel nz = 675 contains a line object which displays its values using only markers.

Show a smaller version of the matrix as a sample.

g = numgrid(R,10)
g = 10×10

     0     0     0     0     0     0     0     0     0     0
     0     1     5     9    13    17    25    33    41     0
     0     2     6    10    14    18    26    34    42     0
     0     3     7    11    15    19    27    35    43     0
     0     4     8    12    16    20    28    36    44     0
     0     0     0     0     0    21    29    37    45     0
     0     0     0     0     0    22    30    38    46     0
     0     0     0     0     0    23    31    39    47     0
     0     0     0     0     0    24    32    40    48     0
     0     0     0     0     0     0     0     0     0     0

Discrete Laplacian

Use delsq to generate the discrete Laplacian. Use the spy function again to get a graphical feel of the matrix elements.

D = delsq(G);
spy(D)
title('The 5-Point Laplacian')

Figure contains an axes object. The axes object with title The 5-Point Laplacian, xlabel nz = 3255 contains a line object which displays its values using only markers.

Determine the number of interior points.

N = sum(G(:)>0)
N = 
675

Dirichlet Boundary Value Problem

Solve the Dirichlet boundary value problem for the sparse linear system. The problem setup is:

delsq(u) = 1 in the interior, u = 0 on the boundary.

rhs = ones(N,1);
if (R == 'N') % For nested dissection, turn off minimum degree ordering.
   spparms('autommd',0)
   u = D\rhs;
   spparms('autommd',1)
else
   u = D\rhs; % This is used for R=='L' as in this example
end

Map the solution onto the L-shaped grid and plot it as a contour map.

U = G;
U(G>0) = full(u(G(G>0)));
clabel(contour(U));
prism
axis square ij

Figure contains an axes object. The axes object contains 15 objects of type contour, line, text.

Now show the solution as a mesh plot.

mesh(U)
axis([0 n 0 n 0 max(max(U))])
axis square ij

Figure contains an axes object. The axes object contains an object of type surface.

See Also