# How to speed up Solve with increasing number of equations

I am numerically solving systems with n number of equations (e.g. like an ecological system with n species), and as I increase the number of equations Solve (or NSolve) grinds to a halt. I suspect it might be how I am building the model because if I code versions where the variables don't use indices, i.e. xi not x[i], then it produces solutions reasonably quickly. But the idea is, of course, to make systems with increasing dimensionality (i.e. various numbers of equations).

Here is an example

n = 2.; (*number of x equations*)
listVs = Flatten[{y, Table[x[i], {i, 1, n}], z}] (*list of variables*)

(*generate equations*)
eqny = a - b y - Sum[c[i] x[i] y, {i, 1, n}];
listeqnx = Table[c[i] x[i] y - d[i] x[i]^2 - (e[i] x[i] z)/(f + z), {i, 1, n}];
eqnz = g Sum[x[i], {i, 1, n}] - h z;

listEqns = Flatten[{eqny, listeqnx, eqnz}];
listEqns0 = Map[# == 0 &, listEqns]; (*set to zero for solver*)
% // MatrixForm

pars1 = {a -> 10.^3, b -> 0.05, f -> 10., g -> 0.1, h -> 0.05 };
pars2 = Flatten[{
Table[c[i] -> First[RandomVariate[UniformDistribution[{0., 1.}], 1] ], {i, 1, n}],
Table[d[i] -> First[RandomVariate[UniformDistribution[{0., 1.}], 1] ], {i, 1, n}],
Table[e[i] -> First[RandomVariate[UniformDistribution[{0., 1.}], 1] ], {i, 1, n}]}];

equils = Solve[listEqns0 /. pars1 /. pars2, listVs, Reals,
WorkingPrecision -> MachinePrecision]; // ByteCount // AbsoluteTiming

eqsM2 = Map[Flatten, listVs /. equils] // MatrixForm (*easier to visualize*)


e.g. for n = 2 it takes 0.1 secs, for n = 4, 0.7 secs, yet when n = 5 it gets stuck (on a 16GB, i7CPU@ 3.4Ghz desktop)

Ultimately the aim is to throw this on a cluster (parallelize, etc.) but first I'd like to make sure I can get the solver to run at reasonable speeds for higher dimensional systems. Please any suggestions on how to speed up this code are greatly appreciated.

Note 1: I am only interested in positive, real solutions, however, I have found that adding conditions beyond 'Reals' to the Solve actually slows it down. Not sure if that might help.

Note 2: Maybe I need to use Compile ? I tried feeding it a list {{y,_Reals},{x[1],_Reals},...etc but it didn't work.

• Try with FindRoot
– zhk
Commented Mar 8, 2018 at 9:19
• Or at least NSolve... Commented Mar 8, 2018 at 9:34
• Hi, thanks for your responses. As I mentioned using NSolve does not solve the problem (it also effectively stalls for n >4), and FindRoot only gives one solution if I implement it in this way: FindRoot[listEqns0 /. pars1 /. pars2, Transpose[{listVs, ConstantArray[1, Length[listVs]]}]] Commented Mar 8, 2018 at 12:34
• just curious, why do you put a decimal point on n = 2.; (*number of x equations*) ? Commented Mar 8, 2018 at 18:23
• I was making values all exact for speed, but I must have gotten carried away :) Commented Mar 9, 2018 at 16:26

If you clear denominators you can get solutions for n=5 case using NSolve. For improved accuracy I would recommend using higher precision by setting input precision to say 400. Could do these modifications as below.

eqs = Expand[Numerator[Together[listEqns /. pars1 /. pars2]]] /.
a_Real :> SetPrecision[a, 400];


From here it proceeds easily.

AbsoluteTiming[sols5Big = NSolve[eqs];]

(* Out[71]= {36.436345, Null} *)

In[72]:= Length[sols5Big]

(* Out[72]= 120 *)

In[73]:= Max[Abs[eqs /. sols5Big]]

(* Out[73]= 2.03*10^-314 *)


If instead one uses machine precision then there are 218 solutions. My suspicion is some arise from an overestimate of multiplicity. The residuals are not bad so at least they are likely to be reasonable.

In case anyone would like to know, here is how I answered my question. I used FindRoot as suggested and looped over many different initial conditions to get several fixed points.

Here is my amended code:

n = 2;(*number of x equations*)
listVs = Flatten[{y, Table[x[i], {i, 1, n}], z}] ;(*list of variables*)
listeqnx = Table[c[i] x[i] y - d[i] x[i]^2 - (e[i] x[i] z)/(f + z), {i, 1, n}];
eqnz = g Sum[x[i], {i, 1, n}] - h z;
listEqns = Flatten[{eqny, listeqnx, eqnz}];
listEqns0 = Map[# == 0 &, listEqns];
jac = Outer[D, listEqns, listVs]; (*find Jacobian, needed for FindRoot*)
(*paramter set*)
pars1 = {a -> 10.^3, b -> 0.05, f -> 10., g -> 0.1, h -> 0.05};
pars2 = Flatten[{Table[c[i] -> First[RandomVariate[UniformDistribution[{0., 1.}], 1]], {i, 1, n}],
Table[d[i] -> First[RandomVariate[UniformDistribution[{0., 1.}], 1]], {i, 1,n}],
Table[e[i] -> First[RandomVariate[UniformDistribution[{0., 1.}], 1]], {i, 1, n}]}];

eqsList = {}; (*for storing equilibria*)
jaci = jac /. pars1 /. pars2; (*for findroot*)
listEqns0b = listEqns0 /. pars1 /. pars2; (*for findroot*)
nfxpts = 500;(*number of tries to find fixedpoints*)
varmax = 1000; (*for uniform distribution to find variable initial guesses [0,1000] *)

Do[
guessgen = Table[First[RandomVariate[UniformDistribution[{0., varmax}], i]],{i, 1, n + 2}];
initg = Transpose[{listVs, guessgen}];
equili = FindRoot[listEqns0b, initg, Jacobian -> jaci, WorkingPrecision -> MachinePrecision] // Quiet;
eqsList = Append[eqsList, equili],
{h, 1, nfxpts}]

eqsM0 = Map[Flatten, listVs /. eqsList];
eqsM = DeleteDuplicates[Round[eqsM0, 0.00000001]]; (*hack for removing duplicates*)
eqsM // MatrixForm


There are a couple of hacks along the way but it gets the job done. Here I numerically integrate to show that what I found using FindRoot are indeed fixed points.

vtrans = {y -> y[t], x[1] -> x1[t], x[2] -> x2[t], z -> z[t]};
eqy = eqny /. vtrans;
eqx1 = listeqnx[[1]] /. vtrans;
eqx2 = listeqnx[[2]] /. vtrans;
eqz = eqnz /. vtrans;
sys1 = {y'[t] == eqy, x1'[t] == eqx1, x2'[t] == eqx2, z'[t] == eqz,
y[0] == eqsM[[1]][[1]], x1[0] == eqsM[[1]][[2]], x2[0] == eqsM[[1]][[3]],
[0] == eqsM[[1]][[4]]}; (*first equil*)
sys2 = {y'[t] == eqy, x1'[t] == eqx1, x2'[t] == eqx2, z'[t] == eqz,
y[0] == eqsM[[2]][[1]], x1[0] == eqsM[[2]][[2]], x2[0] == eqsM[[2]][[3]],
z[0] == eqsM[[2]][[4]]};  (*2nd equil*)

tmax = 5;
sol1 = NDSolve[{sys1 /. pars1 /. pars2}, {y, x1, x2, z}, {t, 0, tmax}];
p1 = Plot[{Evaluate[{y[t]} /. sol1], Evaluate[{x1[t]} /. sol1],
Evaluate[{x2[t]} /. sol1], Evaluate[{z[t]} /. sol1]}, {t, 0, tmax},
PlotLegends -> {y, x1, x2, z}, PlotRange -> {0, All}];

sol2 = NDSolve[{sys2 /. pars1 /. pars2}, {y, x1, x2, z}, {t, 0, tmax}];
p2 = Plot[{Evaluate[{y[t]} /. sol2], Evaluate[{x1[t]} /. sol2],
Evaluate[{x2[t]} /. sol2], Evaluate[{z[t]} /. sol2]}, {t, 0, tmax},
PlotLegends -> {y, x1, x2, z}, PlotRange -> {0, All}];

GraphicsRow[{p1, p2}, ImageSize -> Large]