This example shows how to solve the eigenvalue problem of the Laplace operator on an L-shaped region.
Consider a membrane that is fixed at the boundary of a region in the plane. Its displacement is described by the eigenvalue problem , where is the Laplace operator and is a scalar parameter. The boundary condition is for all .
The Laplace operator is self-adjoint and negative definite, that is, only real negative eigenvalues exist. There is a maximal (negative) discrete eigenvalue, the corresponding eigenfunction is called the ground state. In this example, is an L-shaped region, and the ground state associated with this region is the L-shaped membrane that is the MATLAB® logo.
The simplest approach to the eigenvalue problem is to approximate the Laplacian by a finite difference approximation (a stencil) on a square grid of points with distances hx
in direction and distances hy
in direction. In this example, approximate with a sum S_h
of nine regular grid points around the midpoint . The unknowns are the weights .
syms u(x,y) Eps a11 a10 a1_1 a01 a00 a0_1 a_11 a_10 a_1_1 syms hx hy positive S_h = a_11 * u(x - Eps*hx,y + Eps*hy) +... a01 * u(x,y + Eps*hy) +... a11 * u(x + Eps*hx,y + Eps*hy) +... a_10 * u(x - Eps*hx,y) +... a00 * u(x,y) +... a10 * u(x + Eps*hx,y) +... a_1_1* u(x - Eps*hx,y - Eps*hy) +... a0_1 * u(x,y - Eps*hy) +... a1_1 * u(x + Eps*hx,y - Eps*hy);
Use the symbolic parameter Eps
to sort the expansion of this expression by powers of hx
and hy
. Knowing the weights, you can approximate the Laplacian by setting Eps = 1
.
t = taylor(S_h, Eps, 'Order', 7);
Use the coeffs
function to extract their coefficients of terms with the same powers of Eps
. Each coefficient is an expression containing powers of hx
, hy
, and derivatives of u
with respect to and . Since S_h
represents , the coefficients of all other derivatives of u
must be zero. Extract the coefficients by replacing all derivatives of u
, except and , by 0. Replace and by 1. This reduces the Taylor expansion to the coefficient you want to compute, and leads to the following six linear equations.
C = formula(coeffs(t, Eps, 'All'));
eq0 = subs(C(7),u(x,y),1) == 0;
eq11 = subs(C(6),[diff(u,x),diff(u,y)],[1,0]) == 0;
eq12 = subs(C(6),[diff(u,x),diff(u,y)],[0,1]) == 0;
eq21 = subs(C(5),[diff(u,x,x),diff(u,x,y),diff(u,y,y)],[1,0,0]) == 1;
eq22 = subs(C(5),[diff(u,x,x),diff(u,x,y),diff(u,y,y)],[0,1,0]) == 0;
eq23 = subs(C(5),[diff(u,x,x),diff(u,x,y),diff(u,y,y)],[0,0,1]) == 1;
Since there are nine unknown weights in S_h
, add further equations by requiring that all third-order derivatives of u
are 0.
eq31 = subs(C(4),[diff(u,x,x,x),diff(u,x,x,y),diff(u,x,y,y),diff(u,y,y,y)], [1,0,0,0]) == 0; eq32 = subs(C(4),[diff(u,x,x,x),diff(u,x,x,y),diff(u,x,y,y),diff(u,y,y,y)], [0,1,0,0]) == 0; eq33 = subs(C(4),[diff(u,x,x,x),diff(u,x,x,y),diff(u,x,y,y),diff(u,y,y,y)], [0,0,1,0]) == 0; eq34 = subs(C(4),[diff(u,x,x,x),diff(u,x,x,y),diff(u,x,y,y),diff(u,y,y,y)], [0,0,0,1]) == 0;
Solve the resulting ten equations for the nine unknown weights. Use ReturnConditions
to find all solutions including arbitrary parameters.
[a11,a10,a1_1,a01,a00,a0_1,a_11,a_10,a_1_1,parameters,conditions] = ... solve([eq0,eq11,eq12,eq21,eq22,eq23,eq31,eq32,eq33,eq34], ... [a11,a10,a1_1,a01,a00,a0_1,a_11,a_10,a_1_1], ... 'ReturnConditions',true); expand([a_11,a01,a11;... a_10,a00,a01;... a1_1,a0_1,a_1_1])
ans =
parameters
parameters =
Use the subs
function to replace the weights by their computed values.
C = simplify(subs(C));
The expressions C(7)
, C(6)
, and C(4)
containing the 0th, 1st, and 3rd derivatives of u
vanish.
[C(7), C(6), C(4)]
ans =
The expression C(5)
is the Laplacian of u
.
C(5)
ans =
Thus, with the values of the weights computed above, the stencil S_h
approximates the Laplacian up to order hx^2
, hy^2
for any values of the arbitrary parameter z
, provided that z
is chosen to be of order O(1/hx^2,1/hy^2)
.
Although the solution contains a free parameter z
, the expression C(3)
containing the fourth-order derivatives of u
cannot be turned into zero by a suitable choice of z
. Another option is to turn it into a multiple of the square of the Laplace operator.
syms d
Laplace = @(u) laplacian(u,[x,y]);
expand(d*Laplace(Laplace(u)))
ans(x, y) =
Pick different derivatives of u
in C(3)
, and equate their coefficients with the corresponding terms.
subs(C(3),[diff(u,x,x,x,x),diff(u,x,x,y,y),diff(u,y,y,y,y)],[1,0,0]) == d
ans =
subs(C(3),[diff(u,x,x,x,x),diff(u,x,x,y,y),diff(u,y,y,y,y)],[0,1,0]) == 2*d
ans =
subs(C(3),[diff(u,x,x,x,x),diff(u,x,x,y,y),diff(u,y,y,y,y)],[0,0,1]) == d
ans =
Therefore, you can choose d = hx^2/12 = hy^2/12
and z = 2*d/(hx^2*hy^2)
, implying that hx = hy
and z = 1/(6*hx*hy)
. Hence, the stencil S_h
approximates a modified Laplacian on a square grid with hx = hy = h
.
syms h
hx = h;
hy = h;
d = h^2/12;
Replace hx
and hy
by h
.
C = subs(C);
Replace z
by its value, 1/(6*h^2)
. Because z
does not exist in the MATLAB® workspace, you can access it only as the value stored in the parameters
array.
C = subs(C,parameters,1/(6*h^2));
Verify the formula (1).
simplify(C(3) - d*Laplace(Laplace(u)))
ans(x, y) =
Now, consider the third-order terms in hx
, hy
.
simplify(C(2))
ans =
Since no such terms exist in the expansion of the stencil, the term in (1) is in fact of order . Consider the fourth-order terms of the stencil.
factor(simplify(C(1)))
ans =
Check if these terms can be identified with yet another power of the Laplace operator. However, comparison with
Laplace(Laplace(Laplace(u)))
ans(x, y) =
shows that the expressions of order cannot be identified as some multiple of the third power of the Laplace operator because the coefficients cannot be matched.
For a square grid with distance h
between neighboring grid points and the above choices of the weights, you get:
Use this expansion for the numeric approach to the eigenvalue problem . Add some multiple of to the eigenvalue equation.
The left side of this equation is well approximated by the stencil . Thus, using (2), a numerical eigenvalue of the stencil satisfying must be an approximation of the eigenvalue of the Laplace operator with
For given , solve for to obtain a better approximation of the Laplace eigenvalue. Note that in the solution of the quadratic equation in the correct sign of the square root is given by the requirement that for .
Consider an L-shaped region consisting of three unit squares.
Define the coordinate values of the corners of the region.
xmin =-1; xmax = 1; ymin =-1; ymax = 1;
Consider a square grid consisting of an odd number Nx=2*nx-1
of grid points in x
direction and an odd number Ny=2*ny-1
of grid points in y
direction.
nx = 6; Nx = 2*nx-1; hx = (xmax-xmin)/(Nx-1); ny = 6; Ny = 2*ny-1; hy = (ymax-ymin)/(Ny-1);
Create an Ny
-by- Nx
symbolic matrix . Its entries u(i,j)
represent the values u(xmin + (j - 1)*hx,ymin + (i - 1)*hy)
of the solution u(x,y)
of the eigenvalue problem .
u = sym('u',[Ny,Nx]);
The boundaries of correspond to the following indices:
The left boundary corresponds to (i = 1:Ny, j = 1)
.
The lower boundary corresponds to (i = 1, j = 1:Nx)
.
The right boundary corresponds to (i = 1:ny, j = Nx)
and (i = ny:Ny, j = nx)
.
The upper boundary corresponds to (i = Ny, j =1:nx)
and (i = ny, j = nx:Nx)
.
u(:,1) = 0; % left boundary u(1,:) = 0; % lower boundary u(1:ny,Nx) = 0; % right boundary, upper part u(ny:Ny,nx) = 0; % right boundary, lower part u(Ny,1:nx) = 0; % upper boundary, left part u(ny,nx:Nx) = 0; % upper boundary, right part
The region with and does not belong to . Set the corresponding matrix entries (i = ny + 1:Ny, j = nx + 1:Nx)
to zero. They play no further role and will be ignored.
u(ny + 1:Ny,nx + 1:Nx) = 0;
The unknowns of the problem are the following matrix entries:
u
u =
The interior points of the region correspond to the indices that contain the unknown values of the problem. Collect these unknowns in a vector vars
.
[I,J] = find(u~=0); vars = u(u~=0);
Associate a symbolic expression (given by the stencil derived in the first part of this example) with each index (that is, with each unknown).
n = length(vars); Lu = sym(zeros(n,1)); for k=1:n i = I(k); j = J(k); Lu(k) = 1/6*u(i+1,j-1) + 2/3*u(i+1,j) + 1/6*u(i+1,j+1) ... + 2/3*u( i ,j-1) - 10/3*u( i ,j) + 2/3*u( i ,j+1) ... + 1/6*u(i-1,j-1) + 2/3*u(i-1,j) + 1/6*u(i-1,j+1); end Lu = Lu/hx^2;
Because this expression is linear in the unknown elements of u
(stored in vars
), you can treat it as a matrix acting on the vector vars
.
S_h = jacobian(Lu, vars);
You can treat S_h
as a matrix approximation of the Laplace operator. Compute its eigenvectors and eigenvalues.
[V,D] = eig(vpa(S_h));
The three maximal eigenvalues are given by the first diagonal elements of D.
[D(1,1), D(2,2), D(3,3)]
ans =
Since for this approximation you used a grid with a small number of points, only the leading digits of the eigenvalues are correct.
The third highest eigenvalue of the Laplace operator on the L-shaped region is known exactly. The exact eigenfunction of the Laplace operator is the function associated with the (exact) eigenvalue . Indeed, using equation (3) above, you can derive a better approximation of the Laplace eigenvalue from the stencil eigenvalue :
mu = D(3,3)
mu =
lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)
lambda =
Plot the eigenfunction associated with the third highest eigenvalue.
v = V(:,3); for k=1:n u(I(k),J(k)) = v(k); end u = double(u); surf(xmin:hx:xmax,ymin:hy:ymax,u'); view(125, 30);
When you use symbolic matrices, increasing the number of grid points drastically is not recommended because symbolic computations are significantly slower than numerical computations with MATLAB double-precision matrices. This part of the example demonstrates how to use sparse double arithmetic which allows to refine the numerical grid. The L-shaped region is set up the same way as before. Instead of denoting the interior points by symbolic unknowns, initialize the grid values u
by ones and define by setting the values of boundary points and exterior points to zero. Instead of defining a symbolic expression for each interior point and computing the stencil as the Jacobian, set up the stencil matrix directly as a sparse matrix.
xmin =-1; xmax = 1; ymin =-1; ymax = 1; nx = 30; Nx = 2*nx-1; hx = (xmax-xmin)/(Nx-1); ny = 30; Ny = 2*ny-1; hy = (ymax-ymin)/(Ny-1); u = ones(Ny,Nx); u(:,1) = 0; % left boundary u(1:ny,Nx) = 0; % right boundary, upper part u(ny:Ny,nx) = 0; % right boundary, lower part u(1,:) = 0; % lower boundary u(Ny,1:nx) = 0; % upper boundary, left part u(ny,nx:Nx) = 0; % upper boundary, right part u(ny + 1:Ny,nx + 1:Nx) = 0; [I,J] = find(u ~= 0); n = length(I); S_h = sparse(n,n); for k=1:n i = I(k); j = J(k); S_h(k,I==i+1 & J==j+1)= 1/6; S_h(k,I==i+1 & J== j )= 2/3; S_h(k,I==i+1 & J==j-1)= 1/6; S_h(k,I== i & J==j+1)= 2/3; S_h(k,I== i & J== j )=-10/3; S_h(k,I== i & J==j-1)= 2/3; S_h(k,I==i-1 & J==j+1)= 1/6; S_h(k,I==i-1 & J== j )= 2/3; S_h(k,I==i-1 & J==j-1)= 1/6; end S_h = S_h./hx^2;
Here, S_h
is the (sparse) stencil matrix. Use eigs
that handles sparse matrices to compute the three largest eigenvalues.
[V,D] = eigs(S_h,3,'la');
The three maximal eigenvalues are the first diagonal elements of D.
[D(1,1), D(2,2), D(3,3)]
ans = 1×3
-9.6493 -15.1742 -19.7006
D(3,3)
approximates the exact eigenvalue . Using the equation (3) above, derive a more accurate approximation of the Laplace eigenvalue from the stencil eigenvalue .
mu = D(3,3)
mu = -19.7006
lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)
lambda = -19.7393
Plot the eigenfunction associated with the third highest eigenvalue.
v = V(:,3); for k=1:n u(I(k),J(k)) = v(k); end surf(xmin:hx:xmax, ymin:hy:ymax,u'); view(125, 30);
Note that the MATLAB membrane
function computes the eigenfunctions of the Laplace operator by different methods.
membrane(3, nx - 1, 8, 8);