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 and upper bounds for the leading 9 exact eigenvalues.

Last updated by Xuefeng LIU, July. 13, 2017

Step 1 : Mesh generation and FEM space definition

In [1]:
from dolfin import *

N=16; h=1.0/N
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CG", 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 [4]:
# 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: Extract the eigenvalue and eigenfunction

In [5]:
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:   2.019*pi^2 (exact one:   2*pi^2)
The 2th approximate eigenvalue:   5.083*pi^2 (exact one:   5*pi^2)
The 3th approximate eigenvalue:   5.130*pi^2 (exact one:   5*pi^2)
The 4th approximate eigenvalue:   8.305*pi^2 (exact one:   8*pi^2)
The 5th approximate eigenvalue:  10.381*pi^2 (exact one:  10*pi^2)
The 6th approximate eigenvalue:  10.390*pi^2 (exact one:  10*pi^2)
The 7th approximate eigenvalue:  13.572*pi^2 (exact one:  13*pi^2)
The 8th approximate eigenvalue:  13.983*pi^2 (exact one:  13*pi^2)
The 9th approximate eigenvalue:  18.042*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 [8]:
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=0
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()

Step 6. Lower eigenvalue bounds

Calculate the lower bounds for the leading 9 exact eigenvalues.

Notice

This method only holds for convex domain with homogeneous Dirichlet boundary condition.

In [9]:
Ch=0.493*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.55969  <   19.73921  (Difference:    0.17952) 
The lower bound of 2 theigenvalue:   47.88567  <   49.34802  (Difference:    1.46236) 
The lower bound of 3 theigenvalue:   48.31052  <   49.34802  (Difference:    1.03750) 
The lower bound of 4 theigenvalue:   76.05259  <   78.95684  (Difference:    2.90424) 
The lower bound of 5 theigenvalue:   93.37696  <   98.69604  (Difference:    5.31908) 
The lower bound of 6 theigenvalue:   93.44742  <   98.69604  (Difference:    5.24862) 
The lower bound of 7 theigenvalue:  118.83436  <  128.30486  (Difference:    9.47050) 
The lower bound of 8 theigenvalue:  122.01551  <  128.30486  (Difference:    6.28934) 
The lower bound of 9 theigenvalue:  152.31428  <  177.65288  (Difference:   25.33860) 

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

In [11]:
from dolfin import *

%matplotlib inline

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

mesh = Mesh("L_uniform.xml")
mesh = refine(mesh)
nodes = mesh.coordinates();
elements = tri.Triangulation(nodes[:, 0], nodes[:, 1], mesh.cells());
plt.triplot(elements)
Out[11]:
[<matplotlib.lines.Line2D at 0x7fc0febc9110>,
 <matplotlib.lines.Line2D at 0x7fc0febc92d0>]
In [12]:
N=16; h=1.0/N

mesh = Mesh("L_uniform.xml")
for k in range(0,N-1):
    mesh = refine(mesh)

V = FunctionSpace(mesh, "CR", 1)

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

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

ExceptionTraceback (most recent call last)
<ipython-input-12-ead6ba98bbcd> in <module>()
      3 mesh = Mesh("L_uniform.xml")
      4 for k in range(0,N-1):
----> 5     mesh = refine(mesh)
      6 
      7 V = FunctionSpace(mesh, "CR", 1)

/usr/lib/python2.7/dist-packages/dolfin/mesh/refinement.pyc in refine(mesh, cell_markers, redistribute)
     75     # Call C++ refinement
     76     if cell_markers is None:
---> 77         cpp.refine(refined_mesh, mesh, redistribute)
     78     else:
     79         cpp.refine(refined_mesh, mesh, cell_markers, redistribute)

/usr/lib/python2.7/dist-packages/dolfin/cpp/mesh.pyc in refine(*args)
   6911 
   6912     """
-> 6913     return _mesh.refine(*args)
   6914 
   6915 def p_refine(*args):

Exception: std::bad_alloc
In [ ]: