solve() gives incorrect answer

I am having trouble with the solve() command in Matlab R2012a 64 bit. I will supply an example to show what I mean. I stated the specificiation of my computer as well, I do not know if that is considered useful.

----------------------------

COMPUTER SPECIFICATIONS

----------------------------

Windows 7 Enterprise SP 1

64-bit OS

Intel Core i7-2670QM CPU 2.2 GHz

8,00 GB RAM

----------------------------

EXAMPLE

----------------------------

I have a function

G = (- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2)).

I am interested in the intersection with y = -1.

Both

solve(G+1)

and

solve('(- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2)) == -1')

give me

ans = 
1

This is only one of the intersections, as can be seen from the figure below. Does anyone have an idea what I might be overlooking?

Thanks in advance,

Jori

[Walter Roberson]

The other solution is the rather messy

69/400 + (1/1200)*3^(1/2)*((27883*(127918+30*315201^(1/2))^(1/3) + 400*(127918+30*315201^(1/2))^(2/3) + 1009600)/(127918 + 30*315201^(1/2))^(1/3))^(1/2) + (1/1200)*6^(1/2)*((27883*(127918 + 30*315201^(1/2))^(1/3) - 200*(127918 + 30*315201^(1/2))^(2/3) - 504800 + 1193127*3^(1/2)*((127918 + 30*315201^(1/2))^(1/3)/(27883*(127918 + 30*315201^(1/2))^(1/3) + 400*(127918 + 30*315201^(1/2))^(2/3) + 1009600))^(1/2)*(127918 + 30*315201^(1/2))^(1/3))/(127918 + 30*315201^(1/2))^(1/3))^(1/2)

... for whatever that is worth ;-)

Answer 1

[Jonathan Epperl]

I am astounded... So I am using 2011b on a Mac, but that shouldn't account for the differences. I am getting different results from the ones you are getting.

syms x
G = (- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2))
solve(G+1)

yields

ans =

                                                                            1
                                           0.80485291557706605237097971857745
                                           0.29354887914157893712269217217225
 - 0.20420089735932249474683594537485 + 0.025051677816927281916779859288826*i
 - 0.20420089735932249474683594537485 - 0.025051677816927281916779859288826*i

Which is very funny, because the last 3 results correspond to G==1, so they're wrong.

If I change the last line

solve(G+1,'MaxDegree',4)

Then I get

    ans =

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        1
   %snipped a very long expression

vpa(ans)

    ans =

                                  1.0
   0.80485291557706605237097971857745

 So now those are the correct expressions.

Then what you posted yields

solve('(- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2)) == -1')
Error using solve>getEqns (line 292)
' (- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2)) == -1 ' is not a valid expression or equation.

Error in solve (line 141)
[eqns,vars,options] = getEqns(varargin{:});

And this here ('==' changed to '=')

 solve('(- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2)) = -1')

again yields the 5 solutions from above, including the three wrong ones.

Did you maybe assign some value to x at any point, or did anything else that could lead to such weird behavior? If you issue clear all and type my exact commands, do you get different results?

In any case, two solutions:

1. Add 'MaxDegree',4 and hope for the best
2. Use Mathematica or even just Wolfram|Alpha, which solved your example perfectly when I tried it: http://wolfr.am/UvC8yp

[Walter Roberson]

Recognition of the == for symbolic expressions was new in R2012a I believe (might have been R2012b, but after R2011b for sure.)

[Jori]

Adding 'MaxDegree',4 works.

>> clear all
>> solve('(- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2))+1','MaxDegree',4)

ans =

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        1
and a very long symbolic expression

>> vpa(ans)

ans =

                                  1.0
   0.80485291557706605237097971857745

Thank you for your help. The documentation of the solve() function states for 'MaxDegree'

solve(...,'MaxDegree',n) controls the maximum degree of polynomials
    for which explicit formulas will be used during the computation.
    n must be a positive integer smaller than 5. The default is 3.

If I understand it correctly, a higher order polynomial allows for computation of 'harder' intersections? Do you have a better interpretation of this option?

Thanks! Jori

[Jonathan Epperl]

I think your interpretation is pretty accurate, even though I don't know what Matlab does if it encounters a polynomial with degree greater than MaxDegree. Your problem indicates that it goes ahead and seeks one numerical solution and continues with that. I think Mathematica leaves stuff like Root[x^6+3 x^4-x+1] in the output. Anyway, glad it works!

[Jori]

It seems to be a recurring problem, say I have the function

>> clear all
>> syms x
>> G = (x - (23*x^2)/2 + 1/2)/(10*(x^5)^(1/2))

And I compute

>> m = 0; s = 4; double(solve(G == (2*m-s)/s,'MaxDegree',4))

ans =

      1.0000
      0.4388

Which is correct. But if I try

>> m = 1; s = 4; double(solve(G == (2*m-s)/s,'MaxDegree',4))

ans =

      5.0935

While it should give me another intersection at about x = 0.3.

Would it be a good idea to replace my installation of Matlab with an earlier version? Should I report this recurring problem to Mathworks? Anyone with any ideas on what is going on?

Should I make a new post, such that others might see it as well?

[Jori]

And when I try

>> m = 3; s = 4; solve(G == (2*m-s)/s,'MaxDegree',4)
Error using mupadmex
Internal error with symbolic engine. Please quit and
restart MATLAB.

Error in mupadengine/evalin (line 97)
                [res,status] = mupadmex(statement);

Error in mupadengine/feval (line 150)
            [S,err] = evalin(engine,stmt);

Error in solve (line 160)
  sol = eng.feval('symobj::solvefull',eqns,vars);|

[Jonathan Epperl]

Hm, that is mysterious... Here is what my Matlab does if I enter your commands:

>> syms x real
>> G = (x - (23*x^2)/2 + 1/2)/(10*(x^5)^(1/2))

G =

(x - (23*x^2)/2 + 1/2)/(10*(x^5)^(1/2))

>> m = 1; s = 4; (solve(G - ((2*m-s)/s),'MaxDegree',4))

ans =

                                             5.0935065938427671851921446479193
                                            0.30296905750943425877783588487726
                                            0.22908035767961411609832320482884
   - 0.1677780045159077800341518688127 + 0.011755832620508898266757829430043*i
   - 0.1677780045159077800341518688127 - 0.011755832620508898266757829430043*i

>> m = 0; s = 4; (solve(G - ((2*m-s)/s),'MaxDegree',4))

ans =

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 1
   % Snip

>> vpa(ans)

ans =

                                  1.0
   0.43879264431511722778976907661477

>> m = 3; s = 4; (solve(G - ((2*m-s)/s),'MaxDegree',4))

ans =

                                             5.0935065938427671851921446479193
                                            0.30296905750943425877783588487726
                                            0.22908035767961411609832320482884
   - 0.1677780045159077800341518688127 + 0.011755832620508898266757829430043*i
   - 0.1677780045159077800341518688127 - 0.011755832620508898266757829430043*i

So I am not seeing the error you are seeing, and instead of yielding too few answers, solve returns wrong ones. And complex ones even after I declared x to be real. (You'll notice I changed your syntax from G == ... to G - ..., my version of the symbolic toolbox doesn't like the '==' and it seems better practice doing that anyway).

Reporting that to The MathWorks seems a good idea. I have 2012b at home, I'll try tonight and see whether that yields yet other results... Apart from hoping the MathWorks to fix anything (or us finding an error in what we are doing), if you have access to Mathematica, or Maple or any other symbolic math software you could try that. Like I said, Wolfram|Alpha (which is free) solved your first problem accurately.

Answer 2

[Matt Fig]

Try this way:

>> solve('(- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2))+1')

ans =

                                      1
                         0.804852915577
                         0.293548879142
   - 0.204200897359 + 0.0250516778169*i
   - 0.204200897359 - 0.0250516778169*i

[Jori]

I tried it, but it does not work. This seems pretty weird.

>> clear all
>> solve('(- 13*x^2 + 2*x + 1)/(10*(x^5)^(1/2))+1')

ans =

1

Could it have to do with the version of Matlab I am using?