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

Question

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

Answer 1

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[1] *)
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.

Answer 2

`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.