1
$\begingroup$

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: enter image description here

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

$\endgroup$
2
  • $\begingroup$ 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}}];? $\endgroup$ Feb 17, 2020 at 23:31
  • $\begingroup$ Yes, I think that is the correct formula, sorry for the ambiguity above. $\endgroup$
    – Gomboc
    Feb 17, 2020 at 23:42

1 Answer 1

1
$\begingroup$

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}]

Figure 1

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

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.