Symbolic solver yields wrong solution
2 次查看(过去 30 天)
显示 更早的评论
Hello,
I am not sure if I am overlooking something very obvious or whether the MATLAB symbolic solver is giving me a wrong solution.
I have the following code:
syms p t
formula = 84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 + (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
solutions = solve(formula == 0, t)
vpasolutions = vpa(solutions)
Now MATLAB 2019a gives me the following result:
formula =
84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 + (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
solutions =
0
-p*2i
p*2i
-(p*2i)/3
(p*2i)/3
-(p*(272/27 - (64*13^(1/2))/27)^(1/2))/2
(p*(272/27 - (64*13^(1/2))/27)^(1/2))/2
-(p*((64*13^(1/2))/27 + 272/27)^(1/2))/2
(p*((64*13^(1/2))/27 + 272/27)^(1/2))/2
vpasolutions =
0
-p*2.0i
p*2.0i
-p*0.66666666666666666666666666666667i
p*0.66666666666666666666666666666667i
-0.61797697405792080417600009731734*p
0.61797697405792080417600009731734*p
-2.1575776918969228428670635119892*p
2.1575776918969228428670635119892*p
So altogether there seem to be 5 real and 4 complex solutions for any given 'p'. I am only interested in the real ones. Now I plot the graph for p = 1 using the following code:
t_values = linspace(-2.2,2.2,1000);
formula_values = subs(formula,{p,t},{1,t_values});
plot(t_values,formula_values)
grid on
I get the attached figure. Here one can see that my function should not be zero around -0.618 and +0.618, which seems to contradict the answers from the symbolic solver. Now I know that one cannot always trust a graph with discrete values but even if you decrease the step size, this does not seem to change, and wolfram alpha gives me all results except for those two ominous ones.
Furthermore, if I add the assumption that p must be positive, these ominous results also vanish:
syms p t
assume(p>0)
formula = 84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 + (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
solutions = solve(formula == 0, t)
vpasolutions = vpa(solutions)
yields
formula =
84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 + (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
solutions =
0
-p*2i
p*2i
-(p*2i)/3
(p*2i)/3
-(p*((64*13^(1/2))/27 + 272/27)^(1/2))/2
(p*((64*13^(1/2))/27 + 272/27)^(1/2))/2
vpasolutions =
0
-p*2.0i
p*2.0i
-p*0.66666666666666666666666666666667i
p*0.66666666666666666666666666666667i
-2.1575776918969228428670635119892*p
2.1575776918969228428670635119892*p
In my understanding, this should not happen since before I made an assumption about p, MATLAB should give me the general solution that is valid for all p (even the positive ones). Now that I made an assumption, the number of solutions should not change.
Can you spot a mistake? Did I not understand the behavior of 'solve' correctly? Or is the symbolic solver really faulty?
2 个评论
Alex Mcaulley
2019-6-18
编辑:Alex Mcaulley
2019-6-18
Really interesting!
syms p t
formula = 84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 + (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
solutions = solve(formula == 0, t)
solutions =
0
-p*2i
p*2i
-(p*2i)/3
(p*2i)/3
-(p*(272/27 - (64*13^(1/2))/27)^(1/2))/2
(p*(272/27 - (64*13^(1/2))/27)^(1/2))/2
-(p*((64*13^(1/2))/27 + 272/27)^(1/2))/2
(p*((64*13^(1/2))/27 + 272/27)^(1/2))/2
solutions 6 and 7 are not solutions in fact:
simplify(subs(formula,t,solutions(7)))
ans =
-(32768*p^6*(25*13^(1/2) - 86))/729 %This is only solution for p == 0
And then, when you assume that p can not be 0, this solution disappears (because for p ~= 0 is not a solution)
John D'Errico
2019-6-18
Well, they are technically solutions, if you recognize that terms like
(4*p^2 + t^2)^(1/2)
have TWO branches.
There is probably a way to tell MATLAB to follow only the positive branch.
回答(1 个)
John D'Errico
2019-6-18
编辑:John D'Errico
2019-6-18
It seems a subtle thing. We need to "branch" out from our preconceptions. ;-)
pretty(solutions)
/ 0 \
| |
| -p 2i |
| |
| p 2i |
| |
| p 2i |
| - ---- |
| 3 |
| |
| p 2i |
| ---- |
| 3 |
| |
| -#1 |
| |
| #1 |
| |
| -#2 |
| |
\ #2 /
where
/ 272 64 sqrt(13) \
p sqrt| --- - ----------- |
\ 27 27 /
#1 == ---------------------------
2
/ 64 sqrt(13) 272 \
p sqrt| ----------- + --- |
\ 27 27 /
#2 == ---------------------------
2
So #1 contains the questionable solutions. Are they? Are they spurious? It looks like case #1 converts to case #2, if we take the negative branch of the sqrt for sqrt(13).
It looks like solutions [1 2 3 4 5 8 9] all trivially satisfy the original problem, with 6 and 7 a problem.
simplify(subs(formula,t,solutions))
ans =
0
0
0
0
0
-(32768*p^6*(25*13^(1/2) - 86))/729
-(32768*p^6*(25*13^(1/2) - 86))/729
0
0
But suppose we look carefully at formula?
formula
formula =
84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 + (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
In there, I see two square roots. When you take a square root, you can take the positive or the negative branch. Either is equally valid. So for any values of p and t, if I write Q as
Q = (4*p^2 + t^2)^(1/2)
then we might have intended -Q also, and that expression raised to the 1/2 power becomes ambiguous.
Let me now change formula, swapping the sign in front of the square rooted terms. This allows either of those 1/2 powers to act as if we took the negative branch.
>> formulamod = 84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 - (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
formulamod =
84*p^2*t^4 - 112*p^4*t^2 - 64*p^6 + 27*t^6 - (4*p^2 + t^2)^(1/2)*(4*p^2 + 9*t^2)^(1/2)*(16*p^4 - 88*p^2*t^2 + 9*t^4)
Now subs in the solutions.
>> simplify(subs(formulamod,t,solutions))
ans =
-128*p^6
0
0
0
0
0
0
(32768*p^6*(25*13^(1/2) + 86))/729
(32768*p^6*(25*13^(1/2) + 86))/729
So solutions #6 and #7 trivially solve formulamod. Effectively, solutions #6 and #7 solve the problem where EXACTLY one of the terms in formula (where a 1/2 power was taken) follows the negative branch of the square root. So all 9 solutions are technically solutions to the problem, as long as either branch in those squre roots is acceptable.
Does this get to the "root" of the matter? Perhaps. Don't just "woodenly" think a square root is always a positive number. Should I just "leaf" things alone? ;-) (Sorry about that.)
4 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Assumptions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!