Documentation

This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English verison of the page.

Note: This page has been translated by MathWorks. Please click here
To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

Finite Difference Laplacian

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

The Domain

For this example, NUMGRID numbers points within an L-shaped domain. The SPY function is a very useful tool for visualizing the pattern of non-zero elements in a given matrix.

R = 'L'; % Other possible shapes include S,N,C,D,A,H,B
% Generate and display the grid.
n = 32;
G = numgrid(R,n);
spy(G)
title('A finite difference grid')

% Show a smaller version as sample.
g = numgrid(R,12)
g = 

     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     6    11    16    21    26    36    46    56    66     0
     0     2     7    12    17    22    27    37    47    57    67     0
     0     3     8    13    18    23    28    38    48    58    68     0
     0     4     9    14    19    24    29    39    49    59    69     0
     0     5    10    15    20    25    30    40    50    60    70     0
     0     0     0     0     0     0    31    41    51    61    71     0
     0     0     0     0     0     0    32    42    52    62    72     0
     0     0     0     0     0     0    33    43    53    63    73     0
     0     0     0     0     0     0    34    44    54    64    74     0

The Discrete Laplacian

Use DELSQ to generate the discrete Laplacian. The SPY function gives a graphical feel of the population of the matrix.

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

% Number of interior points
N = sum(G(:)>0)
N = 675

The Dirichlet Boundary Value Problem

Finally, we solve the Dirichlet boundary value problem for the sparse linear system. The problem is setup as follows:

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

The Solution

Map the solution onto the grid and show it as a contour map.

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

Now show the solution as a mesh plot.

colormap((cool+1)/2);

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

Was this topic helpful?