1
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – Rico
    Feb 7, 2017 at 23:26
  • $\begingroup$ @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! $\endgroup$
    – larry
    Feb 8, 2017 at 0:14
  • 1
    $\begingroup$ According to the documentation, NSolve[expr,vars,Reals] finds solutions over the domain of real numbers. $\endgroup$
    – Felix
    Feb 8, 2017 at 0:48
  • $\begingroup$ @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. $\endgroup$
    – larry
    Feb 8, 2017 at 1:06
  • 2
    $\begingroup$ 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"] $\endgroup$
    – Daniel Lichtblau
    Feb 8, 2017 at 16:10

1 Answer 1

Reset to default
6
$\begingroup$ 
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

enter image description here

EDIT: Using NSolve

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

enter image description here

EDIT 2: Or use high precision calculations

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

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.