# Solve is being too slow - any way to solve this system of equations?

I have 3 equations with 3 variables and I want a solution for the 3 of them. Moreover, I'm not really concerned for solutions outside of the range between 0 and 1. I tried using Solve to do this, but it's taking way too long. Any suggestions? Each of the 3 equations needs to be solved for 0. I need to do this several times (say, 20-30) so it would be optimal to have a way to do this that it wouldn't take the whole day and without having mathematica crashing in the middle of it.

Solve[{(-0.8702809089282518 + 1.2424405081409913 Y2^0.5 +
1.7541599570585076 Y2 + 1.2424405081409913 Y2^1. -
1.2424405081409913 Y2^1.5 - 1.1056986938629452 Y2^2 -
1.2424405081409913 Y2^2. + 0.3397298264448023 Y2^2.5 +
0.2218196457326892 Y2^3 + 0.3397298264448023 Y2^3. +
Y1^2 (-0.016997674002504914 + 0.02773304705671855 Y2^0.5 +
0.016997674002504914 Y2 + 0.02773304705671855 Y2^1.) +
Y0^2 (-0.10623546251565573 + 0.17333154410449098 Y2^0.5 +
0.10623546251565573 Y2 + 0.17333154410449098 Y2^1.) +
Y1 (0.2447665056360708 - 0.35498300232599744 Y2^0.5 -
0.3705492932546072 Y2 - 0.35498300232599744 Y2^1. +
0.16085167292896763 Y2^1.5 + 0.1257827876185364 Y2^2 +
0.16085167292896763 Y2^2.) +
Y0 (0.6119162640901772 - 0.8874575058149936 Y2^0.5 -
0.926373233136518 Y2 - 0.8874575058149936 Y2^1. +
0.402129182322419 Y2^1.5 + 0.314456969046341 Y2^2 +
0.402129182322419 Y2^2. +
Y1 (-0.08498837001252461 + 0.13866523528359276 Y2^0.5 +
0.08498837001252461 Y2 +
0.13866523528359276 Y2^1.)))/((1. +
Y2^0.5) (1.7777777777777777 - 0.5555555555555556 Y0 -
0.2222222222222222 Y1 - 1. Y2)^2 Y2^0.5 (-8. + 5. Y0 +
2. Y1 + 1. Y2)) == 0,
1/(3 + Y0^0.5 + 2 Y1^0.5)^2 ((
0.03668478260869564 (1 + Y0^0.5) (3 + Y0^0.5 +
2 Y1^0.5) (3 - Y0 - 2 Y1))/(3 - 1.4375 Y0 - 1.375 Y1 -
0.1875 Y2)^2 + (
0.07336956521739128 (1 + Y1^0.5) (3 + Y0^0.5 +
2 Y1^0.5) (3 - Y0 - 2 Y1))/(3 - 1.4375 Y0 - 1.375 Y1 -
0.1875 Y2)^2 + (
0.16843100189035914 (1. + Y0^0.5) (3. + Y0^0.5 +
2. Y1^0.5) (-1.9393939393939392 + 1. Y0 +
0.7878787878787878 Y1 +
0.15151515151515152 Y2))/((-2.0869565217391304 + 1. Y0 +
0.9565217391304348 Y1 + 0.13043478260869565 Y2) (-1.6 +
1. Y0 + 0.4 Y1 + 0.2 Y2)) - (
0.08421550094517957 (1. + Y0^0.5) (-3. + Y0 +
2. Y1) (-1.9393939393939392 + 1. Y0 +
0.7878787878787878 Y1 + 0.15151515151515152 Y2))/(
Y1^0.5 (-2.0869565217391304 + 1. Y0 +
0.9565217391304348 Y1 + 0.13043478260869565 Y2) (-1.6 +
1. Y0 + 0.4 Y1 + 0.2 Y2)) + (
0.21521739130434786 (3. + Y0^0.5 + 2. Y1^0.5) (-3. + Y0 +
2. Y1) (-1.9393939393939392 + 1. Y0 +
0.7878787878787878 Y1 + 0.15151515151515152 Y2))/(
Y1^0.5 (-2.0869565217391304 + 1. Y0 +
0.9565217391304348 Y1 + 0.13043478260869565 Y2) (-1.6 +
1. Y0 + 0.4 Y1 + 0.2 Y2)) + (
0.693675889328063 (1. + Y1^0.5) (3. + Y0^0.5 +
2. Y1^0.5) (-2.4615384615384612 +
1.2692307692307692 Y0 + 1. Y1 +
0.1923076923076923 Y2))/((-2.1818181818181817 +
1.0454545454545454 Y0 + 1. Y1 +
0.13636363636363635 Y2) (-4. + 2.5 Y0 + 1. Y1 +
0.5 Y2)) - (
0.3468379446640315 (1. + Y1^0.5) (-3. + Y0 +
2. Y1) (-2.4615384615384612 + 1.2692307692307692 Y0 +
1. Y1 + 0.1923076923076923 Y2))/(
Y1^0.5 (-2.1818181818181817 + 1.0454545454545454 Y0 +
1. Y1 + 0.13636363636363635 Y2) (-4. + 2.5 Y0 + 1. Y1 +
0.5 Y2))) ==
0, -((0.007561436672967865 (1. + Y0^0.5) (-3. + Y0 +
2. Y1))/((3. + Y0^0.5 +
2. Y1^0.5) (2.0869565217391304 - 1. Y0 -
0.9565217391304348 Y1 - 0.13043478260869565 Y2)^2)) - (
0.01512287334593573 (1. + Y1^0.5) (-3. + Y0 + 2. Y1))/((3. +
Y0^0.5 + 2. Y1^0.5) (2.0869565217391304 - 1. Y0 -
0.9565217391304348 Y1 - 0.13043478260869565 Y2)^2) + (
0.031558185404339245 -
0.031558185404339245 Y0)/(1.2307692307692308 - 1. Y0 -
0.15384615384615385 Y1 -
0.07692307692307693 Y2)^2 + (-0.3282051282051282 +
0.2512820512820513 Y0 + 0.05128205128205128 Y1 +
0.02564102564102564 Y2)/((-1.2307692307692308 + 1. Y0 +
0.15384615384615385 Y1 + 0.07692307692307693 Y2) (-1.6 +
1. Y0 + 0.4 Y1 + 0.2 Y2)) + (
0.2512820512820513 (-1. + Y0) (-1.306122448979592 + 1. Y0 +
0.20408163265306123 Y1 + 0.10204081632653061 Y2))/((1. +
Y0^0.5) Y0^0.5 (-1.2307692307692308 + 1. Y0 +
0.15384615384615385 Y1 + 0.07692307692307693 Y2) (-1.6 +
1. Y0 + 0.4 Y1 + 0.2 Y2)) + (
0.07173913043478261 (1. + Y0^0.5) (-1.9393939393939394 +
1. Y0 + 0.787878787878788 Y1 +
0.15151515151515155 Y2))/((3. + Y0^0.5 +
2. Y1^0.5) (-2.0869565217391304 + 1. Y0 +
0.9565217391304348 Y1 + 0.13043478260869565 Y2) (-1.6 +
1. Y0 + 0.4 Y1 + 0.2 Y2)) + (
0.14347826086956522 (1. + Y1^0.5) (-1.9393939393939394 +
1. Y0 + 0.787878787878788 Y1 +
0.15151515151515155 Y2))/((3. + Y0^0.5 +
2. Y1^0.5) (-2.0869565217391304 + 1. Y0 +
0.9565217391304348 Y1 + 0.13043478260869565 Y2) (-1.6 +
1. Y0 + 0.4 Y1 + 0.2 Y2)) - (
0.035869565217391305 (1. + Y0^0.5) (-3. + Y0 +
2. Y1) (-1.9393939393939394 + 1. Y0 +
0.787878787878788 Y1 + 0.15151515151515155 Y2))/(
Y0^0.5 (3. + Y0^0.5 + 2. Y1^0.5)^2 (-2.0869565217391304 +
1. Y0 + 0.9565217391304348 Y1 +
0.13043478260869565 Y2) (-1.6 + 1. Y0 + 0.4 Y1 +
0.2 Y2)) - (
0.07173913043478261 (1. + Y1^0.5) (-3. + Y0 +
2. Y1) (-1.9393939393939394 + 1. Y0 +
0.787878787878788 Y1 + 0.15151515151515155 Y2))/(
Y0^0.5 (3. + Y0^0.5 + 2. Y1^0.5)^2 (-2.0869565217391304 +
1. Y0 + 0.9565217391304348 Y1 +
0.13043478260869565 Y2) (-1.6 + 1. Y0 + 0.4 Y1 +
0.2 Y2)) + (
0.10760869565217392 (-3. + Y0 + 2. Y1) (-1.9393939393939394 +
1. Y0 + 0.787878787878788 Y1 + 0.15151515151515155 Y2))/(
Y0^0.5 (3. + Y0^0.5 + 2. Y1^0.5) (-2.0869565217391304 +
1. Y0 + 0.9565217391304348 Y1 +
0.13043478260869565 Y2) (-1.6 + 1. Y0 + 0.4 Y1 +
0.2 Y2)) == 0, 0 <= Y2 <= 1, 0 <= Y1 <= 1, 0 <= Y0 <= 1}, {Y2,
Y1, Y0}]

• Drop 0 <= Y2 <= 1, 0 <= Y1 <= 1, 0 <= Y0 <= 1}, use FindRoot instead of Solve, and replace {Y2, Y1, Y0} with {{Y2, 0.3}, {Y1, 0.2}, {Y0, 0.1}}. You should get as the answer {Y2 -> 0.169988, Y1 -> 0.209275, Y0 -> 0.290306} (which satisfies your constraints).
– JimB
Nov 15, 2018 at 16:55
• The search starting values for FindRoot are not critical in this instance. Look at SeedRandom[0]; Table[ FindRoot[eqns, {{Y2, RandomReal[]}, {Y1, RandomReal[]}, {Y0, RandomReal[]}}], {10}] Nov 15, 2018 at 19:39

Call your expression e0. I removed the inequalities. Then we can make this into a system of bona fide polynomial equations in a few steps.

First rationalize to make the exponents into rationals. This also happenstance makes all coefficients into rationals.

e1 = MapAll[Expand, Rationalize[e0]];


There are a few variables in radicals. We replace them by squares (works fine since we only want positive solutions anyway).

e2 = Together[PowerExpand[e1 /. {Y0 -> y0^2, Y1 -> y1^2, Y2 -> y2^2}]];


We have denominators. We'll take a look at their square-less factors.

dene2 = Denominator[e2];
dens = Rest[FactorList[Apply[PolynomialLCM, dene2]]][[All, 1]]

(* Out[222]= {y0, y1,
3 + y0 + 2 y1, y2, -8 + 5 y0^2 + 2 y1^2 + y2^2, -16 + 13 y0^2 +
2 y1^2 + y2^2, -48 + 23 y0^2 + 22 y1^2 + 3 y2^2, -16 + 5 y0^2 +
2 y1^2 + 9 y2^2} *)


These are not allowed to vanish. Four our purposes (I checked this) it suffices just to make sure the individual variables do not vanish. We enforce that by adding a polynomial y0*y1*y2*yrecip - 1 (with a new variable yrecip that we really don't care about).

We now have a system that NSolve can handle in reasonable time.

AbsoluteTiming[
sols = NSolve[Join[Numerator[e2], {y0*y1*y2*yrecip - 1}]];]

(* Out[228]= {5.081856, Null} *)


We gather the solutions of interest. First remove complex values, then remove those that fall outside the range of interest, and last take squares to account for the earlier change of variables.

solvals = {y0, y1, y2} /. sols;
realvals = Select[solvals, FreeQ[#, Complex] &];
Select[realvals, AllTrue[#, 0 < # < 1 &] &]^2

(* Out[248]= {{0.290306173455, 0.209275354184, 0.169987697235}} *)


You can use NMinimize instaed of FindRoot[]:

If you set the equations inside of FindRoot (" first argument") by eq the solution is easily evaluated by

NMinimize[{1, eq}, {Y0, Y1, Y2}]
(*{1., {Y0 -> 0.290306, Y1 -> 0.209275, Y2 -> 0.169988}}*)


There is no need to define starting values ...