MATLAB: Solve returns different answers to the same question

solvesymbolicsymbolic toolboxsystem of equations

I am trying to use solve() to solve a system of 4 equations for 4 variables, using the following code.
syms x1 x2 x3 x4
vars = [x1 x2 x3 x4];
eqns = [ x1 == 30/7 - 1/(x3+x4), x2==30/7 - 1/(x3+x4), x3== 10/3 - 1/(x1+x2+x4), x4==10/3 - 1/(x1+x2+x3) ];
C = solve(eqns, vars);
Then vpa(C.x1) returns:
ans =
                                           -1.6666666666666666666666666666667
                                           -1.6666666666666666666666666666667
   4.1316999594422134359804216941161 - 2.5330953622691428389795980612384e-39i
    4.343827578149360853360217735709 + 9.3940138280508998241066580657853e-40i
 0.095901033836997139230789141603439 + 1.6227020204672685655672638821872e-39i
Which is the correct solution. However, when the same problem is reformulated using:
eqns = [ x1 == 10/3 - 1/(x2+x3+x4), x2==10/3 - 1/(x1+x3+x4), x3== 30/7 - 1/(x1+x2), x4==30/7 - 1/(x1+x2) ];
Note that although the names of the variables and the order they are entered is different, this system is identical to the original. Then vpa(C.x4) (which should contain the same solutions as C.x1) instead returns:
ans =
  -1.6666666666666666666666666666667
  -1.6666666666666666666666666666667
    4.343827578149360853360217735709
 0.095901033836997139230789141603439
   4.1316999594422134359804216941161
This is presenting a problem for me since somehow it seems to be losing the imaginary portion of each of the complex solutions, only in the second scenario. Why is this happening and how can I fix it so that it is returning the full complex solutions?

Best Answer

It is possible to demonstrate that one of the roots you are seeing imaginary portions for has a closed form value of
 10/63 * sqrt(3) * sqrt(499) * cos(atan(63/1596949 * sqrt(3) * sqrt(203226157))/3 + pi/6) - 110/63 - 10/63 * sqrt(499) * sin(atan(63/1596949 * sqrt(3) * sqrt(203226157))/3 + pi/6)
That involves two arctan of the same expression 63*sqrt(3)*sqrt(203226157)*(1/1596949) . You can evaluate that expression and see that it is a number between 0 and 1; therefore the arctan is real-valued. Therefore the expression as a whole is real valued.
The imaginary components you see do not exist and are due to round-off error. If you increase the number of digits in the vpa then the imaginary components will get smaller and smaller.
... actually I need to revise that closed form solution, as I was not paying attention to the way the variables get renamed. I will come back to this with a correction; I need to do something else for now.
Related Question