Solve[ ] likes extra variables

In[1]:= (* This problem was posted on LinkedIn by Songxin (Derek)Liang
16^(x^2-y)+16^(y^2-x) == 1 *)

In[2]:= (* Solve and Reduce cannot get a solution *)
Solve[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, y}]

During evaluation of In[2]:= Solve::nsmet: This system cannot be solved with the methods available to Solve.

Out[2]= Solve[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}]

In[3]:= Solve[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, y}, Reals]

During evaluation of In[3]:= Solve::nsmet: This system cannot be solved with the methods available to Solve.

Out[3]= Solve[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}, Reals]

In[4]:= Reduce[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, y}]

During evaluation of In[4]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce.

Out[4]= Reduce[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}]

In[5]:= Reduce[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}, Reals]

During evaluation of In[5]:= Reduce::nsmet: This system cannot be solved with the methods available to Reduce.

Out[5]= Reduce[16^(x^2 - y) + 16^(-x + y^2) == 16, {x, y}, Reals]

In[6]:=  (*  Adding two extra variables to the problem gives a very \
complex solution (not show here) which contains b as a parameter *)

In[7]:= sln = 
  Solve[{x^2 - y == a, y^2 - x == b, 16^a + 16^b == 1}, {x, y, a, b}];

During evaluation of In[7]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[7]:= Solve::svars: Equations may not give solutions for all "solve" variables.

In[8]:=  (* Setting b = -1/4 gives the solution when the two terms \
are equal  since 16^(-1/4) = 1/2*)
sln1 = sln /. b -> -1/4 // N

Out[8]= {{x -> -0.5 + 1. I, y -> -0.5 - 1. I, 
  a -> -0.25}, {x -> -0.5 - 1. I, y -> -0.5 + 1. I, 
  a -> -0.25}, {x -> 0.5, y -> 0.5, a -> -0.25}, {x -> 0.5, y -> 0.5, 
  a -> -0.25}}

In[9]:=  (* Testing shows the solutions are correct*)
16^(x^2 - y) + 16^(-x + y^2) - 1 /. sln1


Out[9]= {4.44089*10^-16 + 0. I, 4.44089*10^-16 + 0. I, 6.66134*10^-16,
  5.55112*10^-16}

In[10]:=  (* Setting b = -1 gives another solution *)
sln2 = sln /. b -> -1 // N

Out[10]= {{x -> -0.719587 + 0.942605 I, y -> -0.347422 - 1.35657 I, 
  a -> -0.0232774}, {x -> -0.719587 - 0.942605 I, 
  y -> -0.347422 + 1.35657 I, 
  a -> -0.0232774}, {x -> 0.719587 - 0.440069 I, 
  y -> 0.347422 - 0.633336 I, 
  a -> -0.0232774}, {x -> 0.719587 + 0.440069 I, 
  y -> 0.347422 + 0.633336 I, a -> -0.0232774}}

In[11]:=  (* Testing shows these solutions are also correct*)
16^(x^2 - y) + 16^(-x + y^2) - 1 /. sln2


Out[11]= {-1.66533*10^-16 + 6.34877*10^-16 I, 
 5.55112*10^-17 + 6.34877*10^-16 I, 5.55112*10^-17 - 2.59722*10^-16 I,
  1.17961*10^-16 - 2.69342*10^-16 I}
1 Like

Interesting. However, the equation is underdetermined. You get individual solutions by adding constraints, for example

Solve[16^(x^2 - y) + 16^(y^2 - x) == 16 && x == 0 && Abs[y] < 2,
 {x, y}]
Solve[16^(x^2 - y) + 16^(y^2 - x) == 16 && x + 2 y == 0 &&
  Abs[y] < 2, {x, y}]

On my computer I can get one solution in 43 mSec by your method and all the solutions in 208 mSec. by my method

In[18]:= AbsoluteTiming[
 Solve[16^(x^2 - y) + 16^(y^2 - x) == 16 && x == 0 && Abs[y] < 2, {x, 
    y}];]

Out[18]= {0.0426189, Null}

In[17]:= AbsoluteTiming[
 sln = Solve[{x^2 - y == a, y^2 - x == b, 16^a + 16^b == 1}, {x, y, a,
      b}];]

During evaluation of In[17]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[17]:= Solve::svars: Equations may not give solutions for all "solve" variables.

Out[17]= {0.207841, Null}

How do you know that your method gives ALL complex solutions? It seems that it is based on a partial solution of

Solve[16^a + 16^b == 16, b]

which involves taking one single branch of the inverse of the complex exponential function. You may miss any solution on the other branches.

Yes, you are probably correct. I tried using Reduce rather than Solve but never got a solution in a reasonable amount of time. However, setting b equal to -1/4 gave all the solutions posted in LinkedIn, replying to the question.

This two-steps variation has a chance of getting all solutions:

Solve[{x^2 - y == a, y^2 - x == b}, {x, y}] /. 
 Solve[16^a + 16^b == 16, a][[1]]

in terms of the continuous parameter b and the discrete parameter C[1].

Yes, Solve sometimes can use some human help.

Here is a way to recover the real solutions:

Clear[x1, x2, x3, x4, y1, y2, y3, y4, x, y, a, b, n];
{{x1[b_, n_], y1[b_, n_]}, {x2[b_, n_], y2[b_, n_]}, {x3[b_, n_], 
    y3[b_, n_]}, {x4[b_, n_], 
    y4[b_, n_]}} = {x, y} /. 
     Solve[{x^2 - y == a, y^2 - x == b}, {x, y}] /. 
    Solve[16^a + 16^b == 16, a][[1]] /. C[1] -> n;
ParametricPlot[{Chop@{x3[b, 0], y3[b, 0]}, 
  Chop@{x4[b, 0], y4[b, 0]}}, {b, -2, 2}]

The result can be checked against this:

ContourPlot[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, -1, 2}, {y, -1, 2}]
1 Like

The right-hand side is 1, not 16. In my first post I had 16 for the Solve and Reduce tests that didn’t work, but had 1 as the right-hand side in the solution method I proposed. Sorry about that.

If the right-hand side is 1 there seems to exist only one real solution:

Clear[x1, x2, x3, x4, y1, y2, y3, y4, x, y, a, b, n];
{{x1[b_, n_], y1[b_, n_]}, {x2[b_, n_], y2[b_, n_]}, {x3[b_, n_], 
    y3[b_, n_]}, {x4[b_, n_], 
    y4[b_, n_]}} = {x, y} /. 
     Solve[{x^2 - y == a, y^2 - x == b}, {x, y}] /. 
    Solve[16^a + 16^b == 1, a][[1]] /. C[1] -> n;
{x3[-1/4, 0], y3[-1/4, 0]} // FullSimplify

More than merely adding variables, maybe Mathematica likes the symmetry structure that you introduced with a and b.

I wonder if there are other problems where this approach works.

I found the real solution using Minimize

Minimize[(16^(x^2 - y) + 16^(y^2 - x) - 1)^2, {x, y}]

{0, {x -> 1/2, y -> 1/2}}

Hmmm, why doesn’t this work?

Minimize[16^(x^2 - y) + 16^(y^2 - x), {x, y}]

NMinimize does, though.

I minimized the square of the difference of the left-hand side and right-hand side, since the square should be a minimum of zero when the equation is satisfied. Why should the minimum of the difference of the first power of the difference gave the correct solution, in general?

In this particular case we know it is the same. In general, the necessary condition for a minimum of (f[x, y] - 1)^2 is implied by the necessary condition for a minimum of f[x, y]:

D[(f[x, y] - 1)^2, {{x, y}}] == {0, 0}
D[f[x, y], {{x, y}}] == {0, 0}

I wonder why Mathematica can solve one and not the other.

Another example of this technique

In[1]:= (* Consider the equation Sin[x^2-y^2]+ Cos[x^2+y^2]==1,
with x and y between 0 and 1.
It is one equation with two variables 
and is clearly underdetermined.  
An obvious solution is x -> 0 and y -> 0,             since Sin[0]==0 \
and Cos[0]= 1.    
However Solve[] cannot find a solution *)

In[2]:= Solve[{Sin[x^2 - y^2] + Cos[x^2 + y^2] == 1, 0 <= x <= 1, 
  0 <= y <= 1}, {x, y}]

During evaluation of In[2]:= Solve::nsmet: This system cannot be solved with the methods available to Solve.

Out[2]= Solve[{Cos[x^2 + y^2] + Sin[x^2 - y^2] == 1, 0 <= x <= 1, 
  0 <= y <= 1}, {x, y}]

In[3]:= Solve[{Sin[x^2 - y^2] + Cos[x^2 + y^2] == 1, 0 <= x <= 1, 
  0 <= y <= 1}, {x, y}, Reals]

During evaluation of In[3]:= Solve::nsmet: This system cannot be solved with the methods available to Solve.

Out[3]= Solve[{Cos[x^2 + y^2] + Sin[x^2 - y^2] == 1, 0 <= x <= 1, 
  0 <= y <= 1}, {x, y}, Reals]

In[4]:= (* Mathematica can find a solution 
by minimizing the square of the difference
 of the left-hand side and the right-hand side *)
slnm = Minimize[{(Sin[x^2 - y^2] + Cos[x^2 + y^2] - 1)^2, 0 <= x <= 1,
    0 <= y <= 1}, {x, y}]

Out[4]= {0, {x -> Sqrt[\[Pi]]/2, y -> Sqrt[\[Pi]/3]/2}}

In[5]:=  (* Extra variables are added giving 24 solutions 
with a as a parameter *)
Length[slnab = 
  Solve[{Sin[a] + Cos[b] == 1, a == x^2 - y^2, b == x^2 + y^2}, {a, b,
     x, y}]]

During evaluation of In[5]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[5]:= Solve::svars: Equations may not give solutions for all "solve" variables.

Out[5]= 24

In[6]:= (* Setting a -> 0 gives the solutions for x and y equal 0 *)
Union[slnab /. a -> 0]

Out[6]= {{b -> 0, x -> 0, y -> 0}}

In[7]:= (* What about the solution found by Minimize ?
First calculate a *)

In[8]:= x^2 - y^2 /. slnm[[2]]

Out[8]= \[Pi]/6

In[9]:= (* The 12th solution gives the result found by Minimize *)
(slnab /. a -> \[Pi]/6)[[12]]

Out[9]= {b -> \[Pi]/3, x -> Sqrt[\[Pi]]/2, y -> Sqrt[\[Pi]/3]/2}

In[10]:= slnab[[12]]

Out[10]= {b -> ArcCos[1 - Sin[a]], 
 x -> Sqrt[ArcCos[1 - Sin[a]] + ArcSin[Sin[a]]]/Sqrt[2], 
 y -> Sqrt[ArcCos[1 - Sin[a]] - ArcSin[Sin[a]]]/Sqrt[2]}

In[11]:= (* This is the solution *)
(Sin[x^2 - y^2] + Cos[x^2 + y^2]) /. slnab[[12]] // FullSimplify

Out[11]= 1