Time-dependent Schrodinger equation implementation in FEniCS

mathematics
1

For our Bachelors thesis we’re trying to solve the Schrodinger equation i\partial_tu = -\nabla^2u+Vu in FEniCS. Given the domain [-5, 5]^2 with an initial value of u_0(x, y)=e^{(-2(x^2+y^2))} and with V(x,y)=\frac{x^2}{4}+9y^2 at time T=2 the solution |u|^2 should look like the image in this post, given by our supervisor https://scicomp.stackexchange.com/questions/34443/time-dependent-schrodinger-equation-implementation-in-fenics

This is the result we get:

enter image description here

The equation is dicretized with Crank-Nicolson as follows:
(where u=p+qi)

For the real part:

\frac{\Delta t}{2}(\nabla^2p^{n+1}-Vp^{n+1})-q^{n+1}=-q^n + \frac{\Delta t}{2}(Vp^n-\nabla^2p^n)

for the imaginary part:

\frac{\Delta t}{2}(\nabla^2q^{n+1}-Vq^{n+1})+p^{n+1}=-p^n + \frac{\Delta t}{2}(\nabla^2q^n-Vq^n)

In the code below, u[0] is p and u[1] is q

We’re having trouble creating the correct residual when solving for u.

In FEniCS the equation is written in weak form and when multiplying with the test function we change sign of dot(grad(), grad()) to accommodate for multiplication with i. This is were we probably are at fault because we do not know how to construct the residual. When changing the sign of the gradient for the imaginary case the solution does not converge.

T = 2.0            # final time
num_steps = 200  # number of time steps
dt = T / num_steps # time step size

mesh = RectangleMesh(Point(-5, -5), Point(5, 5), 55, 55)
V = VectorFunctionSpace(mesh, 'Lagrange', 1, dim=2)
# Define initial value
u_0 = Expression(  ( 'exp(-2*(pow(x[0], 2)+pow(x[1], 2)))', '0'), 
degree=1)
u_n = interpolate(u_0, V)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V,(0, 0), boundary)

u = Function(V)
v = TestFunction(V)
f = Expression(('(1/4)*pow(x[0], 2) + 9*pow(x[1], 2)', '0'), degree=2)

t=0

#Residual
ReVL = -u[1]*v[0]*dx + (dt/2)*(-dot(grad(u[0]), grad(v[0])) - 
f[0]*u[0]*v[0])*dx
ReHL = -u_n[1]*v[0]*dx + (dt/2)*(f[0]*u_n[0]*v[0] + dot(grad(u_n[0]), 
grad(v[0])))*dx
ImVL = u[0]*v[1]*dx + (dt/2)*(-dot(grad(u[1]), grad(v[1])) - 
f[0]*u[1]*v[1])*dx
ImHL = u_n[0]*v[1]*dx + dt/2*(f[0]*u_n[1]*v[1] + dot(grad(u_n[1]), 
grad(v[1])))*dx

FReal = ReVL - ReHL
FIm = ImVL - ImHL

Fny = FReal + FIm

J = derivative(Fny, u)
for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(Fny == 0, u, bc, J=J)
    # Update previous solution
    u_n.assign(u)

What are we doing wrong here?

2

Hi, there are a few things I’ve noticed:

  • there is a sign error on the right-hand side of your equation for the imaginary part (I believe the whole RHS should be -1x what it currently is);

  • if you want your solution to match your supervisor’s, you should begin by swapping the x and y components in your V function to get the vertical, rather than horizontal, alignment. Also, experiment with running for different times T to get agreement.

  • the equation you are solving is linear so you can use the more standard finite element setup with the solve command, rather than obtaining the Jacobian matrix. I think this will be much simpler.

3

Can you please show how this is plotted?