I aim to solve the (non-linear) Schrodinger equation using the Crank-Nicolson method in Python. Here are my two functions.

def diags(N, dt, dx):
    Adiag = sp.empty(N)
    Asup = sp.empty(N)
    Asub = sp.empty(N)
    bdiag = sp.empty(N)
    bsup = sp.empty(N)
    bsub = sp.empty(N)

    Adiag.fill(1 + dt/dx**2)
    bdiag.fill(1 - dt/dx**2)

    A = sp.sparse.spdiags([Adiag, Asup, Asub], [0, 1, -1],N, N)
    b = sp.sparse.spdiags([bdiag, bsup, bsub], [0, 1, -1],N, N)
    return A, b

def crankNicolson(psi, dt, V, A, b, g=0):
    psi = np.exp(-dt*(V + g*np.abs(psi.reshape(-1))**2))*psi
    psi = sp.sparse.linalg.bicg(A, b*psi.reshape(-1).astype(np.complex128))[0]
    return psi

which has been taken from here. However, the given link aims to propagate in imaginary time, therefore I would need to have $ dt \rightarrow -i dt $ in order to convert the code to propagate in real time instead. This can simply be done by putting -1j*dt as an input into the functions instead of just dt.

To put simply, the function diags() is run once to obtain A and b, which are then inputs for the function crankNicolson(), which in turn is run within a loop to iterate the process and propagate it in time.

However when I do this and propagate, my situation does not evolve in time but remains stationary. Can anyone see what's up with it? It absolutely should evolve, as I compare it to an explicit-midpoint method of mine which evolves in time well using the exact same input parameters.

Some additional information:

The equation in dimensionless form is $$ i\frac{\partial \psi}{\partial t} = -\frac{1}{2}\frac{\partial^2}{\partial x^2}\psi + V \psi + g|\psi|^2\psi $$

where the $V$ term is a self-consistent potential which is calculated separately (it is not important for this question as it is definitely correct), and $g$ is the coupling constant. The boundary conditions are that $\psi(0) = 0$ and $\psi(L) = 0$, owing to the fact that this is a decomposition of a spherically symmetric system with wavefunction $\phi$ according to $\phi = \psi/x$ (though this detail is perhaps not important).

The initial condition is a half-Gaussian, given by $\psi = x \alpha \exp(-x^2/\beta)$, which is centred on $x=0$.

  • $\begingroup$ Crossposted from physics.stackexchange.com/q/781911/2451 $\endgroup$
    – Qmechanic
    Sep 26, 2023 at 6:24
  • 1
    $\begingroup$ Does your code work correctly with no nonlinearity or potential? If things aren't working, try to simplify until it does. $\endgroup$
    – whpowell96
    Sep 27, 2023 at 3:03

1 Answer 1


Here the link to a Colab notebook implementing Crank-Nicolson method in a finite difference method discretization of your problem: https://colab.research.google.com/drive/1u1G4Tk7ssk9BJYCksR5GrADkdCGWfC2p?usp=sharing


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.