0
$\begingroup$

I am trying to solve matrix ODE:

(*change global precision*)
$PreRead = (# /.s_String /; StringMatchQ[s, NumberString] && Precision@ToExpression@s == MachinePrecision :> s <> "`44." &);
(*define matrix components*)
ρ[t_] := Table[a[i, j][t], {i, 4}, {j, 4}] 
slnn = Flatten[Table[a[i, j][t], {i, 4}, {j, 4}]];

(*define initial condition for ODE*)
init = ( {{0, 0, 0, 0},{0, 1, 0, 0},{0, 0, 0, 0},{0, 0, 0, 0}} ); γ = 1; θ = Pi/2;
(*now define constants*)
r[i_, j_] := (i - j) di; di = 1/2*10^-7; μ = 10^-10; ω0 =2*10^15; c = 3*10^8; N[2 Pi/(ω0/c)]
σ1p = ( {{0, 0, 1, 0},{0, 0, 0, 1},{0, 0, 0, 0},{0, 0, 0, 0}} ); σ1m = ( {{0, 0, 0, 0},{0, 0, 0, 0},{1, 0, 0, 0},{0, 1, 0, 0}} ); σ2p = ( {{0, 1, 0, 0},{0, 0, 0, 0},{0, 0, 0, 1},{0, 0, 0, 0}} ); σ2m = ( {{0, 0, 0, 0},{1, 0, 0, 0},{0, 0, 0, 0},{0, 0, 1, 0}} );
Λ12 = 3/4 (24 Cos[1/3] + 9 Sin[1/3]);
Λ21 = 3/4 (24 Cos[1/3] + 9 Sin[1/3]);
sh = Sinh[2]; ch = Cosh[2];
Λp21 = -0.47260812;
Λp11 = -0.47260812;
(*now we solve ODE*)
solsq = NDSolve[{ρ'[
      t] == -sh ch (((1/2 + 
              I Λp21) (σ2p.σ1p.ρ[
                t] + ρ[t].σ2p.σ1p) + (1/2 + 
              I Λp21) (σ1p.σ2p.ρ[
                t] + ρ[t].σ1p.σ2p)) + ((1/2 - 
              I Λp21) (σ2m.σ1m.ρ[
                t] + ρ[t].σ2m.σ1m) + (1/2 - 
              I Λp21) (σ1m.σ2m.ρ[
                t] + ρ[t].σ1m.σ2m)) - ((1 + 
              2 I Λp21) σ2p.ρ[
              t].σ1p + (1 + 
              2 I Λp21) σ1p.ρ[
              t].σ2p + (1 + 
              2 I Λp11) σ1p.ρ[
              t].σ1p + (1 + 
              2 I Λp11) σ2p.ρ[
              t].σ2p) - ((1 - 
              2 I Λp21) σ1m.ρ[
              t].σ2m + (1 - 
              2 I Λp21) σ2m.ρ[
              t].σ1m + (1 - 
              2 I Λp11) σ1m.ρ[
              t].σ1m + (1 - 
              2 I Λp11) σ2m.ρ[
              t].σ2m)) - 
      1/2 (sh + 
         1) ((ρ[t].σ1p.σ1m + ρ[
             t].σ1p.σ2m + ρ[
             t].σ2p.σ1m + ρ[
             t].σ2p.σ2m) + \
(σ1p.σ1m.ρ[t] + σ1p.σ2m.ρ[
             t] + σ2p.σ1m.ρ[
             t] + σ2p.σ2m.ρ[t]) - 
     2 (σ1m.ρ[t].σ1p + σ1m.ρ[
          t].σ2p + σ2m.ρ[
          t].σ1p + σ2m.ρ[t].σ2p)) - 
  1/2 sh ((ρ[t].σ1m.σ1p + ρ[
         t].σ1m.σ2p + ρ[
         t].σ2m.σ1p + ρ[
         t].σ2m.σ2p) + \
(σ1m.σ1p.ρ[t] + σ1m.σ2p.ρ[
             t] + σ2m.σ1p.ρ[
             t] + σ2m.σ2p.ρ[t]) - 
         2 (σ1p.ρ[t].σ1m + σ1p.ρ[
              t].σ2m + σ2p.ρ[
              t].σ1m + σ2p.ρ[t].σ2m)) - 
      I Λ21 (σ1p.σ2m.ρ[t] - ρ[
           t].σ1p.σ2m) - 
      I Λ21 (σ2p.σ1m.ρ[t] - ρ[
           t].σ2p.σ1m) , ρ[0] == init}, 
   slnn, {t, 0, 1.6}, WorkingPrecision -> 44];
(*plot*)
Plot[{Evaluate[a[2, 2][t] /. solsq], Evaluate[a[4, 4][t] /. solsq], 
  Evaluate[a[1, 1][t] /. solsq], Evaluate[a[3, 3][t] /. solsq], 
  Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] + 
   Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq]}, {t,
   0, 1.1}, PlotRange -> {-2, 6}, 
 PlotLegends -> {"ud", "dd", "uu", "du", "sum"}]

Now I get the result: enter image description here

which cannot be correct, because we can prove that:

1,"uu","ud","du","dd" must be between 0 and 1

2, "sum" should be a constant, it is the trace of matrix ρ, and it can be easily proved by taking the trace of both sides of our ODE. when we take the trace, we have

0.*10^-43 I a[1, 4][t] + 0.*10^-43 I a[4, 1][t]

If we increase the WorkingPrecision in the first line of the code, we can make it closer to 0. However, the plot still shows a very bad behavior no matter how I change the WorkingPrecision.

I try to pick the value of solution:

Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] +Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq] /. t -> 1
Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] +Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq] /. t -> 1.5
Evaluate[a[2, 2][t] /. solsq] + Evaluate[a[4, 4][t] /. solsq] +Evaluate[a[1, 1][t] /. solsq] + Evaluate[a[3, 3][t] /. solsq] /. t -> 2

They give:

{1.00000000000000000000}
{1.000000000}
no significant digits available to display

That's why the solution seems so badly behaved, but I can't fix it. Increasing the working precision seems not to help. Can anyone help me solve this numerical issue? Thanks!

$\endgroup$
  • $\begingroup$ I noticed that if I change parameters from sh=Sinh[2] to sh=Sinh[1], the solution will be better. However, I literally need the solution for sh=Sinh[2]. $\endgroup$ – Jieyu You Sep 21 '16 at 18:19
  • $\begingroup$ Early in the code is the statement N[2 Pi/(ω0/c)], which appears to do nothing. Why is it there? $\endgroup$ – bbgodfrey Sep 22 '16 at 4:24
  • $\begingroup$ If the components of a are supposed to be less than one in absolute value, then you have a bug in your ODE, because this condition already is violated by t == 0.1, and WorkingPrecision clearly is not at issue. By the way, I suggest you use Λp21 = Rationalize[-0.47260812, 10^-15], and similarly for Λp11. $\endgroup$ – bbgodfrey Sep 22 '16 at 4:57
  • $\begingroup$ @bbgodfrey, Oh, I just wanted to test if the precision was indeed increased, and I forgot to delete it. $\endgroup$ – Jieyu You Sep 22 '16 at 15:36
  • $\begingroup$ I see that my last comment may be unclear. I suggested that you use Λp21 = Rationalize[-0.47260812, 10^-15], because not doing so can interfere with high WorkingPrecision. However, this suggestion will not resolve your problem, which I believe is not associated with precision issues. Instead, I believe that there is a typo in your ODE. $\endgroup$ – bbgodfrey Sep 22 '16 at 15:42