"solve " function returns inaccurate solutions
8 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
moh pouladin
el 21 de En. de 2019
Comentada: moh pouladin
el 22 de En. de 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.
0 comentarios
Respuesta aceptada
John D'Errico
el 21 de En. de 2019
Editada: John D'Errico
el 21 de En. de 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.
Más respuestas (4)
Birdman
el 21 de En. de 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.
6 comentarios
Walter Roberson
el 21 de En. de 2019
For your release, leave out the ==0
syms x real
solx=vpa(solve(x-sqrt(x+1)-sqrt(x-1),x))
madhan ravi
el 21 de En. de 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])
3 comentarios
madhan ravi
el 21 de En. de 2019
Editada: madhan ravi
el 21 de En. de 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.
Walter Roberson
el 21 de En. de 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.
6 comentarios
Walter Roberson
el 21 de En. de 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
el 21 de En. de 2019
Editada: John D'Errico
el 21 de En. de 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.
0 comentarios
Ver también
Categorías
Más información sobre Assumptions en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!