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)
Asup.fill(-dt/(2*dx**2))
Asub.fill(-dt/(2*dx**2))
bdiag.fill(1 - dt/dx**2)
bsup.fill(dt/(2*dx**2))
bsub.fill(dt/(2*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$.
