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 $$\partial \Omega $$ of a region $$\Omega $$ in the plane. Its displacement $$u(x,y)$$ is described by the eigenvalue problem $$\Delta u=\lambda u$$, where $$\Delta u={u}_{xx}+{u}_{yy}$$ is the Laplace operator and $$\lambda $$ is a scalar parameter. The boundary condition is $$u(x,y)=0$$ for all $$(x,y)\in \partial \Omega $$.

The Laplace operator is self-adjoint and negative definite, that is, only real negative eigenvalues $$\lambda $$ exist. There is a maximal (negative) discrete eigenvalue, the corresponding eigenfunction $$u$$ is called the *ground state*. In this example, $$\Omega $$ 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 $$\Delta 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 $$\Delta u$$ with a sum `S_h`

of nine regular grid points around the midpoint $$(x,y)$$. The unknowns are the weights $${a}_{-1-1},\dots ,{a}_{11}$$.

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 $${u}_{xx}+{u}_{yy}$$, the coefficients of all other derivatives of `u`

must be zero. Extract the coefficients by replacing all derivatives of `u`

, except $${u}_{xx}$$ and $${u}_{yy}$$, by 0. Replace $${u}_{xx}$$ and $${u}_{yy}$$ 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 =$$\left(\begin{array}{ccc}z& \frac{1}{{\mathrm{hy}}^{2}}-2\hspace{0.17em}z& z\\ \frac{1}{{\mathrm{hx}}^{2}}-2\hspace{0.17em}z& 4\hspace{0.17em}z-\frac{2}{{\mathrm{hx}}^{2}}-\frac{2}{{\mathrm{hy}}^{2}}& \frac{1}{{\mathrm{hy}}^{2}}-2\hspace{0.17em}z\\ z& \frac{1}{{\mathrm{hy}}^{2}}-2\hspace{0.17em}z& z\end{array}\right)$$

parameters

`parameters = $$z$$`

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 = $$\left(\begin{array}{ccc}0& 0& 0\end{array}\right)$$`

The expression `C(5)`

is the Laplacian of `u`

.

C(5)

ans =$$\frac{{\partial}^{2}}{\partial {x}^{2}}\mathrm{}u\left(x,y\right)+\frac{{\partial}^{2}}{\partial {y}^{2}}\mathrm{}u\left(x,y\right)$$

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) =$$d\hspace{0.17em}\frac{{\partial}^{4}}{\partial {x}^{4}}\mathrm{}u\left(x,y\right)+2\hspace{0.17em}d\hspace{0.17em}\frac{{\partial}^{2}}{\partial {y}^{2}}\mathrm{}\frac{{\partial}^{2}}{\partial {x}^{2}}\mathrm{}u\left(x,y\right)+d\hspace{0.17em}\frac{{\partial}^{4}}{\partial {y}^{4}}\mathrm{}u\left(x,y\right)$$

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 =$$\frac{{\mathrm{hx}}^{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 = $${\mathrm{hx}}^{2}\hspace{0.17em}{\mathrm{hy}}^{2}\hspace{0.17em}z=2\hspace{0.17em}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 =$$\frac{{\mathrm{hy}}^{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`

.

$${S}_{h}=\Delta u+\frac{{h}^{2}}{12}\phantom{\rule{0.2777777777777778em}{0ex}}{\Delta}^{2}u+O({h}^{3}).\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(1)$$

```
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) = $$0$$`

Now, consider the third-order terms in `hx`

, `hy`

.

simplify(C(2))

`ans = $$0$$`

Since no such terms exist in the expansion of the stencil, the term $$O({h}^{3})$$ in (1) is in fact of order $$O({h}^{4})$$. Consider the fourth-order terms of the stencil.

factor(simplify(C(1)))

ans =$$\left(\begin{array}{cccccc}\frac{1}{360}& h& h& h& h& \frac{{\partial}^{6}}{\partial {x}^{6}}\mathrm{}u\left(x,y\right)+5\hspace{0.17em}\frac{{\partial}^{2}}{\partial {y}^{2}}\mathrm{}\frac{{\partial}^{4}}{\partial {x}^{4}}\mathrm{}u\left(x,y\right)+5\hspace{0.17em}\frac{{\partial}^{4}}{\partial {y}^{4}}\mathrm{}\frac{{\partial}^{2}}{\partial {x}^{2}}\mathrm{}u\left(x,y\right)+\frac{{\partial}^{6}}{\partial {y}^{6}}\mathrm{}u\left(x,y\right)\end{array}\right)$$

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) =$$\frac{{\partial}^{6}}{\partial {x}^{6}}\mathrm{}u\left(x,y\right)+3\hspace{0.17em}\frac{{\partial}^{2}}{\partial {y}^{2}}\mathrm{}\frac{{\partial}^{4}}{\partial {x}^{4}}\mathrm{}u\left(x,y\right)+3\hspace{0.17em}\frac{{\partial}^{4}}{\partial {y}^{4}}\mathrm{}\frac{{\partial}^{2}}{\partial {x}^{2}}\mathrm{}u\left(x,y\right)+\frac{{\partial}^{6}}{\partial {y}^{6}}\mathrm{}u\left(x,y\right)$$

shows that the expressions of order $$O({h}^{4})$$ 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:

$${S}_{h}=\Delta u+\frac{{h}^{2}}{12}\phantom{\rule{0.2777777777777778em}{0ex}}{\Delta}^{2}u+O({h}^{4}).\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(2)$$

Use this expansion for the numeric approach to the eigenvalue problem $$\Delta u=\lambda u$$. Add some multiple of $${\Delta}^{2}u={\lambda}^{2}u$$ to the eigenvalue equation.

$$\Delta u+\frac{{h}^{2}}{12}\phantom{\rule{0.2777777777777778em}{0ex}}{\Delta}^{2}u=(\lambda +\frac{{h}^{2}}{12}\phantom{\rule{0.2777777777777778em}{0ex}}{\lambda}^{2})\phantom{\rule{0.2777777777777778em}{0ex}}u.$$

The left side of this equation is well approximated by the stencil $${S}_{h}$$. Thus, using (2), a numerical eigenvalue $$\mu $$ of the stencil satisfying $${S}_{h}u=\mu u$$ must be an approximation of the eigenvalue $$\lambda $$ of the Laplace operator with

$$\mu =\lambda +\frac{{h}^{2}}{12}\phantom{\rule{0.2777777777777778em}{0ex}}{\lambda}^{2}+O({h}^{4}).$$

For given $$\mu $$, solve for $$\lambda $$ to obtain a better approximation of the Laplace eigenvalue. Note that in the solution of the quadratic equation in $$\lambda $$ the correct sign of the square root is given by the requirement that $$\lambda \to \mu $$ for $$h\to 0$$.

$$\lambda =\frac{6}{{h}^{2}}\phantom{\rule{0.2777777777777778em}{0ex}}(\sqrt{1+\frac{\mu {h}^{2}}{3}}-1)=\frac{2\mu}{\sqrt{1+\frac{\mu {h}^{2}}{3}}+1}.\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(3)$$

Consider an L-shaped region $$\Omega $$ consisting of three unit squares.

$$\Omega =\{(x,y);\phantom{\rule{0.2777777777777778em}{0ex}}-1\le x\le 0,\phantom{\rule{0.2777777777777778em}{0ex}}-1\le y\le 0\}\cup \{(x,y);\phantom{\rule{0.2777777777777778em}{0ex}}0\le x\le 1,\phantom{\rule{0.2777777777777778em}{0ex}}-1\le y\le 0\}\cup \{(x,y);\phantom{\rule{0.2777777777777778em}{0ex}}-1\le x\le 0,\phantom{\rule{0.2777777777777778em}{0ex}}0\le y\le 1\}.$$

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 $$\Delta u=\lambda u$$.

`u = sym('u',[Ny,Nx]);`

The boundaries of $$\Omega $$ 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<x\le 1$$ and $$0<y\le 1$$ does not belong to $$\Omega $$. 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 =$$\left(\begin{array}{ccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& {u}_{2,2}& {u}_{2,3}& {u}_{2,4}& {u}_{2,5}& {u}_{2,6}& {u}_{2,7}& {u}_{2,8}& {u}_{2,9}& {u}_{2,10}& 0\\ 0& {u}_{3,2}& {u}_{3,3}& {u}_{3,4}& {u}_{3,5}& {u}_{3,6}& {u}_{3,7}& {u}_{3,8}& {u}_{3,9}& {u}_{3,10}& 0\\ 0& {u}_{4,2}& {u}_{4,3}& {u}_{4,4}& {u}_{4,5}& {u}_{4,6}& {u}_{4,7}& {u}_{4,8}& {u}_{4,9}& {u}_{4,10}& 0\\ 0& {u}_{5,2}& {u}_{5,3}& {u}_{5,4}& {u}_{5,5}& {u}_{5,6}& {u}_{5,7}& {u}_{5,8}& {u}_{5,9}& {u}_{5,10}& 0\\ 0& {u}_{6,2}& {u}_{6,3}& {u}_{6,4}& {u}_{6,5}& 0& 0& 0& 0& 0& 0\\ 0& {u}_{7,2}& {u}_{7,3}& {u}_{7,4}& {u}_{7,5}& 0& 0& 0& 0& 0& 0\\ 0& {u}_{8,2}& {u}_{8,3}& {u}_{8,4}& {u}_{8,5}& 0& 0& 0& 0& 0& 0\\ 0& {u}_{9,2}& {u}_{9,3}& {u}_{9,4}& {u}_{9,5}& 0& 0& 0& 0& 0& 0\\ 0& {u}_{10,2}& {u}_{10,3}& {u}_{10,4}& {u}_{10,5}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)$$

The interior points of the region $$\Omega \setminus \partial \Omega $$ 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); 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 = $$\left(\begin{array}{ccc}-9.5214641572625960021345709535953& -14.431096242107969492574666743957& -18.490392088545609858994660377955\end{array}\right)$$`

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 $$\Omega $$ is known exactly. The exact eigenfunction of the Laplace operator is the function $$u(x,y)=\mathrm{sin}(\pi x)\phantom{\rule{0.16666666666666666em}{0ex}}\mathrm{sin}(\pi y)$$ associated with the (exact) eigenvalue $$-2{\pi}^{2}=-19.7392...$$. Indeed, using equation (3) above, you can derive a better approximation of the Laplace eigenvalue $$\lambda $$ from the stencil eigenvalue $$\mu $$:

mu = D(3,3)

`mu = $$-18.490392088545609858994660377955$$`

lambda = 2*mu / (sqrt(1 + mu*hx^2/3) + 1)

`lambda = $$-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); 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 $$\Omega $$ 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 $$\Omega $$ 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 $$-2{\pi}^{2}=-19.7392088...$$. Using the equation (3) above, derive a more accurate approximation of the Laplace eigenvalue $$\lambda $$ from the stencil eigenvalue $$\mu $$.

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