I am new to the field of computational physics and have a couple of questions regarding solving the non-linear Schrdinger equation using Operator splitting.1) If the hamiltonian is of the form $H=\frac{\partial^{2}}{\partial x^{2}}+\gamma|\psi|^{2}$ then the standard procedure I understand is to exponentiate $-i(\gamma |\psi|^{2})\Delta t/\hbar$ and operate it on the initial value of $\psi$, then take a fourier transform to convert it to momentum space and operate it with exponential of $-ip^{2}\Delta t/\hbar$ and convert the resultant back to position space. We repeat this for each time interval $\Delta t$.Instead, why can't we do everything in the momentum space to begin with? Why this back and forth shifting from position to momentum space?2) Suppose now I have an additional term of $\frac{\partial^{2}}{\partial x^{2}}|\psi|^{2}$ in the Hamiltonian, then how do I accomodate this term in the scheme of split operator method?
Split operator method
nonlinear equations;fourier analysis;operator splitting
The Strang splitting method goes like this. You start with the PDE$$ u_t = (L+B)u, \qquad L = \partial_x^2, $$and you notice that when $L$ and $B$ are independent of $x$, the exact solution after time $\delta t$ to this is$$ u = e^{(L+B)\delta t}u_0 \approx e^{\frac12 L\delta t}e^{B\delta t}e^{\frac12 L\delta t}u_0. $$Because $B$ is just a function of $x$, the operator $e^{B\delta t}$ just multiplies by $e^{B(x)\delta t}$.Because $L = F \hat L F^{-1}$ is diagonalizable by the Fourier transform (assuming the right boundary conditions),$$ e^{L\delta t} = F e^{\hat L \delta t} F^{-1},$$where $e^{\hat L\delta t}$ is the operator that multiplies each Fourier mode by $e^{-k^2\delta t}$ (this depends on choice of normalization).This is why Fourier transforms are done at each step: in the Fourier basis, and only in that basis, is $L$ diagonal, which makes it trivial to compute its exponential.For your equation, to get what $e^{B\delta t}$ would look like, you write out the relevant portion of the equation, with only $B$ present:$$ i\hbar u_t = (\gamma |u|^2 + \alpha (|u|^2)_{xx})u. $$One thing you could do is to approximate $g(x) \approx (|u_0|^2)_{xx}$, so that$$ u(t,x) \approx u(0,x)\exp\left(\frac{\gamma|u(0,x)|^2 + \alpha g(x)}{i\hbar}\,\delta t\right). $$Because of the nonlinearity, it might work, but I think there isn't a guarantee that it willI haven't tried it. But the idea is still the same: split the r.h.s. into two operators, and for each operator solve the corresponding PDE, choosing the operators in a way that makes this step easy.