Asked by moh pouladin
on 21 Jan 2019

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 by John D'Errico
on 21 Jan 2019

Edited by John D'Errico
on 21 Jan 2019

Accepted Answer

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.

Answer by Birdman
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.

moh pouladin
on 21 Jan 2019

Walter Roberson
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))

moh pouladin
on 21 Jan 2019

produced the wrong answer, too.

Sign in to comment.

Answer by madhan ravi
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])

moh pouladin
on 21 Jan 2019

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
on 21 Jan 2019

moh pouladin
on 21 Jan 2019

yes my software is R2011b.

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

Sign in to comment.

Answer by Walter Roberson
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.

moh pouladin
on 21 Jan 2019

moh pouladin
on 21 Jan 2019

R2011b

Walter Roberson
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.

Sign in to comment.

Answer by John D'Errico
on 21 Jan 2019

Edited by John D'Errico
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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.