Solve eigenvalue problem of the Laplacian by finite element method

On a square domain $\Omega$, consider the following eigenvalue problem

$$ -\Delta u = \lambda u \mbox{ in } \Omega, u=0 \mbox{ on } \partial \Omega \:. $$

The variational formulation for the above eigenvalue problem is to find $u\in H_0^1(\Omega)$ and $\lambda \in R$ such that $$ \int_{\Omega} \nabla u \cdot \nabla v dx = \lambda \int_{\Omega} uvdx \mbox{ for all } v \in H_0^1(\Omega) \:. $$

Below, we show how to solve the eigenvalue problem step by step.

Particular, in this speical case, the first 9 eigenvalue are given by $$ \left\{\frac{\lambda_i}{\pi^2} \right\}_{i=1}^{9} = \{2,5,5,8,10,10,13,13,18\} \:. $$

Objective:

  • Calculate the lower bounds for the leading 9 exact eigenvalues.

  • The Crouzeix-Rarviart (CR) finite element method will be used to provide lower eigenvalue bounds; such a method even works for non-convex domain.

  • The lower and upper eigenvalue bounds for L-shaped domain are given at the end of this file.

Last updated by Xuefeng LIU, Sep. 4, 2017

Step 1 : Mesh generation and FEM space definition

In [9]:
from dolfin import *

N=16; h=1.0/N
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CR", 1)

# define Dirichlet boundary conditions
def bdry(x, on_boundary):  return on_boundary

bc = DirichletBC(V, Constant(0.0), bdry)

Step 2: Variational formulation

In [2]:
# Define basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
b = dot(u, v)*dx

L = v*dx # To feed an argument to assemble_system

# Assemble the stiffness matrix and the mass matrix.
A, _ = assemble_system(a, L, bc) #The assemble_system commands make sure the symmetry of A
B = assemble(b)

# set the diagonal elements of B corresponding to boundary nodes to zero to
# remove spurious eigenvalues.
bc.zero(B)
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.

Step 3: Calculate matrix and solve the matrix eigenvalue problem

In [3]:
# downcast to PETSc matrices
MatA = as_backend_type(A)
MatB = as_backend_type(B)

# Create eigensolver
eigensolver = SLEPcEigenSolver(MatA, MatB)
# Specify the part of the spectrum desired
eigensolver.parameters["spectrum"] = "smallest magnitude"
# Specify the problem type (this can make a big difference)
eigensolver.parameters["problem_type"] = "gen_hermitian"
# Use the shift-and-invert spectral transformation
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
# Specify the shift
eigensolver.parameters["spectral_shift"] = 1.0e-10

# Compute the eigenvalues. Here is where we specify the number of eigenvalues.
nreq = 9
eigensolver.solve(nreq)

Step 4: Approximate eigenvalues obtained by CR FEM

In [4]:
import numpy as np

exact_eigvalues = np.array([2,5,5,8,10,10,13,13,18])*pi**2;

# Extract the leading eigenpair from the smallest eigenvalue.
for k in range(0,9):
    r, c, rx, cx = eigensolver.get_eigenpair(k)
    exact_eig = exact_eigvalues[k]    
    print "The %dth approximate eigenvalue:%8.3f*pi^2 (exact one:%4d*pi^2)"%(k+1, r/(pi**2), np.rint(exact_eig/(pi**2)))
The 1th approximate eigenvalue:   1.998*pi^2 (exact one:   2*pi^2)
The 2th approximate eigenvalue:   4.972*pi^2 (exact one:   5*pi^2)
The 3th approximate eigenvalue:   4.972*pi^2 (exact one:   5*pi^2)
The 4th approximate eigenvalue:   7.966*pi^2 (exact one:   8*pi^2)
The 5th approximate eigenvalue:   9.843*pi^2 (exact one:  10*pi^2)
The 6th approximate eigenvalue:   9.843*pi^2 (exact one:  10*pi^2)
The 7th approximate eigenvalue:  12.869*pi^2 (exact one:  13*pi^2)
The 8th approximate eigenvalue:  12.869*pi^2 (exact one:  13*pi^2)
The 9th approximate eigenvalue:  16.482*pi^2 (exact one:  18*pi^2)

Step 5: Draw the eigenfunction

Please choose eig_index to be from 0 to 8 and check the shape of eigenfunction.

In [7]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
%matplotlib inline

#Please choose the eige_index from 0 to 8.
eig_index=1
r, c, rx, cx = eigensolver.get_eigenpair(eig_index)

nodes = mesh.coordinates()
x = nodes[:,0]
y = nodes[:,1]

u = Function(V)
u.vector()[:] = rx 
z = u.compute_vertex_values(mesh)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(x, y, z, cmap=cm.jet, linewidth=0.2)
plt.show()

6. Lower eigenvalue bounds

Calculate the lower bounds for the leading 9 exact eigenvalues. To give lower bound, we need the following error estimation for projection ($P_h$) to CR FEM space

$$ \| u - P_h u \| \le 0.1893h \| \nabla(u - P_h u) \| $$

where $h$ is the maximum edge length of triangulation of the domain.

Notice

This method works even for non-convex domain.

In [8]:
Ch=0.1893*h;

import numpy as np

exact_eigvalues = np.array([2,5,5,8,10,10,13,13,18])*pi**2;

# Extract the leading eigenpair from the smallest eigenvalue.
for k in range(0,9):
    eig_value, c, rx, cx = eigensolver.get_eigenpair(k)
    exact_eig = exact_eigvalues[k]
    lower_bound = eig_value/(1+Ch*Ch*eig_value)
    print "The lower bound of %d theigenvalue: %10.5f  < %10.5f  (Difference: %10.5f) "%(k+1, lower_bound, exact_eig, exact_eig - lower_bound )
The lower bound of 1 theigenvalue:   19.66379  <   19.73921  (Difference:    0.07542) 
The lower bound of 2 theigenvalue:   48.73813  <   49.34802  (Difference:    0.60989) 
The lower bound of 3 theigenvalue:   48.73813  <   49.34802  (Difference:    0.60989) 
The lower bound of 4 theigenvalue:   77.76226  <   78.95684  (Difference:    1.19458) 
The lower bound of 5 theigenvalue:   95.84619  <   98.69604  (Difference:    2.84985) 
The lower bound of 6 theigenvalue:   95.84619  <   98.69604  (Difference:    2.84985) 
The lower bound of 7 theigenvalue:  124.79299  <  128.30486  (Difference:    3.51187) 
The lower bound of 8 theigenvalue:  124.79299  <  128.30486  (Difference:    3.51187) 
The lower bound of 9 theigenvalue:  159.04894  <  177.65288  (Difference:   18.60394) 

Setp 7. Upper and Lower eigenvalue bounds for L-shaped domain

We apply the method described above to L-shaped domain. Please make sure the L_uniform.xml is uploaded to the same folder.

First, let check the mesh for L-shaped domain. Use L_uniform.xml as inital mesh. By refine this mesh, we have dense uiform mesh.

In [28]:
from dolfin import *

%matplotlib inline

import matplotlib.pyplot as plt
import matplotlib.tri as tri

N=4; h=1.0/N
mesh = Mesh("L_uniform.xml")
for k in range(0,N-1):
    mesh = refine(mesh)
    
nodes = mesh.coordinates();
elements = tri.Triangulation(nodes[:, 0], nodes[:, 1], mesh.cells());
plt.triplot(elements)
Out[28]:
[<matplotlib.lines.Line2D at 0x7fa3f44b9910>,
 <matplotlib.lines.Line2D at 0x7fa3f44b9b90>]

7.1 Lower bound

In [13]:
V = FunctionSpace(mesh, "CR", 1)

# define Dirichlet boundary conditions
def bdry(x, on_boundary):  return on_boundary

bc = DirichletBC(V, Constant(0.0), bdry)

# Define basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
b = dot(u, v)*dx

L = v*dx # To feed an argument to assemble_system

# Assemble the stiffness matrix and the mass matrix.
A, _ = assemble_system(a, L, bc) #The assemble_system commands make sure the symmetry of A
B = assemble(b)

# set the diagonal elements of B corresponding to boundary nodes to zero to
# remove spurious eigenvalues.
bc.zero(B)

# downcast to PETSc matrices
MatA = as_backend_type(A)
MatB = as_backend_type(B)

# Create eigensolver
eigensolver = SLEPcEigenSolver(MatA, MatB)
# Specify the part of the spectrum desired
eigensolver.parameters["spectrum"] = "smallest magnitude"
# Specify the problem type (this can make a big difference)
eigensolver.parameters["problem_type"] = "gen_hermitian"
# Use the shift-and-invert spectral transformation
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
# Specify the shift
eigensolver.parameters["spectral_shift"] = 1.0e-10

# Compute the eigenvalues. Here is where we specify the number of eigenvalues.
nreq = 9
eigensolver.solve(nreq)


import numpy as np

print "\nApproximate eigenvalues by CR method"

# Extract the leading eigenpair from the smallest eigenvalue.
for k in range(0,9):
    r, c, rx, cx = eigensolver.get_eigenpair(k)
    print "The %dth approximate eigenvalue:%8.3f"%(k+1, r)
    

Ch=0.1893*h;

import numpy as np

print "\nLower bound of eigenvalues"

lower_bound_list = []

# Extract the leading eigenpair from the smallest eigenvalue.
for k in range(0,9):
    eig_value, c, rx, cx = eigensolver.get_eigenpair(k)
    lower_bound = eig_value/(1+Ch*Ch*eig_value)
    lower_bound_list.append(lower_bound)
    print "The lower bound of %d theigenvalue: %10.5f"%(k+1, lower_bound)
Approximate eigenvalues by CR method
The 1th approximate eigenvalue:   9.461
The 2th approximate eigenvalue:  15.110
The 3th approximate eigenvalue:  19.655
The 4th approximate eigenvalue:  29.183
The 5th approximate eigenvalue:  31.118
The 6th approximate eigenvalue:  40.214
The 7th approximate eigenvalue:  43.699
The 8th approximate eigenvalue:  48.244
The 9th approximate eigenvalue:  55.485

Lower bound of eigenvalues
The lower bound of 1 theigenvalue:    9.26488
The lower bound of 2 theigenvalue:   14.61512
The lower bound of 3 theigenvalue:   18.82581
The lower bound of 4 theigenvalue:   27.39247
The lower bound of 5 theigenvalue:   29.09059
The lower bound of 6 theigenvalue:   36.89143
The lower bound of 7 theigenvalue:   39.80325
The lower bound of 8 theigenvalue:   43.53951
The lower bound of 9 theigenvalue:   49.35240

7.2 Upper bound

In [12]:
V = FunctionSpace(mesh, "CG", 1)

# define Dirichlet boundary conditions
def bdry(x, on_boundary):  return on_boundary
bc = DirichletBC(V, Constant(0.0), bdry)

# Define basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
b = dot(u, v)*dx

L = v*dx # To feed an argument to assemble_system

# Assemble the stiffness matrix and the mass matrix.
A, _ = assemble_system(a, L, bc) #The assemble_system commands make sure the symmetry of A
B = assemble(b)

# set the diagonal elements of B corresponding to boundary nodes to zero to
# remove spurious eigenvalues.
bc.zero(B)

# downcast to PETSc matrices
MatA = as_backend_type(A)
MatB = as_backend_type(B)

# Create eigensolver
eigensolver = SLEPcEigenSolver(MatA, MatB)
# Specify the part of the spectrum desired
eigensolver.parameters["spectrum"] = "smallest magnitude"
# Specify the problem type (this can make a big difference)
eigensolver.parameters["problem_type"] = "gen_hermitian"
# Use the shift-and-invert spectral transformation
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
# Specify the shift
eigensolver.parameters["spectral_shift"] = 1.0e-10

# Compute the eigenvalues. Here is where we specify the number of eigenvalues.
nreq = 9
eigensolver.solve(nreq)

import numpy as np

print "\nApproximate eigenvalues by Lagrange method, which are also the exact upper bounds"

upper_bound_list = []

# Extract the leading eigenpair from the smallest eigenvalue.
for k in range(0,9):
    r, c, rx, cx = eigensolver.get_eigenpair(k)
    upper_bound_list.append(r)
    print "The upper bound of %dth eigenvalue:%8.3f"%(k+1, r)
Approximate eigenvalues by Lagrange method, which are also the exact upper bounds
The upper bound of 1th eigenvalue:   9.923
The upper bound of 2th eigenvalue:  15.545
The upper bound of 3th eigenvalue:  20.422
The upper bound of 4th eigenvalue:  30.898
The upper bound of 5th eigenvalue:  33.880
The upper bound of 6th eigenvalue:  44.122
The upper bound of 7th eigenvalue:  47.609
The upper bound of 8th eigenvalue:  53.088
The upper bound of 9th eigenvalue:  62.525

7.3 Draw the lower bound and upper bound in graph

In [26]:
%matplotlib inline
import matplotlib.pyplot as plt

upper_, = plt.plot( range(0,9), upper_bound_list, 'bo-',label="Upper bound")
lower_, = plt.plot( range(0,9), lower_bound_list, 'ro-',label="Lower bound")

plt.legend(handles=[upper_, lower_ ])

plt.grid()
plt.title("Lower and upper eigenvalue bound (L-shaped domain)")
plt.show()
In [ ]: