"solve " function returns inaccurate solutions

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.

Answer 1

[John D'Errico]

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.

[moh pouladin]

Thanks alot.

Answer 2

[Birdman]

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.

[Torsten]

This doesn't answer the question why x=1.11508... is returned as a solution.

[moh pouladin]

I tried the following command:

assume(x,'real')

and I recieved the following error:

"Undefined function 'assume' for input arguments of type 'sym'."

[Walter Roberson]

You are either using R2011b or earlier, or you are using the Maple based symbolic toolbox.

You might be able to use

syms x real

[moh pouladin]

[Walter Roberson]

For your release, leave out the ==0

syms x real
solx=vpa(solve(x-sqrt(x+1)-sqrt(x-1),x))

[moh pouladin]

produced the wrong answer, too.

Answer 3

[madhan ravi]

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])

[moh pouladin]

your first syntax produced the only correct answer.

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

[madhan ravi]

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.

[moh pouladin]

yes my software is R2011b.

I wanted to know the cause behind the wrong answer, which produced by my machine.

Answer 4

[Walter Roberson]

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.

[moh pouladin]

your syntax produced the following error on my machine:

"Error using solve>getEqns (line 292)

' 4 ' is not a valid expression or equation.

Error in solve (line 141)

[eqns,vars,options] = getEqns(varargin{:});"

thanks

[madhan ravi]

post the screenshot of the results that you got

[John D'Errico]

What MATLAB release?

[moh pouladin]

[moh pouladin]

R2011b

[Walter Roberson]

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.

Answer 5

[John D'Errico]

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.