Symbolic solver yields wrong solution

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.

(see https://www.wolframalpha.com/input/?i=solve+84*p%5E2*t%5E4+-+112*p%5E4*t%5E2+-+64*p%5E6+%2B+27*t%5E6+%2B+(4*p%5E2+%2B+t%5E2)%5E(1%2F2)*(4*p%5E2+%2B+9*t%5E2)%5E(1%2F2)*(16*p%5E4+-+88*p%5E2*t%5E2+%2B+9*t%5E4)+%3D+0)

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?

[Alex Mcaulley]

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]

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.

Answer 1

[John D'Errico]

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]

Thanks for the answer but my understanding of math as I learned it is a little different.

This is how I understand the algebra you are trying to describe:

The equation always has two (complex) solutions for x for any given a. Those are and .

However, the equation does not have two solutions for x. It is always the positive "branch" as you would say.

An example:

has the two solutions for x: and . However, the formula always equates to and never !

Now, as I said, this is the algebra that stems from the definitions I was taught at school. It is always possible to make different definitions and see where that takes on. If the results are helpful, one might very well keep those definitions. In this case, one could define that the square root function does indeed have two branches. I have, however, never heard of such a definition.

I just checked wikipedia and discovered that the page about square roots disagrees with itself within the first few lines about that topic by saying "Every nonnegative real number a has a unique nonnegative square root ..." and "Every positive number a has two square roots..." (https://en.wikipedia.org/wiki/Square_root).

The problem I have with your two-branch definition of the square root is that the square root becomes unnecessarily ambiguous. I do not see a need to have it this way. Furthermore, if you plug in a supposed solution to an equation into that equation, it should yield zero, which does not happen here, since the two definitions of the square root are mixed for some reason. It uses your definition for determining solutions but mine for all the arithmetics.

But even if MATLAB used your definition for my formula, it does not do this consistently. Since the square root would have to be ambiguous, MATLAB would have to give me two solutions every time I use it. For example, the following line should yield two results but only yields one:

9^(1/2)

yields only the positive result:

ans =

     3

Even within the symbolic solver, your definition is not used consistently. Take for example the following code:

syms x
formula = x == 9^(1/2);
solve(formula,x)

That gives me:

ans =
 
3

Can you explain to me why it does not consider the negative branch in this example?

Moreover, none of this explains why the ominous solutions disappear when the assumptions is made that p is positive (as described in my first post).

[John D'Errico]

That does not mean your understanding of mathematics is complete.

"In this case, one could define that the square root function does indeed have two branches. I have, however, never heard of such a definition. "

You certainly need to take a complex analysis course if you have never even read about these concepts. At least, if you will be making arguments about mathematics, it helps to have a decent background.

The function y = sqrt(x) is a single valued thing, as it must be, otherwise it is not useful as a function. So we agree that sqrt(x) as a function returns only one value - the principal value.

Anyway, it ia a good idea to test the solutions returned, to verify they are truly meaningful and valid. Spurious solutions can arise. Here, that happens because there is some ambiguity that can creep into an implicit function with those fractional powers. There are multiple places where a square root was apparently necessary in the computation of that solution.

Unfortunately, I cannot trace through the logical path the symbolic toolbox must have taken, to know what was done, and if there was something done that created invalid solutions. All that I can do is as I did. That is, I showed that IF we assume the terms of the form:

(4*p^2 + t^2)^(1/2)

actually have two branches, then solutions #6 and #7 naturally arise, and are solutions. At the same time, they fall into a gray area.

You can feel free to decide they are not valid by your choices, that only the positive branch of the sqrt exists in your world view. That world will be a chaste, limited one.

As for why those solutions disappear when p was constrained to be positive? I'd need to take a guess at that. Constraining p to be a positive number does more than just exclude negative numbers! It also constrains p to be a REAL number. That is just a guess though, without much thought invested.

[Friedrich Tuttas]

I do not quite understand what you mean. You say "The function y = sqrt(x) is a single valued thing..." but also "... there is some ambiguity that can creep into an implicit function with those fractional powers."

To me those two statements disagree. Only an ambiguous function (which, according to the both of us, the square root is not) can lead to ambiguity.

Are you maybe trying to say the following? "While solving an equation, one might take steps that create solutions that are not actually solutions of the original problem." I completely agree with that. A simple example would be: (1)

The first idea could be to square the entire equation, which leads to: (2)

Now this equation has the two solutions: and

However, if we insert these in the original problem we get:

, which is correct, and

, which is false.

I say the latter is false since I used that the square root function "is a single valued thing" and always only yields the positive "branch". We get a false solution here because we squared the equation in the first step. Squaring an equation like (1) as well as multiplying it by a function of x add solutions that do solve the new equation (2) but not necessarily solve (1). This is still a valid way to solve an equation because oftentimes the new equation is easier to understand. One just has to check the results in the end and remove the ones that were artificially added.

If you are saying that MATLAB might have chosen a similar way to solve my equation, that still means that it forgot to check the solutions.

I want to sum this up:

I had an equation and asked MATLAB to find the solutions for t for which this equation is satisfied. Among others, it gave me the solutions and for which MATLAB itself agreed that and . If I ask MATLAB to show me a t that satisfies and it gives me , for which , there is something fishy going on. This, to me, ultimately shows that and are not solutions to my original problem. I have presented the following explanations:

1) MATLAB might use different definitions for the square root function. It might define it as the ambiguous two-branch function within the symbolic solver but defines it as the unambiguous principal branch function for all the arithmetic while trying to figure out what and yield. I have started to disagree with this explanation because, if that were the case, the symbolic solver would have to spit out two solutions every time one uses a square root (which it does not, as can be seen in the example of my last post).

2) and might be solutions to an equation that was created (by squaring or multiplying by a function of x or yet another non-equivalence-transformation) while trying to solve my problem. But MATLAB forgot to check whether they actually solve the original problem.

As for the other problem with the constraint on p. As I said in my starting post: Without any assumption, MATLAB should give me the general solution that is valid for any p. If I now constrain p to the positive (real) numbers, all solutions should still be valid.

EDIT:

I tried the following code

syms x
formula = x - sqrt(2*x + 3)
solve(formula==0,x)

and got the following result

formula =
 
x - (2*x + 3)^(1/2)
 
 
ans =
 
3

So MATLAB seems to agree that there is only one solution to the equation (1). Now, if I turn the 3 in the equation into a symbolic varbiable y, the following happens:

syms x y
formula = x - sqrt(2*x + y)
solve(formula==0,x)

yields

formula =
 
x - (2*x + y)^(1/2)
 
 
ans =
 
 (y + 1)^(1/2) + 1
 1 - (y + 1)^(1/2)

Two solutions! So for an undefined y there are apparently more solutions than for . This should not happen! I should be able to insert into the two solutions and get all solutions for this case. But instead of (as in the first code snippet) I get and ! Here, MATLAB contradicts itself.

[Friedrich Tuttas]

It has been a while but I just found this video which might clarify things for you:

https://www.youtube.com/watch?v=gIkoX8WRH3c

The square root is not ambiguous!