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

Question

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[14];
n = 10; 
d = 0.5; 
G = RandomGraph[{Round[n], Round[n*(n - 1)*d]}, DirectedEdges -> True]; 
aG = AdjacencyMatrix[G]; 
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[14];
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[#[[1]] -> #[[2]]] & /@Thread[Table[Subscript[#, n], {n, n}] & /@ {\[Alpha], \[Eta], \[Gamma], \[Tau], \[Theta]} -> params]];

initialVals = Array[Subscript[x, #][0] == 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?

Please give me some ideas about improving the calculations.

Note:

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

Answer 1

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[14];

 (*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[14];
 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[#[[1]] -> #[[2]]] & /@ 
    Thread[Table[
         Subscript[#, n], {n, 
          n}] & /@ {\[Alpha], \[Eta], \[Gamma], \[Tau], \[Theta]} -> 
      params]];

 initialVals = Array[Subscript[x, #][0] == 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)"