Solving 2D Schrodinger Equation

Question

Hi,

I am trying to solve a basic 2D time-independent Schrodinger equation using FEniCS

http://imgur.com/KHTaTie

where A is a constant. I have found the following example, which has been very helpful:
http://en.wikiversity.org/wiki/User:Tclamb/FEniCS#Quantum_Harmonic_Oscillator

My problem is a bit more general, however (U(x,y) is defined numerically, it would be a Function on some space V). Does anyone have any experience solving a problem of this form. Specifically, I am unsure of how to write the weak variational form. I suspect it should look something like

http://imgur.com/phcBbKL

where v(x,y) is a trial function. In otherwords, I am wondering how to define the 2D general Schrodinger equation in FEniCs language, where the the Schrodinger Equation for a 1D harmonic oscillater has been previously defined as

#define problem
a = (inner(grad(u), grad(v)) \
     + x2*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)

Thanks

Comments

Hi, have you tried running the code? If you define x2 as Function(V) it should run okay. Could you provide some typical x2 and mesh?

---

Hi, I did not realise it would be that simple. I will replace x2 with my function and let you know if it works. Thanks

---

I believe I have it working, except for the boundary conditions, which don't seem to be obeyed. Thanks

---

Are you prescribing any boundary conditions? I assume you want u=0 everywhere on the boundary. You should then add something like this to your code

# V is the function space where you seek solution
bc = DirichletBC(V, Constant(0.), DomainBoundary())
bc.apply(A)
bc.apply(M)

Also, to preserve symmetry, you might want to create a fake right-hand side and then use
assemble_system

L = Constant(0.)*v*dx
A, _ = assemble_system(a, L, bc)  # bcs already applied
M, _ = assemble_system(m, L, bc) # bcs already applied

---

assemble_system is interesting. My calculation seems to run smoothly (with correct results) without it. Would there be a speed advantage in implementing it?

Currently it seems to be throwing an error at the SLEPcEigenSolver phase (I am probably not implementing it correctly)

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

#assemble stiffness matrix
L = Constant(0.)*v*dx
A, _ = assemble_system(a, L, bc)  # bcs already applied
M, _ = assemble_system(m, L, bc) # bcs already applied


#create eigensolver
eigensolver = SLEPcEigenSolver(A,M)

The last line gives an error:
TypeError: in method 'new_SLEPcEigenSolver', argument 1 of type 'dolfin::PETScMatrix const &'

Thanks again

---

Hi, as the error message says-the matrices must be PETScMatrix objects for the eigen-solver to work. With the above code assemble_system returns dolfin.cpp.la.Matrix. Sorry about that, this should work

A = PETScMatrix()
M = PETScMatrix()
_ = PETScVector()

assemble_system(a, L, bc, A_tensor=A, b_tensor=_)
assemble_system(m, L, bc, A_tensor=M, b_tensor=_)

The reason to use assemble_system is to preserve symmetry of A and M after the boundary conditions are applied.

---

Thanks, this worked with the caveat that I don't apply boundary condtions to M.

I.e. Here is a script I wrote to calculate and interactively plot the first five eigenfunctions for a quantum cylindrical infinite well. It works, but if I uncomment line 30 and comment line 31, the solver fails to correctly calculate the 4th and 5th eigenvalues.

from dolfin import *

#define mesh and function space
mesh = CircleMesh(Point(0,0),5,0.7)
V = FunctionSpace(mesh, 'Lagrange', 3)

#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)

Pot = Expression('0.0')

#define problem
a = (inner(grad(u), grad(v)) \
     + Pot*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,bc, A_tensor=M, 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
eigensolver.solve(5)

u = Function(V)
for i in range(0,5):
    #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)

---

Hi, I don't know quantum oscillators well enough to comment on the bcs. In the link you provided there were no bcs at all. Anyways, with your above code, if I set bcs on A and M, use linear instead of cubic polynomials and solve the whole spectrum then first three eigenvalues are same as with no bcs on M case, then there are few eigenvalues with r=1 and then I get correct 4th and 5th eigenvalues. In general, ignoring r=1, the spectrums are very similar. However, some eigenfunctions in the two cases have opposite signs.

---

Where are the eigenfunctions with eigenvalue=1.0 coming from exactly? If one applies boundary conditions to the solution of an eigenfunction, I would guess these need to be applied to both A and M—however only the case where M is left unchanged seems to work? Am I missing something here?