Maintain unitary time evolution for a nonlinear ODE

I want to solve a nonlinear ODE of matrix $$A(t)$$ $$\mathrm{i}\dot A = A(t)M(t),\:\mathrm{with}\: M(t)=A^\dagger(t)H(t)A(t)$$ where $$H(t)$$ and hence $$M(t)$$ are Hermitian. Therefore, I presume the time evolution of $$A(t)$$ is unitary. Is there any algorithm that can maintain this unitarity? I've heard of Crank-Nicolson, but is it for linear ODE and $$M$$ is independent of $$t$$?

Crank-Nicholson is usually used for linear ODE but can also be extended to a nonlinear ODE

$$\dot z = f(z)$$

like so:

$$\frac{z_{n + 1} - z_n}{\delta t} = \frac{f(z_{n + 1}) + f(z_n)}{2}.$$

Another similar alternative is the implicit midpoint rule:

$$\frac{z_{n + 1} - z_n}{\delta t} = f\left(\frac{z_{n + 1} + z_n}{2}\right).$$

The Crank-Nicholson and implicit midpoint rules are the same for linear ODE, but for nonlinear ODE, they're different. The reason I mention the implicit midpoint rule is that it's a symplectic integration scheme, whereas the Crank-Nicholson method is not. If you haven't run across this concept before, I highly recommend reading Geometric Numerical Integration by Hairer, Lubich, and Wanner*.

One blunt approach you could take is to approximate $$M$$ as piecewise constant in each timestep and use the Pade approximant

$$e^{iz} \approx \frac{1 + \frac{1}{2}iz}{1 - \frac{1}{2}iz}$$

for $$e^{-i\delta tM}$$. This approximation is especially nice because it wraps the real line into the complex unit circle. When you then apply it to matrices, it takes things with real eigenvalues (like Hermitian matrices) and maps them into unitary matrices. You then get the following scheme:

\begin{align} A_{n + 1} & = A_n \cdot \left(I - \frac{1}{2}i\delta t M_n\right)\left(I + \frac{1}{2}i\delta tM_n\right)^{-1}\\ M_{n + 1} & = A_{n + 1}^\dagger H(t_{n + 1})A_{n + 1} \end{align}

If $$A_0$$ is a unitary matrix, then every successive iterate will be a unitary matrix (at least in real arithmetic). It's only a first-order accurate integration scheme, but the important part is that it preserves the structure of your system. Hopefully this is a good starting point from which you can derive more accurate schemes if need be.

*As Lutz Lehmann pointed out in the comments, your problem isn't a Hamiltonian system as such, so what I said about symplectic integrators is a bit of a non-sequitur. I mentioned this more just for general knowledge. Your particular problem has non-trivial geometric structure and so do Hamiltonian systems. If you often find yourself solving these kinds of problems then the knowledge you acquire about methods for Hamiltonian systems will transfer well to other types.

• Thanks for the nice and informative answer! I find it quite clarifying and the second part works for sure. But I have a question for the Crank-Nicholson and implicit midpoint rules you wrote: how to see that they preserve unitarity here? And the nonlinearity requires to know $z_{n+1}$ in order to evaluate $f(z_{n+1})$ and then update from $z_n$ from $z_{n+1}$, which looks weird. Do I misunderstand anything? Feb 14 '21 at 6:17
• The order of the factors in $A_{n+1}$ must be the other way around for the given equation $i\dot A=AM$. Also the $M$ in the quotient should have a factor $1/2$. Can you confirm that the scheme has still order 2? See my construction of a midpoint-like method in the cross-post math.stackexchange.com/questions/4019719. (Also, there is no Hamiltonian present in the meaning of classical mechanics, so a method being symplectic does not much to the situation.) Feb 14 '21 at 8:23
• @LutzLehmann Thanks for catching my silly mistake! I've updated the answer. Feb 15 '21 at 2:36
• @xiaohuamao You're correct, in both the Crank-Nicholson and implicit midpoint rules the next state $z_{n + 1}$ appears on both the left- and right-hand side of the discretized equation. You now have to solve a nonlinear equation at every time step. Implicit methods are much more expensive computationally than explicit methods, but in return you often get much better stability and structure-preserving qualities. Feb 15 '21 at 2:39
• Thanks. So your Crank-Nicholson rule doesn't work nor preserve unitarity here? (But yes, I like the general point you brought up.) Feb 15 '21 at 3:00