Why does NSolve give different results in different version of Mathematica?

Question

I am working on a problem with a colleague who uses Mathematica v.10.2 while I have v.11.0 Student Edition. We are solving the exact same equations but get two slightly different results. Why would this be the case? The equations are below

eq1 = 1/8 (3/b + 1/(4 + b) + (4 + a)/(4 + 2 a + (4 + a) b) + (2 (2 + a))/(
     2 b + a (2 + b)) + (4 + a)/(4 (5 + b) + a (6 + b))) == 1

eq2 = 1/(2 + a) + 
   1/2 (4/(8 a + 6 a^2 + a^3) + 1/((4 + a) (4 + 2 a + (4 + a) b)) + 
      2/((2 + a) (2 b + a (2 + b))) + 
      1/((4 + a) (4 (5 + b) + a (6 + b)))) == 1

We are looking for a particular solution (real, positive, less than 1), but for arguments sake we have solved it with NDSolve without restricting the solution. The last solution in the list is the one we want.

NSolve[{eq1,eq2},{a,b}]

Solution with v.11.0 Student Edition:

{a -> 0.454152, b -> 0.603878}

Solution with v.10.2:

{a -> 0.454231 + 0.0000159456 I, b -> 0.604194 + 4.76104*10^{-7} I}

Why do we not agree on this result?

Answer 1

Use Solve and then numericize the result (Root objects)

eq1 = 1/8 (3/b + 
      1/(4 + b) + (4 + a)/(4 + 2 a + (4 + a) b) + (2 (2 + a))/(2 b + 
         a (2 + b)) + (4 + a)/(4 (5 + b) + a (6 + b))) == 1;

eq2 = 1/(2 + a) + 
    1/2 (4/(8 a + 6 a^2 + a^3) + 1/((4 + a) (4 + 2 a + (4 + a) b)) + 
       2/((2 + a) (2 b + a (2 + b))) + 1/((4 + a) (4 (5 + b) + a (6 + b)))) ==
    1;

sol = Solve[{eq1, eq2}, {a, b}] // RootReduce // N[#, 10] &

(* {{a -> -5.847070681, b -> -8.057765676}, {a -> -3.297542697, 
  b -> 0.4665224792}, {a -> -2.660974148, 
  b -> 1.193835974}, {a -> -2.422609025, 
  b -> -3.329443943}, {a -> -1.761191102, 
  b -> 0.4792402029}, {a -> -1.668415568, 
  b -> -3.854562268}, {a -> -1.218484056, 
  b -> -0.4911704356}, {a -> -1.200634769, 
  b -> -4.482336957}, {a -> -0.5235331275, 
  b -> 1.150513108}, {a -> 0.3265558204, 
  b -> -3.886887946}, {a -> 0.3919689993, 
  b -> -1.014016456}, {a -> 0.4541518692, 
  b -> 0.6038775534}, {a -> 0.4743429352, 
  b -> -5.008727663}, {a -> 1.190322971, 
  b -> -0.5685459315}, {a -> -5.941745405 - 0.347087341 I, 
  b -> -3.905214957 + 0.324730348 I}, {a -> -5.941745405 + 0.347087341 I, 
  b -> -3.905214957 - 0.324730348 I}, {a -> -3.784060936 - 0.331813229 I, 
  b -> -3.954766281 + 0.331343393 I}, {a -> -3.784060936 + 0.331813229 I, 
  b -> -3.954766281 - 0.331343393 I}, {a -> -2.934304037 - 0.967402950 I, 
  b -> -3.740284783 + 1.954661627 I}, {a -> -2.934304037 + 0.967402950 I, 
  b -> -3.740284783 - 1.954661627 I}} *)

Verifying the solutions

And @@ (eq1 && eq2 /. sol)

(* True *)

EDIT: For real, positive solutions

solR = Solve[{eq1, eq2, a > 0, b > 0}, {a, b}] // RootReduce // N[#, 10] &

(* {{a -> 0.4541518692, b -> 0.6038775534}} *)