# My Program using NSolve and FindInstance has been running for days without a solution

I have a system of 5 non-linear equations and 5 unknowns (p1,p2,p3,p4,p5 bellow), and a few inequalities. I am trying to understand whether the solution is unique. For this purpose, I tried to run FindInstance asking for 2 solutions, but the system could not provide an answer even after 2 days. I also tried NSolve and had the same problem. (I also tried FindInstance asking for 1 solution and that didn’t work either, even though, analytically I know that a solution must exist). I would really appreciate your help. Here is the code:

SetOptions[EvaluationNotebook[], CellContext -> Notebook]
ClearAll["Global*"];

k1 = RandomInteger[{1, 100}];
k2 = RandomInteger[{1, 100}];
k3 =  RandomInteger[{1, 100}];
k4 = RandomInteger[{1, 100}];
k5 = RandomInteger[{1, 100}];
k12 = RandomInteger[{1, 100}];
k13 = RandomInteger[{1, 100}];
k14 = RandomInteger[{1, 100}];
k15 = RandomInteger[{1, 100}];
k23 = RandomInteger[{1, 100}];
k24 = RandomInteger[{1, 100}];
k25 = RandomInteger[{1, 100}];
k34 = RandomInteger[{1, 100}];
k35 = RandomInteger[{1, 100}];
k45 = RandomInteger[{1, 100}];
k123 = RandomInteger[{1, 100}];
k124 = RandomInteger[{1, 100}];
k125 = RandomInteger[{1, 100}];
k134 = RandomInteger[{1, 100}];
k135 = RandomInteger[{1, 100}];
k145 = RandomInteger[{1, 100}];
k234 = RandomInteger[{1, 100}];
k235 = RandomInteger[{1, 100}];
k245 = RandomInteger[{1, 100}];
k345 = RandomInteger[{1, 100}];
k1234 = RandomInteger[{1, 100}];
k1235 = RandomInteger[{1, 100}];
k1245 = RandomInteger[{1, 100}];
k1345 = RandomInteger[{1, 100}];
k2345 = RandomInteger[{1, 100}];
k = {k1, k2, k3, k4, k5, k12, k13, k14, k15, k23, k24, k25, k34, k35,
k45, k123, k124, k125, k134, k135, k145, k234, k235, k245, k345,
k1234, k1235, k1245, k1345, k2345};
n = k1 + k2 + k3 + k4 + k5 + k12 + k13 + k14 + k15 + k23 + k24 + k25 +
k34 + k35 + k45 + k123 + k124 + k125 + k134 + k135 + k145 +
k234 + k235 + k245 + k345 + k1234 + k1235 + k1245 + k1345 + k2345;
kk = N[{k1/n, k2/n, k3/n, k4/n, k5/n, k12/n, k13/n, k14/n, k15/n,
k23/n, k24/n, k25/n, k34/n, k35/n, k45/n, k123/n, k124/n, k125/n,
k134/n, k135/n, k145/n, k234/n, k235/n, k245/n, k345/n, k1234/n,
k1235/n, k1245/n, k1345/n, k2345/n}];

A = { k1/p1 + k12/(p1 + p2) + k13 /(p1 + p3) + k14/(p1 + p4) +
k15/(p1 + p5) + k123/(p1 + p2 + p3) + k124/(p1 + p2 + p4) +
k125/(p1 + p2 + p5) + k134 /(p1 + p3 + p4) + k135/(p1 + p3 + p5) +
k145/(p1 + p4 + p5) + k1234/(p1 + p2 + p3 + p4) +
k1235/(p1 + p2 + p3 + p5) + k1245/(p1 + p2 + p4 + p5) +
k1345 /(p1 + p3 + p4 + p5) == n,
k2/p2 + k12/(p1 + p2) + k23 /(p2 + p3) + k24 /(p2 + p4) +
k25/(p2 + p5) + k123 /(p1 + p2 + p3) + k124 /(p1 + p2 + p4) +
k125/(p1 + p2 + p5) + k234/(p2 + p3 + p4) + k235 /(p2 + p3 + p5) +
k245/(p2 + p4 + p5) + k1234/(p1 + p2 + p3 + p4) +
k1235/(p1 + p2 + p3 + p5) + k1245/(p1 + p2 + p4 + p5) +
k2345/(p2 + p3 + p4 + p5) == n,
k3/p3 + k13 /(p1 + p3) + k23/(p2 + p3) + k34/(p3 + p4) +
k35 /(p3 + p5) + k123/(p1 + p2 + p3) + k134/(p1 + p3 + p4) +
k135/(p1 + p3 + p5) + k234/(p2 + p3 + p4) +
k235/(p2 + p3 + p5) + k345/(p3 + p4 + p5) +
k1234/(p1 + p2 + p3 + p4) + k1235/(p1 + p2 + p3 + p5) +
k1345/(p1 + p3 + p4 + p5) + k2345/(p2 + p3 + p4 + p5) == n,
k4/p4 + k14/(p1 + p4) + k24/(p2 + p4) + k34/(p3 + p4) +
k45/(p4 + p5) + k124 /(p1 + p2 + p4) + k134/(p1 + p3 + p4) +
k145/(p1 + p4 + p5) + k234 /(p2 + p3 + p4) + k245/(p2 + p4 + p5) +
k345 /(p3 + p4 + p5) + k1234/(p1 + p2 + p3 + p4) +
k1245/(p1 + p2 + p4 + p5) + k1345/(p1 + p3 + p4 + p5) +
k2345/(p2 + p3 + p4 + p5) == n ,
k5/p5 + k15/(p1 + p5) + k25/(p2 + p5) + k35/(p3 + p5) +
k45/(p4 + p5) + k125 /(p1 + p2 + p5) + k135/(p1 + p3 + p5) +
k145/(p1 + p4 + p5) + k235/(p2 + p3 + p5) + k245/(p2 + p4 + p5) +
k345 /(p3 + p4 + p5) + k1235/(p1 + p2 + p3 + p5) +
k1245/(p1 + p2 + p4 + p5) + k1345/(p1 + p3 + p4 + p5) +
k2345/(p2 + p3 + p4 + p5) == n , p1 >= 0, p1 <= 1, p2 >= 0,
p2 <= 1, p3 >= 0, p3 <= 1, p4 >= 0, p4 <= 1, p5 >= 0, p5 <= 1,
p1 + p2 + p3 + p4 + p5 == 1}

Print["===================="];
f = FindInstance[A, {p1, p2, p3, p4, p5}, Reals, 2]
p = NSolve[A, {p1, p2, p3, p4, p5}, Reals]
dimp = Dimensions[p];

• 1) Is there a particular random seed needed or will it always have solutions for any setting of the k? 2) I could show you how to get some close solutions by minimization - how precise do you want the solutions? 3) How did you construct the equations to begin with and how do you know a solution exists? Aug 8, 2020 at 12:33
• 1) It is possible using a fixed point theorem to show that a solution exists for any of k-s that are positive ( k1,…,k5 must be strictly positive). 2) I don’t really care about the solution itself, I just want to know if it is unique. Thanks for your help Aug 8, 2020 at 12:58
• In that case it looks completely intractable to me. If you look at your first equation and put it in ==0 form, Together[A[] /. Equal -> Subtract] // Numerator it's equivalent to finding the roots of a massive multivariate polynomial. Add in the other equations and inequalities that must be satisfied and this problem is infinitely more difficult. It might be possible to show that two 'close' numerical results exist, though numerical results tell us nothing about uniqueness. Again, I think you should show how you constructed these equations because they are too complex in this form. Aug 8, 2020 at 13:30
• 1) What do you mean by "show how you constructed the equations? In a simpler version with only 4 equations and 4 unknowns, I was able to get a solution. Unfortunately, it does not work with 5 2) How can I get close solutions by minimization? Perhaps that could help Aug 8, 2020 at 13:45
• The fixed point theorem only guarantees that the lhs of each equation has a fixed point - it does not mean that that this fixed point needs to be the same for all equations. As far as I can tell there are no guarantees a solution exists. Aug 8, 2020 at 14:20

FindRoot give a close solution, and the evidence suggests there is no other solution in the constrained search space:

(* initialized the k's and A with SeedRandom *)
vars = {p1, p2, p3, p4, p5};
FindRoot[Most@ Cases[A, _Equal],
{{p1, 1/5}, {p2, 1/5}, {p3, 1/5}, {p4, 1/5}, {p5, 1/5}}]

(*
{p1 -> 0.144271, p2 -> 0.246401, p3 -> 0.25002, p4 -> 0.219996,
p5 -> 0.142424}
*)


It is close, but not very close, to satisfying the 6th constraint p1 + p2 + p3 + p4 + p5 == 1. There is naturally some doubt whether six equations in five unknowns even has a solution.

vars /. % // Total

(*  1.01634  *)


Randomly choosing starting points in the unit cube all led to the same solution. Here's an attempt to find a better solution:

Do[
If[0.985 <
Total[vars /.
FindRoot[Most@Cases[A, _Equal],
Transpose@{vars, #, 0 vars + \$MachineEpsilon, 0 vars + 1.}]] <
1.015,
Return[#, Do]
] &@RandomReal[1, 5],
{100}]


It returns nothing in all my trials. If it finds a different, better solution, it will return the starting points. One thousand iterations takes less than 2.5 seconds.

• Thanks, this is very helpful Aug 9, 2020 at 7:12
A = {p1 >= 0, p1 <= 1, p2 >= 0, p2 <= 1, p3 >= 0, p3 <= 1, p4 >= 0,   p4 <= 1, p5 >= 0, p5 <= 1, p1 + p2 + p3 + p4 + p5 == 1}
p = Reduce[A, {p1, p2, p3, p4, p5}, Reals]
dimp = Dimensions[p]


((p1 == 0 && ((p2 ==
0 && ((p3 == 0 && 0 <= p4 <= 1) || (0 < p3 < 1 &&
0 <= p4 <= 1 - p3) || (p3 == 1 && p4 == 0))) || (0 < p2 <
1 && ((0 <= p3 < 1 - p2 &&
0 <= p4 <= 1 - p2 - p3) || (p3 == 1 - p2 &&
p4 == 0))) || (p2 == 1 && p3 == 0 && p4 == 0))) || (0 <
p1 <= 1 && ((0 <= p2 <
1 - p1 && ((0 <= p3 < 1 - p1 - p2 &&
0 <= p4 <= 1 - p1 - p2 - p3) || (p3 == 1 - p1 - p2 &&
p4 == 0))) || (p2 == 1 - p1 && p3 == 0 && p4 == 0)))) &&
p5 == 1 - p1 - p2 - p3 - p4


p1==0 implies the is not satisfiable.

This takes 0.8 s on my computer.

The other part of A is not solvable with Mathematica. Neither Reduce, Solve, NSolve, FindRoot, ... will make a solution. Five equations require five parameters for the solution. No chance.

It is on occasion meaningful to find solutions of subsets of the equation first. The most concrete clue has Roots itself. Roots can solve equations of such high order but that must not be. Mathematica is finished by order 3 or 4 if it gets complicated. For FindInstance higher order polynomials are solvable but under the restriction that only up to 3 or 4 parameters are in the set of equations and that these are not overdetermined.

There by be chance that certain combination of the k's can be solved, but at random this is low, very low.

Many of these higher dimensional problems are ill-conditioned. If built-ins run suspiciously long terminate the run. Many conditions lead to even more ill-conditions. Try to reduced the dimensionality first and then run the built-ins in lower dimensions for speed and power and enhanced success chances.

I personally prefer to work with Reduce only. A is in the form more suitable for Reduce. FindInstance prefers in this cases the logical conjunctions. If this was not overdetermined the RandomSeeding option of FindInstance may give an time reduce.

A really good starting point is calculating so called degrees of freedom. If there are 5 variable only 5 equations can be satisfied with each additional equation limits without an extra degree of freedom the chances to get a good solution or a solution at all. The mathematical methodology is called lagrangian multipliers.

Mathematica is not strong in theory but in allowing to solve problems in applied mathematics. Theory of solvability is a step child of most mathematica authors and professors. It is thought to be an advanced chapter in mathematics. That is compensated in some generations of students by passing them criteria to consider if a problem has a solution. Mathematica nowadays includes more and more tools to judge whether solutions exist. FindInstance, Solve or Reduce do not.

This knowledge is to decide when Solve is suitable. Getting always instructions and answer is boring. Making discoveries on oneself behalve is great. So tried a walkthrough the knowledge conserved in mathematica.stackexchange for example on th comparison between Solve and Reduce: Search for Solve and Reduce on mathematia.stackexchange.com.

The is the article that got most votes: what-is-the-difference-between-reduce-and-solve.

The main idea to keep in mind is:

Solve[a*x == 0, x]
(* Out: {{x -> 0}} *)

Reduce[a*x == 0, x]
(* Out: a == 0 || x == 0 *)


So Reduce is ad hoc much more flexible and versatile. Indeed Reduce is a so-called master built-in aggregating all the other built-ins of Mathematica into one big picture scope built-in that can be used to benchmark against sets of problems or other software comparable to such challenging sets. So Backsubstitution and check for the validity. It knows LogicalExpand. The advertising sentence from the Mathematica documentation is:

When expr involves only polynomial equations and inequalities over real or complex domains, then Reduce can always in principle solve directly for all the Subscript[x, i].

So this, on the other hand, states the limits of Reduce to cope with polynomials and inequalities on Reals or Complexes and subsets.

It solves but it does not check for solvability!

When expr involves only polynomial conditions, Reduce[expr,vars,Reals] gives a cylindrical algebraic decomposition of expr`.

State about the methodology implemented and applied. A cylindrical algebraic decomposition is a very sophisticated and advanced methodology that is usually not taught to students just in specialized courses for phd students and more advanced one. This methodology has a branch dealing with the solvability of problems.

• So you think Reduce will work better than Nsolve and FindInstance? Aug 8, 2020 at 16:03
• This is not really addressing the question that was asked. I realize there are suggestions, but without actually trying them it might be better to offer them in comments. Aug 8, 2020 at 22:10
• @gabi-gayer: No, no, no. I think Reduce is nicer in the postprocessing of the output and better readable with that. All built-ins to find solution have their methods and accompany together to give better performances in certain applications. Sep 20, 2020 at 17:54