Conservation of probability issue when solving ODE in Mathematica

Click For Summary
The discussion focuses on solving a two-level Schrödinger equation using Mathematica, with specific initial conditions and parameters. The equation involves time-dependent terms and is set up to calculate the probabilities of two states, represented by functions x(t) and y(t). The user notes that the expected total probability, given by the sum of the squares of the magnitudes of x(t) and y(t), should remain constant at 1. However, the output from the Mathematica code shows deviations from this expected value, likely due to numerical errors. The user seeks advice on how to improve the accuracy of the solution, particularly in relation to the probability of state y(t), which is currently yielding small values and significant error. Suggestions include adjusting the AccuracyGoal and PrecisionGoal parameters in the NDSolve function to enhance the numerical stability and accuracy of the results over time.
kelly0303
Messages
573
Reaction score
33
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_0sin(\omega t)\\
-iW+dE_0sin(\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.
 

Attachments

  • schrodinger.png
    schrodinger.png
    3.6 KB · Views: 136
Technology news on Phys.org
In NDSolve, you should modify AccuracyGoal and/or PrecisionGoal to try and get better results.