# Homogeneous linear ODEs with time-dependent coefficient matrix, optimalisation of NDSolve algorithm

I would like to ask for advice on the following (probably not too complicated) problem:

I have a linear system of ODEs of the following form:

$$\dot{X}(t) = M(t)X(t)$$

where $$X(t)$$ is a vector, $$M(t)$$ is a matrix of $$n \times n$$ size. The algorithm for the equations itself is:

n = 2 (1 + 2*size)^2;
X[t_] = Table[xi[i, t], {i, 1, n}];
X0 = Table[0, {n}];
X0[[Round[n/2]]] = 1;
diffEq = Table[
X'[t][[i]] == -0.5 I*(Evomatrix[t].X[t])[[i]], {i, 1, n}];
initialCond = Table[X[0][[i]] == X0[[i]], {i, 1, n}];
sol = NDSolve[{Join[diffEq, initialCond]}, X[t], {t, 0, 1},
Method -> "BDF",
PrecisionGoal -> 0, AccuracyGoal -> 0, MaxStepSize -> 0.00001]

Plot[Abs[Evaluate[Table[xi[i, t] /. sol, {i, 1, 2}]]], {t, 0, 1},
PlotRange -> All]
Plot[Evaluate[Abs[xi[1, t] /. sol]^2 + Abs[xi[2, t] /. sol]^2 ], {t,
0, 1}, PlotRange -> All]


which works perfectly with only two variable, but does not finish running even under an hour with 18 variables.

The time-dependent matrix itself is defined in a rather complicated way, also involving matrix-inversion (I tried using LinearSolve but I did not notice any difference). Even still, the matrix is evaluated for concrete values in about a second, so I imagine that with proper stepsize it should work out. Anyway here are the codes:

gamma = -0.005;
omega = 1;
size = 0;
eltolas = 200;
Szorz[elso_, masod_] :=
Exp[-0.5 (Abs[elso])^2 - 0.5 (Abs[masod])^2 + Conjugate[elso]*masod ]


Here "size" is the important (integer) parameter. With 0, the algorithm runs quickly, but if it is 1, I already have a problem.

l := Sqrt[Pi]*(Mod[y - 1, 1 + 2*size] - size + Re[eltolas]);
m := Sqrt[Pi]*Floor[(y - 1)/(1 + 2*size) - size + Im[eltolas]];
s := Sqrt[Pi]*(Mod[x - 1, 1 + 2*size] - size + Re[eltolas]);
p := Sqrt[Pi]*Floor[(x - 1)/(1 + 2*size) - size + Im[eltolas]];
BRAalfaP[t_] := gamma + (l + I*m - gamma) Exp[-I*omega*t];
KETalfaP[t_] := gamma + (s + I*p - gamma) Exp[-I*omega*t];
BRAalfaN[t_] := -gamma + (l + I*m + gamma) Exp[-I*omega*t];
KETalfaN[t_] := -gamma + (s + I*p + gamma) Exp[-I*omega*t];
BRAGammaP[t_] := gamma*Im[l + I*m - (l + I*m - gamma) Exp[-I*omega*t]];
KETGammaP[t_] := gamma*Im[s + I*p - (s + I*p - gamma) Exp[-I*omega*t]];
BRAGammaN[t_] := -gamma*
Im[l + I*m - (l + I*m + gamma) Exp[-I*omega*t]];

KETGammaN[t_] := -gamma*
Im[s + I*p - (s + I*p + gamma) Exp[-I*omega*t]]

PP[t_] :=
Table[Szorz[BRAalfaP[t], KETalfaP[t]]*
Exp[I*(KETGammaP[t] - BRAGammaP[t])] , {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

NN[t_] :=
Table[Szorz[BRAalfaN[t], KETalfaN[t]]*
Exp[I*(KETGammaN[t] - BRAGammaN[t])] , {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

PN[t_] :=
Table[Szorz[BRAalfaP[t], KETalfaN[t]]*
Exp[I*(KETGammaN[t] - BRAGammaP[t])] , {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

NP[t_] :=
Table[Szorz[BRAalfaN[t], KETalfaP[t]]*
Exp[I*(KETGammaP[t] - BRAGammaN[t])] , {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

PtoN[t_] := Inverse[PP[t]].PN[t];
NtoP[t_] := Inverse[NN[t]].NP[t]


And finally, I define the matrix through ArrayFlatten. I can't copy&paste it here, so I give it as a picture:

My question is, how could I make this algorithm run in a reasonably short time with the "size" parameter being larger than 0?

• I can not understand from the figure how the matrix is arranged. This is true m0 = SparseArray[{}, {(1 + 2 size)^2, (1 + 2 size)^2}]; Evomatrix[t_] := ArrayFlatten[{{m0, PtoN[t]}, {NtoP[t], m0}}];? Commented Feb 17, 2020 at 23:31
• Yes, I think that is the correct formula, sorry for the ambiguity above. Commented Feb 17, 2020 at 23:42

We can use the fomal solution of the matrix equation and the numerical integration scheme of 1-8 order. First, write a simple code to compare with NDSolve. This piece of code works for size = 0,1,2... (note, I changed the definition of the matrix Evomatrix[t]):

size = 0; n = 2 (1 + 2*size)^2;

initialCond = X == X0 /. t -> 0;
gamma = -0.005;
omega = 1;
eltolas = 200;
Szorz[elso_, masod_] :=
Exp[-0.5 (Abs[elso])^2 - 0.5 (Abs[masod])^2 + Conjugate[elso]*masod]
l := Sqrt[Pi]*(Mod[y - 1, 1 + 2*size] - size + Re[eltolas]);
m := Sqrt[Pi]*Floor[(y - 1)/(1 + 2*size) - size + Im[eltolas]];
s := Sqrt[Pi]*(Mod[x - 1, 1 + 2*size] - size + Re[eltolas]);
p := Sqrt[Pi]*Floor[(x - 1)/(1 + 2*size) - size + Im[eltolas]];
BRAalfaP[t_] := gamma + (l + I*m - gamma) Exp[-I*omega*t];
KETalfaP[t_] := gamma + (s + I*p - gamma) Exp[-I*omega*t];
BRAalfaN[t_] := -gamma + (l + I*m + gamma) Exp[-I*omega*t];
KETalfaN[t_] := -gamma + (s + I*p + gamma) Exp[-I*omega*t];
BRAGammaP[t_] := gamma*Im[l + I*m - (l + I*m - gamma) Exp[-I*omega*t]];
KETGammaP[t_] := gamma*Im[s + I*p - (s + I*p - gamma) Exp[-I*omega*t]];
BRAGammaN[t_] := -gamma*
Im[l + I*m - (l + I*m + gamma) Exp[-I*omega*t]];

KETGammaN[t_] := -gamma*Im[s + I*p - (s + I*p + gamma) Exp[-I*omega*t]]

PP[t_] :=
Table[Szorz[BRAalfaP[t], KETalfaP[t]]*
Exp[I*(KETGammaP[t] - BRAGammaP[t])], {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

NN[t_] :=
Table[Szorz[BRAalfaN[t], KETalfaN[t]]*
Exp[I*(KETGammaN[t] - BRAGammaN[t])], {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

PN[t_] :=
Table[Szorz[BRAalfaP[t], KETalfaN[t]]*
Exp[I*(KETGammaN[t] - BRAGammaP[t])], {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

NP[t_] :=
Table[Szorz[BRAalfaN[t], KETalfaP[t]]*
Exp[I*(KETGammaP[t] - BRAGammaN[t])], {x, (1 +
2 size)^2}, {y, (1 + 2 size)^2}];

PtoN[t_] := Inverse[PP[t]].PN[t];
NtoP[t_] := Inverse[NN[t]].NP[t];
m0 = SparseArray[{}, {(1 + 2 size)^2, (1 + 2 size)^2}];
Evomatrix[t_] := -0.5 I*ArrayFlatten[{{m0, PtoN[t]}, {NtoP[t], m0}}];
Evomatrix[2] // MatrixForm

tau = .001; X = Table[x[i][t], {t, 0, 1, tau}, {i, n}];
Table[X[[1, i]] = 0, {i, n}];
X[[1, Round[n/2]]] = 1;
T = Table[t, {t, 0, 1, tau}];
Do[X[[k]] = MatrixExp[tau Evomatrix[tau (k - 1/2)], X[[k - 1]]];, {k, 2, Length[X]}]


Compare the solution with NDSolve[]:

    X1[t_] = Table[xi[i, t], {i, 1, n}];
X0 = Table[0, {n}];
X0[[Round[n/2]]] = 1;
diffEq = Table[X1'[t][[i]] == (Evomatrix[t].X1[t])[[i]], {i, 1, n}];
initialCond = Table[X1[0][[i]] == X0[[i]], {i, 1, n}];
sol = NDSolve[{Join[diffEq, initialCond]}, X1[t], {t, 0, 1},
Method -> "BDF", MaxStepSize -> tau];
Table[Show[Plot[Abs[xi[i, t] /. sol], {t, 0, 1}, PlotRange -> All],
ListPlot[Transpose[{T, Abs[Chop[Flatten[X[[All, i]]]]]}],
PlotStyle -> Orange]], {i, 2}]


We see that there are slight differences, but the code can be used for size = 1:

For size=2 and tau=.01 it takes 24 sec on my ASUS

• Thank you! It seems to be working much more quickly with your algorithm. Commented Feb 18, 2020 at 18:52
• @Gomboc You're welcome. What is the maximum system size? Commented Feb 18, 2020 at 20:58