Why does NDSolveValues give different results?

Question

I'm trying to solve this differential equation to obtain the density matrix elements $\rho_{22}$.

$$\frac{d\rho}{dt}=-i\Big(H_{0}(t)\rho(t)-\rho(t)H^{\dagger}_{0}(t)\Big)$$,

for non-Hermitian Hamiltonian $$H_{0}(t)=\begin{pmatrix}-v t & -{\it i}A \\{\it i}B & v t \\\end{pmatrix},$$

with $$H^{\dagger}_{0}(t)=\begin{pmatrix}-v t & -{\it i}B \\{\it i}A & v t \\\end{pmatrix},$$ my code works properly

        A=0.1;
        B=1;
        v=1;
        t1=-100;
        t2=0.5;    

    H0[t_] = {{-v*t , -I*A}, {I*B, v*t}};

    H0dag[t_] = {{-v*t , -I*B}, {I*A, v*t}};

Solnh = NDSolveValue[{D[rho[t], t] == -I*(H0[t].rho[t] - rho[t].H0dag[t]),rho[t1] == {{1, 0}, {0, 0}}}, rho, {t, t1, t2}];

    Sol1=Solnh[t2];    

    Print[Sol1[[2,2]]]

the results is $\rho_{22}=1.53105 $

I know that it is possible to transform the non-Hermitian Hamiltonain $H_0(t)$ to the Hermitian Hamiltonian by using a non unitary transformation

$$H_h=\begin{pmatrix}\sqrt{\frac{B}{A}} & 0 \\0 & 1 \\\end{pmatrix}.\begin{pmatrix}-v t & -{\it i}A \\{\it i}B & v t \\\end{pmatrix}.\begin{pmatrix}\sqrt{\frac{A}{B}} & 0 \\0 & 1 \\\end{pmatrix}=\begin{pmatrix}-v t & -{\it i}\sqrt{AB} \\{\it i}\sqrt{AB} & v t \\\end{pmatrix}$$.

By applying the same non-unitary transformation on the density matrix we have

$$\rho_h=\begin{pmatrix}\sqrt{\frac{b}{a}} & 0 \\0 & 1 \\\end{pmatrix}.\begin{pmatrix}\rho_{11} & \rho_{12} \\\rho_{21} & \rho_{22} \\\end{pmatrix}.\begin{pmatrix}\sqrt{\frac{a}{b}} & 0 \\0 & 1 \\\end{pmatrix}=\begin{pmatrix}\rho_{11} & \sqrt{\frac{B}{A}}\rho_{12} \\\sqrt{\frac{A}{B}}\rho_{21} & \rho_{22} \\\end{pmatrix}$$.

So $\rho_{22}$ does not change under transformation. Therefore if we solve the master equation for Hermitian Hamiltonian $$H_h=\begin{pmatrix}-v t & -{\it i}\sqrt{AB} \\{\it i}\sqrt{AB} & v t \\\end{pmatrix}=H_h^{\dagger}$$, i.e.,

$$\frac{d\rho_h}{dt}=-i\Big(H_{h}(t)\rho_h(t)-\rho_h(t)H_{h}(t)\Big)$$,

$\rho_{22}$ should be the same as the previous result.

My code is

            A=0.1;
            B=1;
            v=1;
            t1=-100;
            t2=0.5;
    
        Hh[t_] = {{-v*t , -I*Sqrt[A*B]}, {I*Sqrt[A*B], v*t}};
    
    Solh = NDSolveValue[{D[rhoh[t], t] == -I*(Hh[t].rhoh[t] - rhoh[t].Hh[t]),
rhoh[t1] == {{1, 0}, {0, 0}}}, rhoh, {t, t1, t2}];

        Sol2=Solh[t2];  
  
        Print[Sol2[[2,2]]]

which gives $\rho_{h,22}=0.153116 $. Why the results are different in both methods? If we check $\rho_{11}$ which does not changes under transformation, it remains intact and in both approach $\rho_{11}=0.8468$.

I was wondering if someone could help me. Thank you for your kind attention.

Answer 1

Too long for comment.

Well, your new ode is different after transformation. Your claim is that the solution should remain the same due to how the math is supposed to work. This is really a math question. This is your ode in first case

ClearAll["Global`*"]
A = 0.1;
B = 1;
v = 1;
t1 = -100;
t2 = 0.5;

H0[t_] = {{-v*t, -I*A}, {I*B, v*t}};

H0dag[t_] = {{-v*t, -I*B}, {I*A, v*t}};
ode = D[rho[t], t] == -I*(H0[t] . rho[t] - rho[t] . H0dag[t]);
TraditionalForm[ode]
ic = rho[t1] == {{1, 0}, {0, 0}};

Solnh = NDSolveValue[{ode, ic}, rho, {t, t1, t2}];
Sol1 = Solnh[t2];
Print[Sol1[[2, 2]]]

And this is your second case

ClearAll["Global`*"]
A = 0.1;
B = 1;
v = 1;
t1 = -100;
t2 = 0.5;

Hh[t_] = {{-v*t, -I*Sqrt[A*B]}, {I*Sqrt[A*B], v*t}};
ode = D[rhoh[t], t] == -I*(Hh[t] . rhoh[t] - rhoh[t] . Hh[t]);
TraditionalForm[ode]
ic = rhoh[t1] == {{1, 0}, {0, 0}};

Solh = NDSolveValue[{ode, ic}, rhoh, {t, t1, t2}];
Sol2 = Solh[t2];
Print[Sol2[[2, 2]]]

You see the ode system is different in both cases.

Your claim is that the output should be the same, because the "math" should not have caused this to change.

So either NDSolve is wrong, or your math has some faulty assumptions in it.

Are your sure your math is correct and the result should be the same? May be you could check on this again since the difference is not negligible, I do not think this is due to round off errors, or may be it is? it is of order 0.0001 difference.

Did you try DSolve ?

Answer 2

If you are sure your math is correct, you may try to increase accuracy and precision:

AccuracyGoal -> 24
PrecisionGoal -> 10

See more details here.

Having said that, really make sure that your analytical solution (algebraic transformation) is identical.