Setting Boundary Condition on Eigenvalue Problem.

Question

Hi,

Below is code to plot the first three eigenvalues in an infinite potential well
(adapted from http://en.wikiversity.org/wiki/User:Tclamb/FEniCS#Quantum_Harmonic_Oscillator )

While the boundary condition (solution = 0 at the ends) is defined, it doesn't seem to be implemented. Indeed, when the eigenfunctions are plotted by the code, not all of them obey the boundary condition.

from dolfin import *

#define mesh and function space
mesh = IntervalMesh(100, -10, 10)
V = FunctionSpace(mesh, 'Lagrange', 3)

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

#apply essential boundary conditions
bc = DirichletBC(V, 0, boundary)

#define functions
u = TrialFunction(V)
v = TestFunction(V)

Pot = Expression('0.0')

#define problem
a = (inner(grad(u), grad(v)) \
     + Pot*u*v)*dx
m = u*v*dx

#assemble stiffness matrix
A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)

#create eigensolver
eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.parameters['solver']   = 'lapack'
eigensolver.parameters['tolerance'] = 1.e-15

#solve for eigenvalues
eigensolver.solve()

u = Function(V)
for i in range(0,3):
    #extract next eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    print 'eigenvalue:', r

    #assign eigenvector to function
    u.vector()[:] = rx

    #plot eigenfunction
    plot(u)
    interactive()

Answer 1

Hi, you have created object for Dirichlet bcs but haven't really used it.

#assemble stiffness matrix
A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)
bc.apply(A)          # apply the boundary conditions
bc.apply(M)

Comments

Thanks for this.

---

P.S. It looks like the solver would fail for higher eigenvalues.

If I only apply the boundary conditions to A, however, it works fine.