How to symbolically solve a coupled system of equations with Mathematica?

Question

The symbolic solvers of Mathematica are strong, however, the ones fail in some cases, for example, for the system {x^(x + y) == y^12, y^(x + y) == x^3} over the reals. Indeed,

Solve[{x^(x + y) == y^12, y^(x + y) == x^3}, {x, y}, Reals]

Solve::nsmet: This system cannot be solved with the methods available to Solve. Try Reduce or FindInstance instead.

and

Reduce[{x^(x + y) == y^12, y^(x + y) == x^3}, {x, y}, Reals]

and

FindInstance[{x^(x + y) == y^12, y^(x + y) == x^3}, {x, y}, Reals]

are running without any response for hours in 14.1 on Windows 10.

The plot

ContourPlot[{x^(x + y) == y^12, y^(x + y) == x^3}, {x, -2, 5}, {y, -2,  5}]

suggets the solutions {x -> 1, y -> 1} and {x -> 4, y -> 2} , but the plot does not show a stand alone solution {x -> 1, y -> -1}. All three solutions can be found numerically by

NSolve[{x^(x+y)==y^12,y^(x+y)==x^3,Abs[x]<=5,Abs[y]<=5},{x,y},Complexes]

{{x -> 0.633359 - 0.238871 I, y -> 0.468235 - 0.809188 I}, {x -> 4.87806 + 0.0317266 I, y -> 1.80523 - 1.61064 I}, {x -> 1., y -> 1.}, {x -> 4.87806 - 0.0317266 I, y -> 1.80523 + 1.61064 I}, {x -> 0.633359 + 0.238871 I, y -> 0.468235 + 0.809188 I}, {x -> 0.536597 + 0.0284836 I, y -> -0.455244 + 0.88115 I}, {x -> 1., y -> -1.}, {x -> 2.58416 - 2.28479 I, y -> -0.776203 + 0.780143 I}, {x -> -0.0133607 + 0.0251473 I, y -> -3.55574 - 2.17787 I}, {x -> 4., y -> 2.}, {x -> 0.816233 + 0.244729 I, y -> 0.81747 + 0.50487 I}, {x -> 4.57655 - 1.68878 I, y -> 1.42345 + 1.68878 I}, {x -> 0.542606 + 0.153983 I, y -> 0.037173 + 0.947603 I}, {x -> 0.693224 - 0.106439 I, y -> -0.870088 + 0.513775 I}, {x -> 2.86014 - 1.25472 I, y -> -1.1218 + 0.18289 I}, {x -> 2.86014 + 1.25472 I, y -> -1.1218 - 0.18289 I}, {x -> 0.693224 + 0.106439 I, y -> -0.870088 - 0.513775 I}, {x -> 0.542606 - 0.153983 I, y -> 0.037173 - 0.947603 I}, {x -> 4.57655 + 1.68878 I, y -> 1.42345 - 1.68878 I}, {x -> 0.816233 - 0.244729 I, y -> 0.81747 - 0.50487 I}, {x -> -0.382865 + 0.00701507 I, y -> -1.0598 - 0.987606 I}, {x -> 4.20838 - 2.08547 I, y -> 2.07887 + 0.0389557 I}, {x -> -3.3177 + 1.10786 I, y -> -0.412412 - 0.398653 I}}

, but there is a chance that there exist other sporadic real solutions. So how to symbolically solve that system with Mathematica?

Answer 1

The system is transcendental, not algebraic and this is the first point why we should take a closer look at the system. It is reasonable to rewrite the system to a more comprehensible form. $$x^{(x+y)}=\exp\left((x+y)\log\left( x\right) \right)$$ Prescribing a condition $x>0\;$ we encounter the problem that $\log\left(x \right)\;$ is not bounded from below, moreover Log[z] has a branch cut discontinuity in the complex z plane running from -Infinity to 0. Taking the above into account we can easily find

Reduce[{ Exp[(x + y) Log[x]] == y^12, Exp[(x + y) Log[y]] == x^3, 
         10^-3 <= x <= 5, 10^-3 <= y <= 5}, {x, y}]

 (x == 1 && y == 1) || (x == 4 && y == 2)

Setting e.g. 10^-6 as a lower bound for x we might expect a longer symbolic processing. We need not point out the expected domain Reals unless for negative arguments Log becomes complex.

Edit

Searching for larger class of solutions we could transform appropriately the original system to another form in which it might appear easier to solve. At first sight we couldn't find a more convenient form, noetheless there might be an especially suitable transformation of variables. Another aspect is exploiting appropriate tools in Mathematica. Our system of equations appears to be very sensitive on varous ways we could use powerful symbolic solvers.

Reduce[{ Exp[(x + y) Log[x]] - y^12 == 0, 
         Exp[(x + y) Log[y]] - x^3 == 0, 
         -10 <= x <= 10, -10 <= y <= -1/10}, 
       {x, y}, Reals, Cubics -> True] // TraditionalForm

N @ %

( x == 1. && y == -1.) || ( x == 9. && y == -3.) || 
( x == 0.0275247 && y == -6.02752)

More diligent investigation of the system can provide further real solutions, nonetheless we need a kind of mathematical reasoning for deciding how many real solutions exist. Let's demostrate contours of appropriate functions near the irrational solution.

ContourPlot[{Re[Exp[(x + y) Log[x]] - y^12] == 0, 
  Re[Exp[(x + y) Log[y]] - x^3] == 0, Im[Exp[(x + y) Log[y]] - x^3] == 0}, 
  {x, -1/20, 1/20}, {y, -5, -12}, 
  Evaluated -> True, PlotPoints -> 100, MaxRecursion -> 5, 
  PlotLegends -> "Expressions", ContourStyle -> Thick, 
  Epilog -> {Red, PointSize[0.02], Point[{0.02752, -6.02752}]}]

Another problem remaining is a constructive searching for complex symbolic solutions. We expect there are infinitely many complex solutions however it appears a hard task to find even a few of them. Nonetheless we can find a simple example:

With[{ x = I Exp[I Pi/6], y = -I Exp[I Pi/6]}, 
  { x^(x + y) - y^12, y^(x + y) - x^3} // Simplify]

{0, 0}