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

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?

• I've been using Mathematica so long, it's hard for me to put myself in your frame of mind. Why would they return the same result? Different versions are likely to try using different numerical methods. And the newer versions tend to chose better numerical methods for solving the problem. Oct 29, 2018 at 15:39
• Just to be clear, NSolve is not a well defined operation or algorithm. NSolve tries to find the best way to solve the problem that it can. There's no way for it to know beforehand what "the best way" to solve the problem is, so it has to use heuristics and these will change between versions as it gets better. Oct 29, 2018 at 15:42
• I can't explain the different result, but in v10.1 FindRoot[{eq1, eq2}, {{a, 1}, {b, 1}}] yields {a -> 0.454152, b -> 0.603878}, and does so two orders of magnitude faster than NSolve. Oct 29, 2018 at 15:54
• I obtain identical results in versions 7.0 and 11.3. Oct 29, 2018 at 19:38
• If I am not mistaken, those versions used different default methods. Oct 30, 2018 at 16:21

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}} *)