Hi,
The eigenfunctions for a particle in a 1D infinite well are the usual sin functions of the form
$U_1 = \sin(\pi x)$
$U_2 = \sin(2 \pi x)$
$U_3 = \sin(3 \pi x)$
etc.
However, when I solve this problem numerically FEniCS will sometimes produce eigenfunctions with the wrong sign
The code below, for example, produces the eigenfunctions
$U_1 = -\sin(\pi x)$
$U_2 = \sin(2 \pi x)$
even though it reports positive eigenvalues. Is there a way to insist on eigenfunctions that produce positive eigenvalues?
from dolfin import *
#define mesh and function space
mesh = IntervalMesh(200,-5,5)
V = FunctionSpace(mesh, 'Lagrange', 1)
#build essential boundary conditions
def u0_boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V,Constant(0.0) , u0_boundary)
#define functions
u = TrialFunction(V)
v = TestFunction(V)
#define problem
a = (inner(grad(u), grad(v)) \
+ Constant('0.0')*u*v)*dx
m = u*v*dx
A = PETScMatrix()
M = PETScMatrix()
_ = PETScVector()
L = Constant(0.)*v*dx
assemble_system(a, L, bc, A_tensor=A, b_tensor=_)
assemble_system(m, L, A_tensor=M, b_tensor=_)
#create eigensolver
eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.parameters['tolerance'] = 1.e-15
#solve for eigenvalues
neig = 2
eigensolver.solve(neig)
u = Function(V)
for i in range(neig):
#extract next eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
print 'eigenvalue:', r
#assign eigenvector to function
u.vector()[:] = rx
plot(u,interactive=True)
Is there a way to insist only one form is calculated?