Why does the SOLVE function return the wrong answers for some equations?

Why does the SOLVE function return the wrong answers for some equations?

SOLVE returns incorrect symbolic solutions to some symbolic equations. For example, attempting to solve the equation z^6 = i results in incorrect answers:

solve('z^6-i')

 ans =

 [  (1/4*5^(1/2)-1/4+1/4*i*2^(1/2)*(5+5^(1/2))^(1/2))*(-i)^(1/2)]

 [ (-1/4*5^(1/2)-1/4+1/4*i*2^(1/2)*(5-5^(1/2))^(1/2))*(-i)^(1/2)]

 [ (-1/4*5^(1/2)-1/4-1/4*i*2^(1/2)*(5-5^(1/2))^(1/2))*(-i)^(1/2)]

 [  (1/4*5^(1/2)-1/4-1/4*i*2^(1/2)*(5+5^(1/2))^(1/2))*(-i)^(1/2)]

 [                                     1/2*2^(1/2)-1/2*i*2^(1/2)]

 [                                    -1/2*2^(1/2)+1/2*i*2^(1/2)]

Checking the accuracy of these values by performing

double(ans.^6)

shows that the first 4 answers are invalid:

 ans =

   -0.9511 + 0.3090i

   -0.5878 - 0.8090i

    0.5878 - 0.8090i

    0.9511 + 0.3090i

         0 + 1.0000i

         0 + 1.0000i

Answer 1

[MathWorks Support Team]

This issue has been fixed in Symbolic Math Toolbox 3.0.1 (R13+). If you are using an earlier versions:

This bug may be occurring in certain types of polynomials with complex coefficients, so you should watch for further such errors using these polynomials. As a workaround, you can use SOLVE with a similar problem involving a real equation whose roots include those of the polynomial with complex coefficients you are attempting to solve.

If you were working with the polynomial z^6-i, to find the appropriate polynomial you would multiply by the conjugate, z^6+i, to get z^12+1. The solutions to this equation will include those of z^6-i.

The following code finds the six solutions for this problem:

x = solve(z^12+1);  % x holds 12 solution strings

y = double(x.^6);   % contents of y are now y = i or y = -i

ix = (imag(y) > 0); % find indices where y = i

soln = x(ix);       % soln holds the six solution strings where y = i

If you are looking for approximate answers to a polynomial, you can also use SYM2POLY and ROOTS:

syms z

f = z^6-i;  % your polynomial here

p = sym2poly(f);

roots(p)

 ans =

   -0.9659 - 0.2588i

   -0.7071 + 0.7071i

   -0.2588 - 0.9659i

    0.2588 + 0.9659i

    0.7071 - 0.7071i

    0.9659 + 0.2588i