# Reduce and Solve are taking forever

I am running the following script to solve a system of equations. However, it has been running for 12 hours by now with no output. Is this typical? What is happening? Is there anything wrong with how I set up the problem or something to make it run faster? How can I check the progress?

eq1 = {b1, b2} . {1, 1};
eq2 = {b1, b2} . {{a11 + a12}, {a21 + a22}};
eq3 = {b1, b2} . {{(a11 + a12)^2}, {(a21 + a22)^2}};
eq4 = {b1, b2} . {{a11*(a11 + a12) + a12 (a21 + a22)}, {a21*(a11+a12) + a22 (a21 + a22)}};

system = {eq1 == 1, eq2 == 1/2, eq3 == 1/3, eq4 == 1/6,{{r11,r12},{r21,r22}} ==
Transpose[{{r11,r12},{r21,r22}}], Eigenvalues[{{2 a11 r11 + a21 r12 + a21 r21, a12
r11 + a11 r12 + a22 r12 + a21 r22, b1},{a12 r11 + a11 r21 + a22 r21 +a21 r22, a12 r12
+ a12 r21 + 2 a22 r22, b2},{b1, b2, 1}}]>=0}

Reduce[system, {b1, b2, a11, a12, a21, a22, r11, r12, r21, r22}]
Solve[system, {b1, b2, a11, a12, a21, a22, r11, r12, r21, r22}]

• You have at least one typo: 4==1/6 should be eq4==1/6.
– JimB
Commented Apr 15, 2022 at 18:00
• @JimB Thanks that's not in the actual program I ran though Commented Apr 15, 2022 at 18:09
• I think you will need to reformulate your problem, reducing the number of variables involved. At first glance, I doubt that the equation would give unique values for your r variables (but I could be wrong). Commented Apr 15, 2022 at 20:09
• @mikado Are you suggesting that FindInstance might give a solution? I tired it and it took forever too Commented Apr 15, 2022 at 20:15
• I don't think that FindInstance is likely to help either. I don't have a definite suggestion, but I would be analysing this as a matrix/vector problem, rather than decomposing it into matrix/vector elements. I'm not sure that Mathematica can help much there, but I haven't looked very closely. Commented Apr 15, 2022 at 22:08

This is just an extended comment that suggests that either this is just too complex for Mathematica or that there really isn't much of a "reduction" in the system to be had.

Consider solving the equalities first:

sol = Solve[{eq1 == 1, eq2 == 1/2, eq3 == 1/3, eq4 == 1/6}, {b1, b2, a11, a12}] // Flatten


Now construct the matrix (with r12 == r21  as that is essentially what your Transpose statement does, i.e., the 2x2 "r" matrix is symmetric):

m = {{2 a11 r11 + a21 r12 + a21 r21, a12 r11 + a11 r12 + a22 r12 + a21 r22, b1},
{a12 r11 + a11 r21 + a22 r21 + a21 r22, a12 r12 + a12 r21 + 2 a22 r22, b2},
{b1, b2, 1}} /. r21 -> r12 /. sol // FullSimplify;


Now the system is

system = Eigenvalues[m] > 0


Looking at that result, I'm not so sure any more "reduction" can be achieved.

You can take the extented hints of @JimB, proceed, and get analytical solutions for at least some special cases.

You get different solutions sets for different subsets of variables. (Suppose, you want only real variables)

eq1 = {b1, b2}.{1, 1};
eq2 = {b1, b2}.{{a11 + a12}, {a21 + a22}};
eq3 = {b1, b2}.{{(a11 + a12)^2}, {(a21 + a22)^2}};
eq4 = {b1,
b2}.{{a11*(a11 + a12) + a12 (a21 + a22)}, {a21*(a11 + a12) +
a22 (a21 + a22)}};

vars14 = Variables[Flatten[{eq1, eq2, eq3, eq4}]]

ss = Subsets[vars14, {4}]


Get very general solution as @JimB or some specific solutions.

sol = Solve[{eq1 == 1, eq2 == 1/2, eq3 == 1/3, eq4 == 1/6}, {b1, b2,
a11, a12}, Reals, MaxExtraConditions -> All] // Simplify

(*   {{b1 -> ConditionalExpression[
1 - 1/(4 (1 + 3 a21^2 - 3 a22 + 3 a22^2 + a21 (-3 + 6 a22))),
a21 + a22 != 1/2],
b2 -> ConditionalExpression[1/(
4 (1 + 3 a21^2 - 3 a22 + 3 a22^2 + a21 (-3 + 6 a22))),
a21 + a22 != 1/2],
a11 -> ConditionalExpression[(
1 + 6 a21^2 - 5 a22 + 6 a22^2 + 6 a21 (-1 + 2 a22))/(
3 (-1 + 2 a21 + 2 a22)^2), a21 + a22 != 1/2],
a12 -> ConditionalExpression[-((-1 + a21 + 2 a22)/(
3 (-1 + 2 a21 + 2 a22)^2)), a21 + a22 != 1/2]}}   *)

Manipulate[
solj = Solve[{eq1 == 1, eq2 == 1/2, eq3 == 1/3, eq4 == 1/6}, ss[[j]],
Reals, MaxExtraConditions -> All] // Simplify, {{j, 10}, 1,
Length[ss], 1, Setter}]


Regard the first singular solution, get simple matrix m, denpending on only 3 variables.

solj[[1]]

solj1 = Equal @@@ (solj[[1]]) // Solve // Flatten

(*   {a11 -> 1/6, a12 -> 0, a21 -> 1/2, a22 -> 1/4, b1 -> 3/7, b2 -> 4/7}   *)

(m = {{2 a11 r11 + a21 r12 + a21 r21,
a12 r11 + a11 r12 + a22 r12 + a21 r22,
b1}, {a12 r11 + a11 r21 + a22 r21 + a21 r22,
a12 r12 + a12 r21 + 2 a22 r22, b2}, {b1, b2, 1}} /.
r21 -> r12 /. solj1 // FullSimplify) // MatrixForm

Variables[Flatten@m]

ev = Eigenvalues[m]

(*   {Root[768 r11 + 864 r12 + 1225 r12^2 - 1080 r22 - 1176 r11 r22 -
588 r12 r22 +
1764 r22^2 + (-3600 + 2352 r11 + 7056 r12 - 1225 r12^2 +
3528 r22 + 1176 r11 r22 + 588 r12 r22 -
1764 r22^2) #1 + (-7056 - 2352 r11 - 7056 r12 -
3528 r22) #1^2 + 7056 #1^3 &, 1],
Root[768 r11 + 864 r12 + 1225 r12^2 - 1080 r22 - 1176 r11 r22 -
588 r12 r22 +
1764 r22^2 + (-3600 + 2352 r11 + 7056 r12 - 1225 r12^2 +
3528 r22 + 1176 r11 r22 + 588 r12 r22 -
1764 r22^2) #1 + (-7056 - 2352 r11 - 7056 r12 -
3528 r22) #1^2 + 7056 #1^3 &, 2],
Root[768 r11 + 864 r12 + 1225 r12^2 - 1080 r22 - 1176 r11 r22 -
588 r12 r22 +
1764 r22^2 + (-3600 + 2352 r11 + 7056 r12 - 1225 r12^2 +
3528 r22 + 1176 r11 r22 + 588 r12 r22 -
1764 r22^2) #1 + (-7056 - 2352 r11 - 7056 r12 -
3528 r22) #1^2 + 7056 #1^3 &, 3]}   *)

RegionPlot3D[
And @@ Thread[ev >= 0], {r11, 0, 10}, {r12, -5, 5}, {r22, 0, 10},
PlotPoints -> 30]


ev /. Thread[{r11, r12, r22} -> {10, 1, 2}] // N

(*   {0.253483, 1.13394, 4.94591}   *)


It is sufficient to regard ev[1], which can be solved anylytivally.

RegionPlot3D[ev[[1]] >= 0, {r11, 0, 10}, {r12, -5, 5}, {r22, 0, 10},
PlotPoints -> 30]

red31 = Reduce[
ev[[1]] == 0 && 0 < r11 < 10 && -5 < r12 < 5 && 0 < r22 < 10, {r11,
r12, r22}]

(*   (r11 == 279/245 && r12 == -(48/245) &&
r22 == 32/49) || (279/245 < r11 <
10 && ((r12 == 1/147 (27 - 49 r11) &&
r22 == 1/294 (90 + 98 r11 + 49 r12) -
1/147 Sqrt[
2025 - 4998 r11 + 2401 r11^2 - 8379 r12 + 2401 r11 r12 -
14406 r12^2]) || (1/147 (27 - 49 r11) < r12 <
1/98 (-75 + 49 r11) && (r22 ==
1/294 (90 + 98 r11 + 49 r12) -
1/147 Sqrt[
2025 - 4998 r11 + 2401 r11^2 - 8379 r12 + 2401 r11 r12 -
14406 r12^2] ||
r22 == 1/294 (90 + 98 r11 + 49 r12) +
1/147 Sqrt[
2025 - 4998 r11 + 2401 r11^2 - 8379 r12 + 2401 r11 r12 -
14406 r12^2])) || (r12 == 1/98 (-75 + 49 r11) &&
r22 == 1/294 (90 + 98 r11 + 49 r12) -
1/147 Sqrt[
2025 - 4998 r11 + 2401 r11^2 - 8379 r12 + 2401 r11 r12 -
14406 r12^2])))   *)


This was only the example for a singular solution.

Getting general solutions seems to be a monster-job, may be impossible with analytical tools.