Symbolic solver yields wrong solution

5 views (last 30 days)
Friedrich Tuttas on 31 May 2019
Edited: Friedrich Tuttas on 31 Jul 2019
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?
John D'Errico on 18 Jun 2019
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.

John D'Errico on 18 Jun 2019
Edited: John D'Errico on 18 Jun 2019
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.)
Friedrich Tuttas on 31 Jul 2019
It has been a while but I just found this video which might clarify things for you:
The square root is not ambiguous!

R2019a

Community Treasure Hunt

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

Start Hunting!