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.
FindRoot
$\endgroup$NSolve
... $\endgroup$NSolve
does not solve the problem (it also effectively stalls for n >4), andFindRoot
only gives one solution if I implement it in this way:FindRoot[listEqns0 /. pars1 /. pars2, Transpose[{listVs, ConstantArray[1, Length[listVs]]}]]
$\endgroup$n = 2.; (*number of x equations*)
? $\endgroup$