How to solve a coupled system of transcendental equations?

Question

If I try to Solve the following system of three coupled equations, it does not yield solutions, and but gives an error

eqns = {
   x - Exp[-(1 - x - y - z)] == 0,
   y - Exp[-(1 - 0.5 x - y - z)] == 0,
   z - Exp[-(1 - x - 0.8 y - z)] == 0
   };


Solve[ eqns, {x, y, z}]

Solve::nsmet: This system cannot be solved with the methods available to Solve.

How can I solve the following system of equations?

Answer 1

It doesn't look like there are any real solutions.

eqns = {x - Exp[-(1 - x - y - z)] == 0, 
   y - Exp[-(1 - 0.5*x - y - z)] == 0, 
   z - Exp[-(1 - x - 0.8*y - z)] == 0};

ContourPlot3D[Evaluate@eqns, {x, -1.5, 1.5}, {y, -1.5, 1.5}, {z, -1.5, 1.5}, 
  ViewVector -> {5, 5, 5}]

Edit: Of course there may well be complex ones.

Answer 2

There are infinitely many solutions. Here we demonstrate only a small part of the x - solution space in the complex plane.

Let's start from the begining. There are no real solutions:

Reduce[{ x - Exp[-(1 - x - y - z)] == 0, 
         y - Exp[-(1 - x/2 - y - z)] == 0, 
         z - Exp[-(1 - x - 4/5 y - z)] == 0}, {x, y, z}, Reals]

False

Without domain specification (it means we are looking for complex solutions) Reduce has to be supplemented by an appropriate domain specification, otherwise it won't find any solutions.

To demonstrate what is going on here let's take a pedestrian approach, after all it might be automized a bit.

From the first equation we see that $x=a$ where we introduced $\quad a = \exp(-1+x+y+z)$. Now we get from the second equation $y=a \exp(-x/2)$ and in turn from the first one we have $y=x \exp(-x/2) $. From the third equation we get $z=a\exp(-y/5)$. Next we can obtain an equation for $x$ by eliminating $y$ and $z$ i.e. it yields $$x=e^{-1+x+x e^{-\frac{x}{2}} + x e^{-\frac{x}{5} e^{-\frac{x}{2}}}}$$ We need to restrict appropriately the complex domain since knowing properties of exponential functions we expect there are infinitely many solutions. Now we obtain two exact solutions for x nearest zero in the complex plane:

Reduce[ x == Exp[-1 + x (1 + Exp[-x/2] + Exp[-(x/5) Exp[-(x/2)]])] && 
        Abs[x] < 2, x]
sx = {ToRules @ % } // N;

 x == Root[{-E^(
   E^(-(#1/2) - 
     1/5 E^(-(#1/2)) #1) (-E^(#1/2 + 1/5 E^(-(#1/2)) #1) + 
      E^(#1/2) #1 + E^(1/5 E^(-(#1/2)) #1) #1 + 
      E^(#1/2 + 1/5 E^(-(#1/2)) #1) #1)) + #1 &, 
0.055803204968110806464128827202522530324263607292 - 
 0.510206373369089749181274412607230169611834407013 I}] || 
x == Root[{-E^(
   E^(-(#1/2) - 
     1/5 E^(-(#1/2)) #1) (-E^(#1/2 + 1/5 E^(-(#1/2)) #1) + 
      E^(#1/2) #1 + E^(1/5 E^(-(#1/2)) #1) #1 + 
      E^(#1/2 + 1/5 E^(-(#1/2)) #1) #1)) + #1 &, 
0.055803204968110806464128827202522530324263607292 + 
 0.510206373369089749181274412607230169611834407013 I}]

Root is a symbolic representation of exact solutions, see e.g. How do I work with Root objects?. For example taking restriction Abs[x] < 370/100 after a few minutes I got 102 complex solutions for x. One can show that if x + I y is a solution then x - I y is a solution too (for x and y reals)

To demonstrate the structure of the x-solution space in a small part of the complex plane we define:

f[x_, y_] := -E^(-1 + (1 + E^(1/2 (-x - I y)) + 
   E^(-(1/5) E^(1/2 (-x - I y)) (x + I y))) (x + I y)) + x + I y

and we can catch a glimpse of multitude of solutions

sol = {ToRules@
Reduce[x == 
   Exp[-1 + x (1 + Exp[-x/2] + Exp[-(x/5) Exp[-(x/2)]])] && 
  Abs[x] < 370/100, x]};
Length[sol]

102

ContourPlot[{Re[f[x, y]] == 0, Im[f[x, y]] == 0},
            {x, -3.7, -2.85}, {y, 0.65, 1.5}, 
  PlotPoints -> 50, ImageSize -> Large, 
  ContourStyle -> {Orange, Darker@Cyan}, 
  Epilog -> {Red, PointSize[0.015], Point[ReIm[x /. sol]]}]

Solutions are denoted by red points where curves intersect. For a technical reason not all the intersection point are drawn.

The results for y and z might be expressed also in the symbolic form, nontheless for brevity we demonstrate them only related to given two solutions for x in the numerical form

sy = y -> x Exp[-x/2] /. sx

{y -> 0.177717 - 0.466416 I, y -> 0.177717 + 0.466416 I}

z -> x  Exp[-y/5] /. Transpose[{sx // Flatten, sy}]

{z -> 0.0994856 - 0.485233 I, z -> 0.0994856 + 0.485233 I}

Reduce couldn't find solutions without approprate restriction of domains for x, y, z separately. This is crucial point to get any solutions.

Answer 3

Edit: As pointed out by @anderstood, FindRoot can't be used to decisively say if there are no solutions for a set a equations as it depends on the starting values of variables.

You can use FindRoot to solve these transcendental equations. Real and complex solutions are checked, if any, as follows.

Checking for Real solutions:

exps = {x - Exp[-(1 - x - y - z)], y - Exp[-(1 - 0.5 x - y - z)], 
   z - Exp[-(1 - x - 0.8 y - z)]};
sols = FindRoot[Thread[exps == 0], {x, 0}, {y, 0}, {z, 0}] (*FindRoot throws a warning message*)

{x -> -0.00592577, y -> -0.0047965, z -> -0.00555911}

Substituting the above solution in exps gives,

exps /. sols

{-0.367864, -0.367809, -0.367845}

Checking for Complex solutions: To look for complex solutions, add 1.I to the initial values of variables.

exps = {x - Exp[-(1 - x - y - z)], y - Exp[-(1 - 0.5 x - y - z)], 
   z - Exp[-(1 - x - 0.8 y - z)]};
sols = FindRoot[
  Thread[exps == 0], {x, 0 + 1. I}, {y, 0 + 1. I}, {z, 0 + 1. I}]

{x -> 0.0558032 + 0.510206 I, y -> 0.177717 + 0.466416 I, z -> 0.0994856 + 0.485233 I}

Substituting the above solution into exps gives,

exps /. sols

{8.32667*10^-17 + 1.11022*10^-16 I, 2.77556*10^-17 + 1.11022*10^-16 I, 4.16334*10^-17 + 5.55112*10^-17 I}

Looking at these solutions, it can be seen that the equations only have a complex solution.