# Matrix ODE gives badly-behaved numerical solution [duplicate]

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; ch = Cosh;
Λ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) , ρ == 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: 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!

• I noticed that if I change parameters from sh=Sinh to sh=Sinh, the solution will be better. However, I literally need the solution for sh=Sinh. Sep 21, 2016 at 18:19
• Early in the code is the statement N[2 Pi/(ω0/c)], which appears to do nothing. Why is it there? Sep 22, 2016 at 4:24
• 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. Sep 22, 2016 at 4:57
• @bbgodfrey, Oh, I just wanted to test if the precision was indeed increased, and I forgot to delete it. Sep 22, 2016 at 15:36
• 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. Sep 22, 2016 at 15:42