I've devising a Crank-Nicolson scheme for the Schrodinger equation using a time-independent potential, $V(x)$. In devising this scheme the following equation is approximated and THEN discretized, where the following operator $U(t) = \exp\left( -i H t)\right)$ is the unitary evolution operator. \begin{equation} \psi(t+\Delta t,x) = U(\Delta t) \psi(t,x) = \textrm{e}^{-iH\Delta t}\psi(t,x) \approx \frac{1 - i\frac{1}{2}H\Delta t}{1 + i\frac{1}{2}H\Delta t} \psi(t,x) \:\:\:\:\:\:\:\: (\textrm{eq. 1}) \end{equation} The approximation of the operator is called Cayley's approximation, and it designed to maintain the unitary aspect of the evolution. My question pertains now to how to appropriately compute the error of the scheme. If I evolve the system in time by some amount $T$, then I would normally assume that the error of the scheme would be computed as $$ \varepsilon = \sqrt{\Delta x}\left\lVert U(T) \psi_0(x_j) - \phi(T,x_j) \right\rVert_2 $$ Where, $\psi_0(x)$ is the initial condition of the problem, $x_j$ are the grid points in the $x$-domain, and $\phi(t,x)$ would be the solution to the scheme derived from the rightmost expression in the eq. 1. I've been calculating the errors for the scheme, from these errors, I've calculated the order of accuracy (for Crank-Nicolson is said to be 2 in time and 2 in space). However, when I define my computational gird as $\Delta x = \Delta t^2$ (meaning the dominating errors will be due to the error in time), my order of accuracy (in time) roughly oscillates around an order of accuracy of about 0.25 (far from the order predicted). I've hypothesized that is might because, I'm comparing the evolution of the scheme against the non-approximated evolution operator, $U(\Delta t)$. However, I would conjecture, that to find the true order of accuracy, I would need to apply the approximation of the evolution operator to the initial state, $\psi_0(x)$. Meaning the error would be instead $$ \varepsilon = \sqrt{\Delta x}\left\lVert \frac{1 - i\frac{1}{2}HT}{1 + i\frac{1}{2}HT} \psi_0(x_j) - \phi(T,x_j) \right\rVert_2 $$ I'm not sure, however, this is fully justified. Supposing that it is though, then I lose a convenient computational tool, in that the eigenvalue of the evolution operator is $\lambda = \textrm{e}^{-iE_nt}$, where $E_n$ is the energy of the particular state being evolved. Therefore, computation of the error only involves multiplying $\psi_0(x)$ by $\textrm{e}^{-iE_nT}$. I imagine that I'd need to find the eigenvalue for the the approximation to the evolution, $\lambda_a(\Delta t)$ to that would allow me to calculate $$ \varepsilon = \sqrt{\Delta x}\left\lVert \lambda_a(T) \psi_0(x_j) - \phi(T,x_j)\right\rVert_2 $$ Effectively, this all comes down to answering a few questions
Since I'm approximating the evolution operator, does my error need to calculated with that operator applied to the initial sate?
Second, assuming that the first point is true, how would I go about solving for the eigen value?
There's a lot to unpack in this problem, so it's reasonable if you ask about whether or not my scheme is actually implemented correctly. However, I think the legitimacy of the question will remain valid regardless of whether my code is correct (I've been debugging for quite some time which is why I'm pursuing this avenue of inquiry). This is why I'm not posting the code as this has more to do about calculating errors for the discretization of an approximation for an operator.