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 following image given by our supervisor.

enter image description here

This is however 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'), 
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)


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

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

We've also tried our luck doing as the example found in this link without success: Gross-Pitaevskii in FEniCS

What are we doing wrong here?

  • 2
    $\begingroup$ This seems like a question on how to use fenics, rather than a general scicomp question. You will get better answers if you ask on the fenics forums fenicsproject.discourse.group $\endgroup$ – cfdlab Feb 15 at 9:39
  • $\begingroup$ Discretizing the real and imaginary parts of a Schrodinger equation separately is something new to me... This is not how people solve time-dependent Schrodinger equation. Have you heard about time-domain finite difference? Look at here: iopscience.iop.org/article/10.1088/1751-8113/40/8/013 $\endgroup$ – Alone Programmer Feb 15 at 18:17
  • $\begingroup$ Unfortunately FEniCS does not support complex numbers, this is simply a workaround. $\endgroup$ – Anton Scotte Feb 15 at 22:15
  • $\begingroup$ @AntonScotte you don’t need any complex number if you use time-domain finite difference. $\endgroup$ – Alone Programmer Feb 15 at 23:24
  • $\begingroup$ Also, it seems DOLFINX supports complex numbers: igorbaratta.github.io/gsoc_blog/jekyll/update/2018/07/03/… $\endgroup$ – Alone Programmer Feb 16 at 0:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.