## "solve " function returns inaccurate solutions

on 21 Jan 2019
Latest activity Commented on by moh pouladin

on 22 Jan 2019

### John D'Errico (view profile)

Hi.
I tried to solve the following equation with "solve" command:
x-[\sqrt{x+1}+\sqrt{x-1}] = 0
I entered the following command:
solve('x-sqrt(x+1)-sqrt(x-1)','x')
and recieved 4 solution, one of them was, "1.11508....", which is not correct. there was also the correct solution between the answers , which is 3.9343...
what's the problem, how can I trust the answers produced by commands such as 'solve' command ?
Thanks alot.

R2011b

### John D'Errico (view profile)

on 21 Jan 2019
Edited by John D'Errico

### John D'Errico (view profile)

on 21 Jan 2019

Ok. It looks like your real question is where does the spurious root arise? And, since you are using an old enough release of MATLAB that it will only run on a computer old enough that it needs a crank on the side to start it, I can't run that release to compare. But we can see from whence the spurious solution arises.
First, we need to recognize what happens with a square root. The sqrt function actually has TWO branches, even though people only tend to think of one. So if u=sqrt(x), then so is -u a valid square root. Writing your relation as...
x - sqrt(x+1) - sqrt(x-1) == 0
we really need to think of it as 4 possible problems, thus these three others.
x + sqrt(x+1) + sqrt(x-1) == 0
x + sqrt(x+1) - sqrt(x-1) == 0
x - sqrt(x+1) + sqrt(x-1) == 0
So lets plot all 4 relations on one set of axes. In the plot, sqrt will always takes the positive branch only.
ezplot(@(x) x - sqrt(x+1) - sqrt(x-1),[-1,5])
hold on
ezplot(@(x) x + sqrt(x+1) + sqrt(x-1),[-1,5])
ezplot(@(x) x + sqrt(x+1) - sqrt(x-1),[-1,5])
ezplot(@(x) x - sqrt(x+1) + sqrt(x-1),[-1,5])
refline([0,0])
As you should see, there are 4 different curves there. distinguished by the 4 colors plotted. As well, there are TWO solutions, since both the cyan and blue curves cross the reference line at y==0.
So truly, there are two distinct and fully valid real solutions, as long as you agree that the square root function has two branches.
Your version of solve found both of them, one of which is apparently at 1.11508. It arises as the solution to this variation of your problem (I am using the R2018b version):
syms x
vpa(solve(x-sqrt(x+1)+sqrt(x-1)))
ans =
1.1150879946798484383326671302376
So my version of solve seems to be looking at the positive branch of the sqrt function only.
It is simple enough to solve the problem in a way that generates the spurious solution too.
x = sqrt(x-1) + sqrt(x+1)
Square both sides, to get
x^2 = (x-1) + (x+1) + 2*sqrt((x-1)*(x+1))
Collect, then square again.
(x^2 - 2*x)^2 = 4*(x-1)*(x+1)
See that by squaring things, we resolve those +/- ambiguities. But at the same time, we may introduce spurious solutions, since sqrt as a function tends to imply only the positive square root to many people who want to ignore the negative branch.
That 4th degree polynomial has 4 roots, some of which may be complex. They are rather messy to write in full form, as you might expect.
R = solve((x^2 - 2*x)^2 == 4*(x-1)*(x+1),'maxdegree',4)
R =
1 - (3^(1/2)*(24*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 16*3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 3*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 64*3^(1/2)*6^(1/2)*(3^(1/2)*23^(1/2) + 9)^(1/2))^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6)*(36*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 9*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 48)^(1/4)) - (3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6))
(3^(1/2)*(24*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 16*3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 3*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 64*3^(1/2)*6^(1/2)*(3^(1/2)*23^(1/2) + 9)^(1/2))^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6)*(36*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 9*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 48)^(1/4)) - (3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6)) + 1
(3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6)) - (3^(1/2)*(24*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 16*3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 3*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) + 64*3^(1/2)*6^(1/2)*(3^(1/2)*23^(1/2) + 9)^(1/2))^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6)*(36*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 9*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 48)^(1/4)) + 1
(3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6)) + (3^(1/2)*(24*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 16*3^(1/2)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) - 3*3^(1/2)*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3)*(12*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 3*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 16)^(1/2) + 64*3^(1/2)*6^(1/2)*(3^(1/2)*23^(1/2) + 9)^(1/2))^(1/2))/(6*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/6)*(36*((32*3^(1/2)*23^(1/2))/9 + 32)^(1/3) + 9*((32*3^(1/2)*23^(1/2))/9 + 32)^(2/3) + 48)^(1/4)) + 1
We can resolve them into decimal form using vpa.
vpa(R)
ans =
- 0.52470257992985177015834395726155 - 0.79777765915011117379630925220011i
- 0.52470257992985177015834395726155 + 0.79777765915011117379630925220011i
1.1150879946798484383326671302376
3.9343171651798551019840207842855
So there are two real solutions to the problem once converted into an associated 4th degree polynomial. But only one of them is a solution to the original problem as posed, IF you take only the positive branch of the sqrt.
In the end, you should see that solve did not return an inaccurate solution, but just a surprising one, because you did not take into account the two branches of a sqrt. Is it a spurious solution? I supose you can argue that either way.

on 22 Jan 2019
Thanks alot.

### Birdman (view profile)

on 21 Jan 2019

Try this:
syms x
assume(x,'real');
solx=vpasolve(x-sqrt(x+1)-sqrt(x-1)==0,x)
or
syms x
assume(x,'real');
solx=vpa(solve(x-sqrt(x+1)-sqrt(x-1)==0,x))
You can tell MATLAB that you want x variable to be assumed as a real number.

on 21 Jan 2019
Walter Roberson

### Walter Roberson (view profile)

on 21 Jan 2019
For your release, leave out the ==0
syms x real
solx=vpa(solve(x-sqrt(x+1)-sqrt(x-1),x))

on 21 Jan 2019

on 21 Jan 2019

fzero(@(x)x-sqrt(x+1)-sqrt(x-1),[1 5])
% ^^^----- domain
%or
syms x
sol=vpasolve(x-sqrt(x+1)-sqrt(x-1)==0,x,[1 5])

on 21 Jan 2019
your second suggested syntax produced the following error on my machine:
"Undefined function 'vpasolve' for input arguments of type 'sym'."
however I rather wanted to know, what's the problem with "solve" command, which produces wrong answers?
thanks

on 21 Jan 2019
You seem to be using 2011b version and vpasolve() was introduced in 2012b , plus in later release the correct is being observed . Try clear all and clear global at the very beginning and try your original code again.

on 21 Jan 2019
yes my software is R2011b.
I wanted to know the cause behind the wrong answer, which produced by my machine.

### Walter Roberson (view profile)

on 21 Jan 2019

syms x
solve(x-sqrt(x+1)-sqrt(x-1),'maxdegree', 4)
You will get the solution,
((9*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(2/3) - 3*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/3) - 20)^(1/2)/(6*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/6)) + (20*(9*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(2/3) - 3*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/3) - 20)^(1/2) - 6*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/3)*(9*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(2/3) - 3*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/3) - 20)^(1/2) - 9*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(2/3)*(9*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(2/3) - 3*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/3) - 20)^(1/2) + 6*2^(1/2)*6^(1/2)*(3*3^(1/2)*23^(1/2) + 11)^(1/2))^(1/2)/(6*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/6)*(9*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(2/3) - 3*((4*3^(1/2)*23^(1/2))/9 + 44/27)^(1/3) - 20)^(1/4)) + 1/2)^2 + 1
simplify(expand()) will give you a slightly more compact version of it.

on 21 Jan 2019

on 21 Jan 2019
R2011b
Walter Roberson

### Walter Roberson (view profile)

on 21 Jan 2019
Possibly you need
syms x
solve(x-sqrt(x+1)-sqrt(x-1), x, 'maxdegree', 4)
'maxdegree' is a valid option for your release, provided you are not using the Maple based symbolic engine.

### John D'Errico (view profile)

on 21 Jan 2019
Edited by John D'Errico

### John D'Errico (view profile)

on 21 Jan 2019

Works for me:
syms x
S = solve(x-sqrt(x+1)-sqrt(x-1))
S =
root(z^4 - 2*z^3 + 2*z^2 - 2*z - 1, z, 4)^2 + 1
vpa(S)
ans =
3.9343171651798551019840207842855
If you use 'maxDegree' as 4, as Walter suggests, you get an analytical expression, as he shows. But either way does give you the correct solution.
If you get something else, what MATLAB release are you using? My guess it it is an old release, since it accepts string input for the equation. One of the things we have been asking for is a flag when you post an answer, that tells readers which MATLAB release you are using. That would more easily help to resolve such issues.