How to speed up Solve with increasing number of equations

Question

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.

Answer 1

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.

Answer 2

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]