# Finding correct solution to this equation

Consider the equation

eq = x == y - Log[1 + y];


I would like to solve this equation for y

yofx = y /. Solve[eq, y]


{-1 - ProductLog[-E^(-1 - x)]}

Unfortunately, this solution seems to be incorrect, as can be verified by inspecting the following two plots

{Plot[y - Log[1 + y], {y, 0, 10}], Plot[yofx, {x, 0, 10}]}


On the left, the plot correctly rises and asymptotes to a straight line for large y and x. On the right, the plot should have been a flipped image of the left plot, with x- and y-axes interchanged, but the curve looks nothing like it.

Is it possible to get the correct solution for this equation in Mathematica? If so, how can it be done?

• I get a warning from Solve that the answer may be wrong. Why not follow the advice in the message? Apr 12 at 3:22
• @MichaelE2 Using Reduce produces an output that still contains an unsolved equation for a constant specifying a branch cut. Instead of trying to solve that equation, I just found by trial and error that the correct branch constant value is -1, since the constant must be an integer. Thanks for the suggestion! Apr 12 at 3:29
• If you use Reduce[x == y - Log[1 + y] , y] you'll get a generalized solution. Since it seems you are interested in positive values of x and y, you could try Reduce[x == y - Log[1 + y] && x > 0 && y > 0, y], which narrows your solution to y==-1 - ProductLog[-1, -E^(-1 - x)]. If you then plot it using Plot[-1 - ProductLog[-1, -E^(-1 - x)], {x, 0, 10}], you'll get what looks like the axes-flipped version of your left plot. Apr 12 at 3:30
• @theorist Thank you, that's even better! Apr 12 at 3:32
• It's worth noting that the inverse function Mathematica provided is the inverse of $y - \ln(1 + y)$ restricted to the domain $(-1, 0]$ rather than $[0, \infty)$. So it's not wrong per se, it's just that Mathematica made an assumption different from yours. Apr 12 at 15:49

To plot an equation use ContourPlot or, in this case, ParametricPlot

eq = x == y - Log[1 + y];

ContourPlot[Evaluate@eq, {x, 0, 10}, {y, 0, 10},
FrameLabel -> (Style[#, 16] & /@ {x, y})]


ParametricPlot[{eq[[-1]], y}, {y, 0, 10},
Frame -> True,
PlotRange -> {{0, 10}, {0, 10}},
FrameLabel -> (Style[#, 16] & /@ {x, y})]


(* same plot *)

yofx = SolveValues[{eq, y > 0}, y, Reals][[1]] // Quiet

Plot[Evaluate@yofx, {x, 0, 10},
Frame -> True,
PlotRange -> {{0, 10}, {0, 10}},
FrameLabel -> (Style[#, 16] & /@ {x, y}),
AspectRatio -> 1]


(* same plot *)

• Adding a condition y > 0 and requirement Reals to Solve is a very nice way to get the result directly, without looking at Reduce! Thank you! Apr 12 at 3:48

This is a typical case where conditional and branch cut check routines are not sophisticated enough to get result, although you give all neccesary conditions (x>0 and y>0 or y<0).

But Solve knows the underlying rules, as can be seen, when you insert numbers for variables (here EulerGamma for x). Of course you have to check validity of that trick, e.g. with a graphic or insert solution in eq as shown below.

eq = x == y - Log[1 + y];

ContourPlot[Evaluate@eq, {x, 0, 6}, {y, -3, 3},
GridLines -> {{EulerGamma}, {-1 -
ProductLog[-E^(-1 - EulerGamma)], -1 -
ProductLog[-1, -E^(-1 - EulerGamma)]}}]


{sol1 = Solve[eq /. x -> EulerGamma, y, Reals], sol1 // N}

(*   {{{y -> -1 - ProductLog[-E^(-1 - EulerGamma)]}, {y -> -1 -
ProductLog[-1, -E^(-1 - EulerGamma)]}}, {{y -> -0.729215}, {y ->
1.48916}}}   *)


Get the two solutions

{y1[x_], y2[x_]} = y /. sol1 /. EulerGamma -> x

(*   {-1 - ProductLog[-E^(-1 - x)], -1 - ProductLog[-1, -E^(-1 - x)]}   *)


Test it applying Exp to both sides of equations

Map[Exp[#] &, {eq /. y -> y1[x], eq /. y -> y2[x]}, {2}] //
FullSimplify[#, x > 0] &

(*   {True, True}   *)


Edit

Most simple way to get solutions, is just apply Exp  to both sides of equation.

neweq = Map[Exp[#] &, eq, {1}]

(*   E^x == E^y/(1 + y)   *)

Solve[neweq && x > 0, y, Reals]

(*   {{y -> ConditionalExpression[-1 - ProductLog[-E^(-1 - x)],
x > 0]}, {y ->
ConditionalExpression[-1 - ProductLog[-1, -E^(-1 - x)], x > 0]}}   *)