# Evolving nonlinear Schrodinger equation with higher-order algorithms?

First I will give the relevant information for my question, and then I'll ask the question.

$\large{\textrm{Background}}$

For evolving the nonlinear Schrodinger equation (NLS), one typically uses [a variant of] the Split-Operator, in which one split the nonlinear and linear components of the differential equation and treats them seperately, one after another. The version of the NLS I'm concerned with is

$$\frac{\partial \psi}{\partial t}=\frac{i}{2}\frac{\partial^2 \psi}{\partial x^2}+i|\psi|^2\psi=[\hat{D}+\hat{N}(x,t)]\psi$$

where we have defined

$$\psi=\psi(x,t),\,\,\hat{D}=\frac{i}{2}\frac{\partial^2}{\partial x^2},\,\,\hat{N}(x,t)=i|\psi|^2$$

As outlined quite succinctly in the previously cited wiki page, the general evolution $\psi(x,t)\to\psi(x,t+\Delta t)$ for sufficiently small step size $\Delta t$ is given by

$$\psi(x,t+\Delta t)=e^{\Delta t(\hat{D}+\hat{N})}\psi(x,t)$$

To actually evolve this with a computer, we can split/factor the exponential operator in the line above to a obtain an algorithm that is correct to order $(\Delta t)^2$.

$$\psi(x,t+\Delta t)\approx e^{\Delta t(\hat{D})}e^{\Delta t(\hat{N})}\psi(x,t)$$

Realizing that the second derivative in the $\hat{D}$ operator turns into a multiplication by $-k^2$ in Fourier-space, we can perform this evolution on a computer by calculating, in a discretized fashion of course,

$$\psi(x,t+\Delta t)\approx\mathcal{F}^{-1}\left(e^{-\frac{i}{2}\Delta t (k^2)}\mathcal{F}\left(e^{i\Delta t |\psi^2|}\psi(x,t)\right)[k,t]\right)[x,t]$$

This can be interpreted as follows: we first "evolve" the nonlinear part of the equation, and then we take a time step $\Delta t$ by the "evolution" of the linear part of the equation (in Fourier space). With this scheme, one can start with an initial wave-form/function $\psi(x,t_0)$ and evolve forward in time.

$\large{\textrm{Question}}$

How does one use higher order algorithms on the NLS? By using the Baker-Campell-Hausdorff formula for factorizing exponentials of operators, one can come up with a simple third order algorithm,

$$\psi(x,t+\Delta t)\approx e^{\frac{\Delta t}{2}\hat{D}}e^{\Delta t(\hat{N})}e^{\frac{\Delta t}{2}\hat{D}}\psi(x,t)$$

However, if we try to attack this in the same way we did before, we now have to evaluate $\hat{N}(x,t+\frac{\Delta t}{2})$, which involves evaluating $\psi(x,t+\frac{\Delta t}{2})$, which is the function we are trying to evolve! More precisely, if we only start off with an initial wave-form/function $\psi(x,t)$, then we don't have access to $\psi(x,t+\frac{\Delta t}{2})$. How can I resolve this problem.

Edit: After writing this I realized that by evaluating $e^{\frac{\Delta t}{2}\hat{D}}\psi(x,t)$, it's possible we are obtaining an approximation for $\psi(x,t+\frac{\Delta t}{2})$, which would answer my question affirmatively. However I'm going to keep this question up in case I'm incorrect.

Your added edit to the question points exactly the way this needs to be done: You first evolve your field using a half-step with the $\hat D$ operator, then a full step with the $\hat N$ operator, and then another half step with $\hat D$. These are done one after the other, i.e., the result of one operation is the input for the next.
• The evolution is then inevitably explicit (i.e. only function at time $t$ determines function at time $t+\Delta t$), which would seem to cause problems. Of course, this comes with the advantage of avoiding any difference-relations, making for a much "easier" and faster evolution. Do you know of any articles that discuss this? (or correct me if I'm under a misapprehension) May 2, 2016 at 22:39
• You see this wrong. In essence, you need to do one half-step with operator 1, then a full step with operator 2, then one half-step with operator 1. You can't think of this as going from $t$ to $t+\Delta t/2$, then $t+\Delta t/2$ to $t+3t+\Delta t/2$, and then $t+3t+\Delta t/2$ to $t+2\Delta t$. It's all part of just one time step of length $\Delta t$. May 3, 2016 at 1:56