High-precision eigenvalue bounds of the Laplacian

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 $\lambda_i$ are given by $$ \left\{\frac{\lambda_i}{\pi^2} \right\}_{i=1}^{9} = \{2,5,5,8,10,10,13,13,18\} \:. $$

Question:

Show the high-precision lower and upper bounds for the leading 9 exact eigenvalues.

Last updated by Xuefeng LIU, Sep. 4, 2017

Step 1 : Mesh generation and FEM space definition

In [1]:
from dolfin import *

degree=2

N=16*2; h=1.0/N
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CG", degree)

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

bc = DirichletBC(V, Constant(0.0), bdry)
Calling FFC just-in-time (JIT) compiler, this may take some time.

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 = 5
eigensolver.solve(nreq)

Step 4: Extract the eigenvalue and eigenfunction

In [4]:
#-------------------------------- solve the eigenvalue problem of matrices --------------------

eig_vec_list=[]
eig_list=[]


# Check the number of eigenvalues that converged.
# Extract the eigenvalues (ignore the imaginary part) and compare with the exact values
nconv = eigensolver.get_number_converged()
print "Number of eigenvalues successfully computed: ", nconv
nconv = 5
eigenvalues = []
for i in range(nconv):
   r, c, rx, cx = eigensolver.get_eigenpair(i) # real and complex part of eigenvalue
   eigenvalues.append(r)
   eigenfun = Function(V)
   eigenfun.vector()[:] = rx
   eig_vec_list.append(eigenfun)
   eig_list.append(r)
   print "eignvalue", r

#-------------------------------- Apply Lehmann-Goerish theorem to obtain high-precision bounds -------------------
Number of eigenvalues successfully computed:  5
eignvalue 19.7392265967
eignvalue 49.3481880371
eignvalue 49.3483252128
eignvalue 78.957967741
eignvalue 98.6976497114

Step 5: Draw the eigenfunction

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

In [5]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

#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()
/usr/lib/python2.7/dist-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Lehmann-Goerisch's method:

In [7]:
import  numpy as np

eig_index=0
r, c, rx, cx = eigensolver.get_eigenpair(eig_index)

eig_f = Function(V)
eig_f.vector()[:] = rx 


def domain_bd(x, on_boundary):
    return near(x[1], 1.0) or near(x[1], 0.0) or near(x[0], 1.0) or near(x[0], 0.0)

# Define function G such that G \cdot n = eig_f
deg = degree

RT = FiniteElement("RT", mesh.ufl_cell(), deg)
DG = FiniteElement("DG", mesh.ufl_cell(), deg-1)

W = FunctionSpace( mesh , MixedElement ([RT, DG]) )

(u,rho) = TrialFunctions(W)
(v,phi) = TestFunctions(W)

a = inner(u,v)*dx + div(v)*rho*dx + div(u)*phi*dx
L = -phi*eig_f*dx

# Compute solution
sol = Function(W)
solve(a == L, sol)

(w,psi) = sol.split()

#Save w into pvd file. (opened by Paraview)
File("sol_vec.pvd") << w

#print "div norm", np.power(1.0*norm(w, "Hdiv"),2)
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.
In [14]:
#Take a upper bound for the first eigenvalue and lower bound for the second one.
rho = 49.3

A0 = assemble( inner(grad(eig_f),grad(eig_f))*dx )
A1 = assemble( eig_f*eig_f*dx )
A2 = assemble( inner(w,w)*dx  )
AL = A0 - rho*A1;
BL = A0 - 2*rho*A1 + rho*rho*A2;
mu = AL/BL;
lower_bound = rho - rho/(1-mu)

eig_v = eig_list[0]

print "Upper bound is", eig_v
print "Lower bound is", lower_bound
print "Gap is ", eig_v - lower_bound
print "Error of lower bound is ", 2*np.pi*np.pi - lower_bound
Upper bound is 19.7392265967
Lower bound is 19.7392172051
Gap is  9.39160794644e-06
Error of lower bound is  -8.40294629967e-06
In [ ]: