Asking advice for implementation of Conservative Finite Difference Scheme for numerically solving Gross-Pitaevskii equation

Question

I am trying to numerically solve the Gross-Pitaevskii equation for an impurity coupled with a one-dimensional weakly-interacting bosonic bath, given by (in dimensionless units):

\begin{align} i \frac{\partial \phi_0(x,t)}{\partial t} &= \left\{-\frac{1 + m}{2} \partial_{x}^2 + |\phi_0(x,t)|^2 + G\delta(x) + i V(t) \partial_{x} \right\}\phi_0(x,t)\nonumber\\ V(t) &= V_0 + i\frac{m}{\sqrt{\gamma}} \int dx \phi^*_0(x,t) \partial_{x} \phi_0(x,t) \end{align}

Where $m$,$\gamma$, $V_0$, and $G$ are external parameters which are real numbers. I would like to implement a conservative finite difference method, i.e., a finite difference method that discretizes the differential equation in both time $t$ and space $x$ and solve it while conserving the $L_2$ norm (number of particles), Energy and total momentum.

From studying different literatures (https://doi.org/10.3846/1392-6292.2009.14.109-126, https://www.mdpi.com/518624), I have found out that there is a finite discretization scheme (a ''modified version'' of the Crank-Nicholson scheme) that one can apply on a Gross-Pitaevskii equation which abides by the conservation rules. Without the drift term $i V(t)$, this scheme (derived using equation (3.1) of the first paper mentioned before) is given by:

\begin{align} i \frac{\phi_j^{n+1} - \phi_j^{n}}{\tau} = -\frac{1+m}{2} \left[ \frac{\phi_{j+1}^{n+1} - 2\phi_j^{n+1} + \phi_{j-1}^{n+1}}{2h^2} + \frac{\phi_{j+1}^n - 2\phi_j^n + \phi_{j-1}^n}{2h^2} \right] + \frac{1}{2}\left[\frac{G \delta(j)}{h} + \frac{1}{2} \left(|\phi_j^{n+1}|^2 + |\phi_j^n|^2 \right) \right] \left(\phi_j^{n+1} + \phi_j^n \right) \end{align}

Where $j$ and $n$ are discretized $x$ and $t$ respectively; and $h$ and $\tau$ denotes the values of space and time-steps, respectively. The initial conditions is given by $\phi^0_{j} = 1 \forall j \in [-L/2,L/2]$ and also periodic boundary conditions need to be implemented, given by $\phi_{L/2}^n = -\phi_{-L/2}^n \forall n \leq t_{\text{final}}$. I have two queries:

1. I have seen the Crank-Nicholson scheme to be implemented in the case of the heat-flow equation, where the linearity allows one to express the finite discretized equation to be written in the form of tridiagonal matrix and can be thus solved with Gaussian elimination processes (Thomas algorithm) (for a detailed example, check https://sepwww.stanford.edu/sep/prof/bei/fdm/paper_html/node15.html). But in the GPE case, the non-linearity prevents me from doing so. I would like to know how one can actually implement or solve this equation numerically.

2. Due to the physical nature of the problem, I would like the above-mentioned quantities to be also conserved in my case. With the addition of the drift term $iV(t)$, will the conservation laws still be respected with this conservative finite-difference scheme? I have checked that with the periodic boundary condition, the conservation laws are respected when there is no drift term.

Edit: I had a discussion with someone about these equations. Now I am curious if this set of non-linear coupled equations can be solved via Mathematica in some way?

Answer 1

1. Since the equation is nonlinear, you cannot reduce to a tridiagonal system of linear equations, but you can reduce to a tridiagonal system of nonlinear equations. If these are solved via Newton's method, it will reduce to a sequence of tridiagonal linear solves where the matrix at each step is the tridiagonal Jacobian of the function $F$ written in the form $F(\phi^{n+1})=0$.

2. This is a tricky subject. Generally, finite difference methods aren't guaranteed to conserve any quantities exactly. The only way to check is to write down a discretized version of the conserved quantity then check to see if applying one step of the method preserves its value or not. This becomes more problematic when you introduce drift terms, as the stability of the drift discretization can change depending on the sign of $V$ and more advanced techniques are needed. Finite volume methods, which are designed for conservation laws, are one option, and some other options for getting better conservation in your simulation are discussed in Chris Rackaukas's great answer here.