"solve " function returns inaccurate solutions

10 次查看(过去 30 天)
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.

采纳的回答

John D'Errico
John D'Errico 2019-1-21
编辑:John D'Errico 2019-1-21
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.

更多回答(4 个)

Birdman
Birdman 2019-1-21
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.
  6 个评论
Walter Roberson
Walter Roberson 2019-1-21
For your release, leave out the ==0
syms x real
solx=vpa(solve(x-sqrt(x+1)-sqrt(x-1),x))

请先登录,再进行评论。


madhan ravi
madhan ravi 2019-1-21
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])
  3 个评论
madhan ravi
madhan ravi 2019-1-21
编辑:madhan ravi 2019-1-21
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
moh pouladin 2019-1-21
yes my software is R2011b.
I wanted to know the cause behind the wrong answer, which produced by my machine.

请先登录,再进行评论。


Walter Roberson
Walter Roberson 2019-1-21
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.
  6 个评论
Walter Roberson
Walter Roberson 2019-1-21
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
John D'Errico 2019-1-21
编辑:John D'Errico 2019-1-21
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.

类别

Help CenterFile Exchange 中查找有关 Assumptions 的更多信息

标签

产品


版本

R2011b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by