solutions to eigenvalue problem do not solve boundary conditions

Question

Below is code that solves an eigenvalue problem. Unfortunately, the first couple of solutions do not satisfy the boundary conditions. Is there a way to prevent those from occurring? Maybe I'm not imposing the boundary conditions correctly?

from dolfin import *
import numpy as np

# Create mesh and define function space
nx   = 20
L    = 4.0
mesh = IntervalMesh(nx,0.0,L)

# Order of method
p = 0

# Function Spaces
V0 = FunctionSpace(mesh, 'CG', p+2)

zeroVelocity = Constant(0.0)
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def right(x, on_boundary): return x[0] > L - DOLFIN_EPS
bc0 = DirichletBC(V0, zeroVelocity, left)
bc1 = DirichletBC(V0, zeroVelocity, right)
bcs = [bc0, bc1]

# Parameters
beta = Constant(1e0)
k2   = Constant((2.*np.pi/L)**2)

# Define variational problem
p = TrialFunction(V0)
r = TestFunction(V0)
a = -beta*(r*p)*dx
b = (r.dx(0)*p.dx(0) + k2*r*p)*dx

# Assemble Matrices
A = PETScMatrix()
assemble(a, tensor=A)
[bc.apply(A) for bc in bcs]

B = PETScMatrix()
assemble(b, tensor=B)
[bc.apply(B) for bc in bcs]

# Create eigensolver
eigensolver = SLEPcEigenSolver(A, B)

eigensolver.parameters["spectrum"] = "largest magnitude"

# Compute all eigenvalues of A x = \lambda B x
print("Computing eigenvalues. This can take a minute.")
num = 6
eigensolver.solve(num)

print -1e0/((2.*np.pi/L)**2 + (np.arange(1,5)*np.pi/4.0)**2)

psi = Function(V0)
psi_vec = psi.vector()

# Extract largest (first) eigenpair
print("Before get_eigenpair")
for i in range(num):
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    print("Eigenvalue:", i, r, c)

    plot(Function(V0, rx), title='%d real' % i)

interactive()

Comments

You are more likely to get an answer if you format your code correctly.

---

Thanks for replying.

I tried to follow the format of the eigenvalue problem demo as much as possible. Looks pretty similar to me but I am just a novice and probably messed something up, just not sure what.

If you can point me towards what part of the formatting that is bad that yields incorrect results I am happy to look into it. At the moment I am at a loss as to what part is badly formatted.

Maybe my presentation of the code is done badly as well, I am not sure.

---

I don't mean the code itself, just the way it displays for us to read :-)
If you edit your question, select your code and click 'code sample' (the parenthesis { } button) you will see what I mean.

---

Ah, thanks!

I added this line and tried it again. I am asking for 6 eigenvalues and it says that only 5 have converged. So it seems that the first one is not a problem because it is not a converged eigenvalue.

eigensolver.parameters['verbose'] = True

This leads to another related question, is it possible to ignore the eigenvalues that have not converged?

Thanks again.

---

If you for some reason only expect negative eigenvalues, you can just discard solutions with positive ones (or eigenvalues close to 1.0). Not the best solution, but in my case it helped.

There are some related questions you can check, eg.
http://fenicsproject.org/qa/7956/symmetry-and-dirichlet-eigenvalue-problem?show=7956#q7956
http://fenicsproject.org/qa/6707/assemble-system-for-eigenvalue-problem?show=6707#q6707

---

Thanks for the suggestion and the links. I read them over and they are interesting but maybe not quite what I was looking for.

In this case the spurious eigenvalue is positive, and equal to 1, so easy to spot but in other problems that is not the case.

I thought that there were only 5 converges eigenvalues but I went to the following link

https://answers.launchpad.net/dolfin/+question/210106

and found out how to compute the number of converged eigenvalues and only show those.

It seems that the first mode with eigenvalue 1, which does not satisfy the correct boundary conditions, is considered to be converged.

The fact that the boundary conditions are not satisfied should indicate that it is not a valid solution. I guess I'm back to not knowing why the boundary conditions are not being correctly enforced.

---

A friend pointed out that the link that I posted actually explained what to do. In case other people might benefit from this I am including the revised code below.

Thanks for all the help everyone. I'm happy to say that now I only get converged modes that satisfy the correct Dirichlet boundary conditions.

from dolfin import *
import numpy as np

# Create mesh and define function space
nx   = 200
L    = 4.0
mesh = IntervalMesh(nx,0.0,L)

# Order of method
p = 0

# Function Spaces
V0 = FunctionSpace(mesh, 'CG', p+2)

# Parameters
beta = Constant(1e0)
k2   = Constant((2.*np.pi/L)**2)

# Define variational problem
p = TrialFunction(V0)
r = TestFunction(V0)
a = -beta*(r*p)*dx
b = (r.dx(0)*p.dx(0) + k2*r*p)*dx
tm = r*dx

# boundary
def bdry(x, on_boundary):
    return on_boundary

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

# Assemble Matrices
A,_ = assemble_system(a, tm, bc)
B,_ = assemble_system(b, tm, bc)
bc.zero(B)

A = as_backend_type(A)
B = as_backend_type(B)

# Create eigensolver
eigensolver = SLEPcEigenSolver(A, B)

eigensolver.parameters["spectrum"] = "largest magnitude"
eigensolver.parameters["problem_type"] = "gen_hermitian"
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectral_shift"] = 1.0e-10
eigensolver.parameters["verbose"] = False

# Compute all eigenvalues of A x = \lambda B x
print("Computing eigenvalues. This can take a minute.")
num = 6
eigensolver.solve(num)

nconv = eigensolver.get_number_converged()
print "Number of eigenvalues successfully computed: ", nconv

print -1e0/((2.*np.pi/L)**2 + (np.arange(1,5)*np.pi/4.0)**2)

psi = Function(V0)
psi_vec = psi.vector()

# Extract largest (first) eigenpair
for i in range(nconv):
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    print("Eigenvalue:", i, r, c)

    plot(Function(V0, rx), title='%d real' % i)

interactive()