Eigenvalues of the Laplace Operator

This example shows how to solve the eigenvalue problem of the Laplace operator on an L-shaped region.

Membrane Problem

Consider a membrane that is fixed at the boundary Ω of a region Ω in the plane. Its displacement u(x,y) is described by the eigenvalue problem Δu=λu, where Δu=uxx+uyy is the Laplace operator and λ is a scalar parameter. The boundary condition is u(x,y)=0 for all (x,y)Ω.

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

Nine-Point Finite Difference Approximation

The simplest approach to the eigenvalue problem is to approximate the Laplacian Δu by a finite difference approximation (a stencil) on a square grid of points with distances hx in x direction and distances hy in y direction. In this example, approximate Δu with a sum S_h of nine regular grid points around the midpoint (x,y). The unknowns are the weights a-1-1,,a11.

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 x and y. Since S_h represents uxx+uyy, the coefficients of all other derivatives of u must be zero. Extract the coefficients by replacing all derivatives of u, except uxx and uyy, by 0. Replace uxx and uyy 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], ...

ans = 

(z1hy2-2zz1hx2-2z4z-2hx2-2hy21hy2-2zz1hy2-2zz)[z, 1/hy^2 - 2*z, z; 1/hx^2 - 2*z, 4*z - 2/hx^2 - 2/hy^2, 1/hy^2 - 2*z; z, 1/hy^2 - 2*z, z]

parameters = zz

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 = (000)[sym(0), sym(0), sym(0)]

The expression C(5) is the Laplacian of u.

ans = 

2x2 u(x,y)+2y2 u(x,y)diff(u(x, y), x, 2) + diff(u(x, y), y, 2)

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

Terms Containing Fourth and Higher Order Derivatives

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]);
ans(x, y) = 

d4x4 u(x,y)+2d2y2 2x2 u(x,y)+d4y4 u(x,y)d*diff(u(x, y), x, 4) + 2*d*diff(diff(u(x, y), x, 2), y, 2) + d*diff(u(x, y), y, 4)

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 = 

hx212=dhx^2/12 == d

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 = hx2hy2z=2dhx^2*hy^2*z == 2*d
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 = 

hy212=dhy^2/12 == d

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) = 0sym(0)

Now, consider the third-order terms in hx, hy.

ans = 0sym(0)

Since no such terms exist in the expansion of the stencil, the term O(h3) in (1) is in fact of order O(h4). Consider the fourth-order terms of the stencil.

ans = 

(1360hhhh6x6 u(x,y)+52y2 4x4 u(x,y)+54y4 2x2 u(x,y)+6y6 u(x,y))[sym(1/360), h, h, h, h, diff(u(x, y), x, 6) + 5*diff(diff(u(x, y), x, 4), y, 2) + 5*diff(diff(u(x, y), x, 2), y, 4) + diff(u(x, y), y, 6)]

Check if these terms can be identified with yet another power of the Laplace operator. However, comparison with

ans(x, y) = 

6x6 u(x,y)+32y2 4x4 u(x,y)+34y4 2x2 u(x,y)+6y6 u(x,y)diff(u(x, y), x, 6) + 3*diff(diff(u(x, y), x, 4), y, 2) + 3*diff(diff(u(x, y), x, 2), y, 4) + diff(u(x, y), y, 6)

shows that the expressions of order O(h4) 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 Δu=λu. Add some multiple of Δ2u=λ2u to the eigenvalue equation.


The left side of this equation is well approximated by the stencil Sh. Thus, using (2), a numerical eigenvalue μ of the stencil satisfying Shu=μu 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 h0.


Using Symbolic Matrices to Solve the Eigenvalue Problem

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 u. 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=λu.

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 0<x1 and 0<y1 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 = 

(000000000000u2,2u2,3u2,4u2,5u2,6u2,7u2,8u2,9u2,1000u3,2u3,3u3,4u3,5u3,6u3,7u3,8u3,9u3,1000u4,2u4,3u4,4u4,5u4,6u4,7u4,8u4,9u4,1000u5,2u5,3u5,4u5,5u5,6u5,7u5,8u5,9u5,1000u6,2u6,3u6,4u6,50000000u7,2u7,3u7,4u7,50000000u8,2u8,3u8,4u8,50000000u9,2u9,3u9,4u9,50000000u10,2u10,3u10,4u10,500000000000000000)[sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), u2_2, u2_3, u2_4, u2_5, u2_6, u2_7, u2_8, u2_9, u2_10, sym(0); sym(0), u3_2, u3_3, u3_4, u3_5, u3_6, u3_7, u3_8, u3_9, u3_10, sym(0); sym(0), u4_2, u4_3, u4_4, u4_5, u4_6, u4_7, u4_8, u4_9, u4_10, sym(0); sym(0), u5_2, u5_3, u5_4, u5_5, u5_6, u5_7, u5_8, u5_9, u5_10, sym(0); sym(0), u6_2, u6_3, u6_4, u6_5, sym(0), sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), u7_2, u7_3, u7_4, u7_5, sym(0), sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), u8_2, u8_3, u8_4, u8_5, sym(0), sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), u9_2, u9_3, u9_4, u9_5, sym(0), sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), u10_2, u10_3, u10_4, u10_5, sym(0), sym(0), sym(0), sym(0), sym(0), sym(0); sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0), sym(0)]

The interior points of the region ΩΩ correspond to the indices (i,j) that contain the unknown values u(i,j) 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); 
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 = (-9.5214641572625960021345709535953-14.431096242107969492574666743957-18.490392088545609858994660377955)[-vpa('9.5214641572625960021345709535953'), -vpa('14.431096242107969492574666743957'), -vpa('18.490392088545609858994660377955')]

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 u(x,y)=sin(πx)sin(πy) associated with the (exact) eigenvalue -2π2=-19.7392.... Indeed, using equation (3) above, you can derive a better approximation of the Laplace eigenvalue λ from the stencil eigenvalue μ:

mu = D(3,3)
mu = -18.490392088545609858994660377955-vpa('18.490392088545609858994660377955')
lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)
lambda = -19.796765119155672176257649532142-vpa('19.796765119155672176257649532142')

Plot the eigenfunction associated with the third highest eigenvalue.

v = V(:,3);
for k=1:n
  u(I(k),J(k)) = v(k);
u = double(u);
view(125, 30);

Using Double-Precision Matrices to Solve the Eigenvalue Problem

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;
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 -2π2=-19.7392088.... 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);
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);