Conservation of probability issue when solving ODE

Question

I am trying to solve this two level (Schrodinger) equation as a function of time:

$i\begin{pmatrix} \dot{x}\\ \dot{y} \end{pmatrix} = \begin{pmatrix} 0 & iW+dE_0\sin(\omega t)\\ -iW+dE_0\sin(\omega t) & \Delta \end{pmatrix}\begin{pmatrix} x\\ y \end{pmatrix}$

(I can go into more details about the Hamiltonian if needed). The initial conditions are $x(0)=1$, $y(0)=0$. This is my Mathematica code for it:

W = 10;
OmegaRabi = 1000;
omega = 700000;
delta = 10000;
eqns = {I*x'[t] == (I*W + OmegaRabi*Sin[omega*t])*y[t],
   I*y'[t] == (-I*W + OmegaRabi*Sin[omega*t])*x[t] + y[t]*delta, 
   x[0] == 1, y[0] == 0};
sol = NDSolve[eqns, {x, y}, {t, 0, 4}][[1]]
Plot[{Abs[x[t]]^2 + Abs[y[t]]^2 /. sol}, {t, 0, 1}]

All the variables are as in the original equation except for OmegaRabi which is equal to $dE_0$. The output of this code (which would be the total probability i.e. $|x(t)|^2+|y(t)^2|$) should be constant 1. However I get what is seen in the plot below. I assume this has to do with some rounding numerical errors (?). Is there a way to fix it and ideally have it deviate from 1 much slower as a function of time? I am interested in the $|y(t)|^2$ as a function of time, and if one looks at that, the values for it are around 0.00002, so currently the error after one second seems to be 5 times bigger than the value I am interested in.

Answer 1

omega = 700000; is a rather high frequency, so let's manually limit the step size a bit:

soldense = 
   NDSolve[eqns, {x, y}, {t, 0, 4}, MaxSteps -> Infinity, 
     MaxStepSize -> 1/omega/10][[1]]; // AbsoluteTiming
(* {165.859, Null} *)

Plot[Abs[x[t]]^2 + Abs[y[t]]^2 - 1 /. soldense // Evaluate, {t, 0, 1}]

I doubt if this makes much sense though. The following is a comparison for the default y[t] (note: one probably needs to add MaxSteps -> Infinity to OP's code) and MaxStepSize -> 1/omega/10 y[t]. As we can see, the difference is negligible:

ReImPlot[y[t] /. {sol, soldense} // Evaluate, {t, 4 - 10^-4, 4}]

Answer 2

@xzczd has already pointed out how to deal with the major issue, the high value of omega, by limiting the step size. The main interest in this post is the use of the "Splitting" method. There's rather little explanation of how to construct an effective splitting method. While the basic concept of splitting a vector field into a sum of subfields seems easy and appropriate for this problem, there is some but not much information about the requirements on the subfields and methods for each subfield. Since the attempt below seems to succeed, I thought it worth adding as an example. It's also a bit faster and more accurate than the default method. I don't know if the overall "DifferenceOrder" matters or what effect it has. Changing it to extreme values made no difference. The computed solution has a fixed step size. I don't know if that's a feature of the splitting method or a feature of this problem. The consistent oscillations of the solution makes a fixed-step method seem a reasonable approach.

W = 10;
OmegaRabi = 1000;
omega = 700000;
delta = 10000;
eqns = { (* t <- time[t] to make autonomous *)
   I*x'[t] == (I*W + OmegaRabi*Sin[omega*time[t]])*y[t], 
   I*y'[t] == (-I*W + OmegaRabi*Sin[omega*time[t]])*x[t] + 
     y[t]*delta,
   time'[t] == 1};
ics = {x[0] == 1, y[0] == 0, time[0] == 0};
F1 = eqns;
F1[[2]] = F1[[2]] /. y[t] -> 0;
F1[[3, 2]] = 0;
F2 = eqns;
F2[[1]] = F2[[1]] /. y[t] -> 0;
F2[[2]] = F2[[2]] /. x[t] -> 0;
F1  (* splitting of the vector field *)
F2
PrintTemporary@Dynamic@{Clock@Infinity};
sol = NDSolve[{eqns, ics}, {x, y}, {t, 0, 4}
     , Method -> {"Splitting", "DifferenceOrder" -> 8, 
       "Equations" -> {F1, F2},
       "Method" -> {
         {"ExplicitRungeKutta", "DifferenceOrder" -> 4}, 
         {"ExplicitRungeKutta", "DifferenceOrder" -> 4}}}
     , StartingStepSize -> 1/omega/2, MaxSteps -> 10^8
     ][[1]]; // AbsoluteTiming
Plot @@ {Abs[x[t]]^2 + Abs[y[t]]^2 - 1 /. sol, 
  Flatten@{t, x@"Domain" /. First@sol}}

(*
{I x'[t] == (10 I + 1000 Sin[700000 time[t]]) y[t], 
 I y'[t] == (-10 I + 1000 Sin[700000 time[t]]) x[t], 
 time'[t] == 0}

{I x'[t] == 0,
 I y'[t] == 10000 y[t], 
 time'[t] == 1}

{18.4591, Null}
*)

Update:

Using the submethods "Method" -> {"LocallyExact", "LocallyExact"} in the "Splitting" method above lowers the time to 8 sec. and produces an amazingly identical error plot of Abs[x[t]]^2 + Abs[y[t]]^2 - 1.