# Nsolve gives solutions with a small complex part which can't be eliminated by chop command

I do see that what I'm asking has been asked many times, and I have checked the previous questions, but they didn't work for the specific problem I'm asking: I have six equations with six unknowns:

Sqrt[2] x[2] y[1] + x[1] (x[2] - Sqrt[2] y[2]) == 1

Sqrt[2] x[3] y[2] + x[2] (x[3] - Sqrt[2] y[3]) == 1

x[1] (x[3] + Sqrt[2] y[3]) == 1 + Sqrt[2] x[3] y[1]

x[1]^2 + y[1]^2 == 1

x[2]^2 + y[2]^2 == 1

x[3]^2 + y[3]^2 == 1


and I have tried the following code to calculate the solutions:

sol = NSolve[{eq21, eq22, eq23, eq24, eq25, eq26},
WorkingPrecision -> MachinePrecision] // N // Chop


the problem is, all the solutions includes a complex part of the magnitude of around $10^{-8}$, which can't be chopped off by mathematica. One of the solutions looks like:

{x[1.] -> 1., x[2.] -> 1., x[3.] -> 1.,   y[1.] -> 2.84087*10^-7 - 5.67727*10^-8 I,   y[2.] -> 2.84087*10^-7 - 5.67727*10^-8 I,  y[3.] -> 2.84087*10^-7 - 5.67727*10^-8 I}}


I know that the above system has at least one solution of $[1,1,1,0,0,0]$, so this does not make sense.

For this problem, I could use the Solve command to get the exact solutions, but I will need to take care of more complex coefficients, which is shown to be difficult to be handled by Solve.

So I wonder how to get rid of the complex part, and get the precision of the result to be good? For example, I got some solution that gives me x[i]=1, but y[i] as a small complex number, while it should be 0.

Thanks a lot!

Update I

@Bob, thanks a lot for your code, but somehow, mathematica gives me some different result, especially when I use Nsolve, which returns an empty set. Any idea of what's going on? and could you please give a little bit more explanation about your code to help me understand it?

@ Daniel, even though I do not understand what it means, this method works perfect for this problem, even if I have some complex coefficients. Due to the complexity of the set of equations, I probably need to use Nsolve rather than Solve, and this method fixs the problem perfectly.

• Can't you just tune the tolerance? Like Chop[sol,10^-2] so all your real parts survive while the complex ones which are of the order 10^-7 will get chopped. – Rico Feb 7 '17 at 23:26
• @Rico, thanks for your suggestion, this does work, but I hope to find a more robust way to solve the equations, because the error is not small compared to the machine precision, which makes me not quite confident about how much error it will have when I change the coefficients to more complex form. But your method does work for simple cases, which I will use for now. Thanks a lot! – larry Feb 8 '17 at 0:14
• According to the documentation, NSolve[expr,vars,Reals] finds solutions over the domain of real numbers. – Felix Feb 8 '17 at 0:48
• @Felix, yes, but when I try this approach, it returns an empty set to me, indicating that it thinks that there is no real solution, so this approach does not work here. – larry Feb 8 '17 at 1:06
• It has to do with some amount of inaccuracy in the default method as it wends its way through C^n in path-tracking by homotopy. Sometimes the older method is better at discerning real solutions. So one might try NSolve[...,Method -> "EndomorphismMatrix"] – Daniel Lichtblau Feb 8 '17 at 16:10

eqns = {
Sqrt[2] x[2] y[1] + x[1] (x[2] - Sqrt[2] y[2]) == 1,
Sqrt[2] x[3] y[2] + x[2] (x[3] - Sqrt[2] y[3]) == 1,
x[1] (x[3] + Sqrt[2] y[3]) == 1 + Sqrt[2] x[3] y[1],
x[1]^2 + y[1]^2 == 1,
x[2]^2 + y[2]^2 == 1,
x[3]^2 + y[3]^2 == 1};

vars = Cases[eqns, x[_] | y[_], Infinity] // Union

(*  {x[1], x[2], x[3], y[1], y[2], y[3]}  *)

solns = Solve[eqns, vars];


Verifying solutions

And @@@ (eqns /. solns)

(*  {True, True, True, True, True, True, True, True}  *)

vars /. solns // Column


EDIT: Using NSolve

(solns2 = ((vars /. NSolve[eqns, vars]) /.
z_?NumericQ :> Round[Chop[z, 10^-6], 10^-6] //
Union) /. z_Rational :> N[z]) // Column


EDIT 2: Or use high precision calculations

(soln3 = vars /. NSolve[eqns, vars, Reals, WorkingPrecision -> 30] // Union //
Round[#, 10.^-6] &) // Column


• Very good answer. Any idea why NSolve[eqns, vars, Reals] does not work? – Felix Feb 8 '17 at 4:25
• @Felix - use high precision calculations. See edit 2 – Bob Hanlon Feb 8 '17 at 4:40