# Any idea? why NDSolve is very slow for a very small digraph

I have developed (of course with a lot of help from experts in this forum) the following code to speed up the calculations for large digraphs. However, I ended up with much slower calculations compared to Do loop formulation.

ClearAll[n, d, G, aG, aGsym, Xt, Xt1, alfa, eta, gama, tao, teta, one,
onePlusEta, Aterm, Bterm, Cterm, Dterm, eqs, modelFull];

(* generating a digraph *)
SeedRandom;
n = 10;
d = 0.5;
G = RandomGraph[{Round[n], Round[n*(n - 1)*d]}, DirectedEdges -> True];
aGsym = BitAnd[aG, Transpose[aG]];

(* preliminaries of matrix inputs for the model *)
Xt = DiagonalMatrix[Array[Subscript[x, #][t] &, n]];
Xt1 = Array[Subscript[x, #]'[t] &, n];
alfa = Transpose@Array[Subscript[\[Alpha], #] &, {n, n}];
eta = Array[Subscript[\[Eta], #] &, n];
gama = Transpose[Array[Subscript[\[Gamma], #] &, {n, n}]]*aG;
tao = Array[Subscript[\[Tau], #] &, {n, n}];
teta = Transpose[Array[Subscript[\[Theta], #] &, {n, n}]*aG];
one = Array[1 &, n];

(* building the model *)
onePlusEta = DiagonalMatrix[one + eta]*Xt; Aterm = (gama + teta).Xt;
Bterm = Xt.((alfa*tao)*aG);
Cterm = ((alfa*alfa)*aGsym)*(Xt.Array[Diagonal[Xt] &, {n}])*aGsym;
Dterm = ((gama*gama)*aGsym)*
Transpose[Array[(Subscript[x, #][t])^2 &, {n, n}]*aGsym];
eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm]//FullSimplify;
modelFull = Thread[Xt1 == eqs - Array[Subscript[x, #][t] &, n]];

(* parameter values for numeric solutions *)
ClearAll[params, parameters, initialVals, sol];
SeedRandom;
params = {abs, random, absCost, trsCost, trs} = {
RandomReal[{0.01, 1}, n], RandomReal[{-0.5, 0.5}, n],
RandomReal[{0.01, 0.5}, n], RandomReal[{0.01, 1}, n],
RandomReal[{0.01, 0.5}, n]
};
parameters = Flatten[Thread[#[] -> #[]] & /@Thread[Table[Subscript[#, n], {n, n}] & /@ {\[Alpha], \[Eta], \[Gamma], \[Tau], \[Theta]} -> params]];

initialVals = Array[Subscript[x, #] == RandomReal[{0.01, 1}] &, n];
sol = NDSolve[{modelFull /. parameters, initialVals} // Flatten,
Array[Subscript[x, #][t] &, n], {t, 50},
Method-> "StiffnessSwitching"] // Flatten;

Plot[Evaluate[Array[Subscript[x, #][t] &, n] /. sol], {t, 0, 50},
PlotRange -> All]


Is there any specific reason for the slow process? How can I really improve the performance of the code with larger digraphs?

OR Am I on the WRONG track for speeding up the calculations?

Note:

I can also send the other code using Do and If if somebody is interested in.

• Code works fast. And what is calculated? – Alex Trounev Oct 30 '18 at 21:19
• @AlexTrounev: Thank you very much for checking the speed of the code. The code uses a digraph G to create a system of ODEs and then solve the system sol using the created system. I would appreciate if you run the code with n=50 and n=100 and give me the timing for each one. MyMathematica version is 10.0.2.0 and use Windows 7 (with RAM 4). Also tell me on which machine you run the model. Thank you again. – Tugrul Temel Oct 30 '18 at 21:58
• NDSolve spits out ndsz warning and fails at about t==5.57, tested in v9.0.1 and v11.2. – xzczd Oct 31 '18 at 6:47

I found that the main loss of time and resources occurs while simplifying the equations in line eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm]//FullSimplify;. Moreover, this simplification with a large number of equations leads to errors. I replaced FullSimplify with Simplify[eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm], TimeConstraint -> ts]; with the option TimeConstraint and added to the code a convenient analysis of the influence of the number of equations, the time of symbolic and numerical calculations.

AbsoluteTiming[
ClearAll[n, d, G, aG, aGsym, Xt, Xt1, alfa, eta, gama, tao, teta,
one, onePlusEta, Aterm, Bterm, Cterm, Dterm, eqs, modelFull];
SeedRandom;

(*generating a digraph*)
n = 100; ts = 1; tm = 2; d = 0.5;
G = RandomGraph[{Round[n], Round[n*(n - 1)*d]},
DirectedEdges -> True]; aG = AdjacencyMatrix[G];
aGsym = BitAnd[aG,
Transpose[
aG]];(*picks symmetric cells from "aG"*)(*preliminaries of matrix \
inputs for the model*)
Xt = DiagonalMatrix[
Array[Subscript[x, #][t] &,
n]];(*Subscript[x,i][t] is vertex i at time t*)
Xt1 = Array[Subscript[x, #]'[t] &,
n];(*time derivative of Subscript[x,i][t]*)
alfa = Transpose@Array[Subscript[\[Alpha], #] &, {n, n}];
eta = Array[Subscript[\[Eta], #] &, n];
gama = Transpose[Array[Subscript[\[Gamma], #] &, {n, n}]]*aG;
tao = Array[Subscript[\[Tau], #] &, {n, n}];
teta = Transpose[Array[Subscript[\[Theta], #] &, {n, n}]*aG];
one = Array[1 &, n];
(*building the model*)onePlusEta = DiagonalMatrix[one + eta]*Xt;
Aterm = (gama + teta).Xt;
Bterm = Xt.((alfa*tao)*aG);
Cterm = ((alfa*alfa)*aGsym)*(Xt.Array[Diagonal[Xt] &, {n}])*aGsym;
Dterm = ((gama*gama)*aGsym)*
Transpose[Array[(Subscript[x, #][t])^2 &, {n, n}]*aGsym];
Simplify[eqs = Total[onePlusEta - Aterm + Bterm + Cterm - Dterm],
TimeConstraint -> ts];
modelFull = Thread[Xt1 == eqs - Array[Subscript[x, #][t] &, n]];

(*parameter values for numeric solutions*)
ClearAll[params, parameters, initialVals, sol];
SeedRandom;
params = {abs, random, absCost, trsCost,
trs} = {RandomReal[{0.01, 1}, n], RandomReal[{-0.5, 0.5}, n],
RandomReal[{0.01, 0.5}, n], RandomReal[{0.01, 1}, n],
RandomReal[{0.01, 0.5}, n]};
parameters =
Subscript[#, n], {n,
n}] & /@ {\[Alpha], \[Eta], \[Gamma], \[Tau], \[Theta]} ->
params]];

initialVals = Array[Subscript[x, #] == RandomReal[{0.01, 1}] &, n];
sol = NDSolve[{modelFull /. parameters, initialVals} // Flatten,
Array[Subscript[x, #][t] &, n], {t, tm},
Method -> "StiffnessSwitching"] // Flatten;

Plot[Evaluate[Array[Subscript[x, #][t] &, n] /. sol], {t, 0, tm},
PlotRange -> All,
PlotLabel -> Grid[{{"n", "ts", "tm"}, {n, ts, tm}}, Frame -> All]]] In my version, simplification does not affect the result, but increases the total time by 7-8 times.

\$Version

Out[]= "11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)"

• Thank you so much for the neat analysis and result. With your fix, now the model behaves itself and works smoothly, and takes only 1 minute for n=100. I can live with this speed. – Tugrul Temel Oct 31 '18 at 9:40
• You're welcome! You can do this analysis by simply putting on each suspicious line //AbsoluteTiming – Alex Trounev Oct 31 '18 at 9:55