# Analytical solution of a transcendental equation

I have the next equation

eqh = u == Istim + 100 gE (2/(1 + Exp[(-2 (u - u0))/kv]) - 1) (xD - 0.5);
Istim = 126; gE = 5;kv=20;  u0=150;


and I want to solve it for the variable u in the region of u>u0. Technically, it could be found with some conditions, as it's seen from the plot

xDimp = Solve[eqh, xD];
Plot[Evaluate[xD /. xDimp] , {u, -100, 400}];


But neither Solve nor Reduce don't calculate the explicit expression for u.

Solve[eqh, u]
Solve[eqh && u > u0 , u, Reals]
Reduce[eqh && u > u0 , u, Reals]


arguing it with the next message:

Reduce::inex: Reduce was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Reduce require exact input, providing Reduce with an exact version of the system may help.

• Do you have any reason to believe that the solution is even expressible in terms of elementary functions? IMO questions about analytic solution transcendental equations should always provide some evidence that the solution can be expressed as an explicit formula. – Szabolcs Aug 31 '17 at 9:11
• I'll note that if you give xD an explicit value, Reduce can give you an exact solution in the sense that it can guarantee that all solutions are found and it can compute them to arbitrary precision. It cannot express them using a formula though. – Szabolcs Aug 31 '17 at 9:11
• It may sounds silly, but I got the explicit solution yesterday running the same code. I can't specify the exact value of xD because the problem I'm describing here is an auxiliary one for examining of vector field of the differential equation specified here mathematica.stackexchange.com/questions/154797/… – Artem Zefirov Aug 31 '17 at 9:38
• I don't believe there is an explicit solution to your equation. Even something as simple as $x = e^x$ doesn't have one. That's why Mathematica will return a solution in terms of ProductLog (Lambert W function). ProductLog is just a workaround for not being able to express the solution of an equation otherwise: let's give the solution a name then, prove various properties about it, and come up with a way to calculate it numerically. This is what you can do with your equation as well. – Szabolcs Aug 31 '17 at 9:42
• I see. Thank you for clarifications. – Artem Zefirov Aug 31 '17 at 10:32

This is not really an answer, but a comment. I write it in the window for the answers since it requires illustrations.

In addition to what @Szabolcs wrote in his comment, not only that this equation (a nonlinear Fermi-like distribution) has no exact solution, but it also exhibits a bifurcation. To make sure let us plot the left- and the right-hand parts of your equations for different values of the parameter xD. That can be done with Manipulate:

    Manipulate[
Plot[{u,
126 + 500 (-1 + 2/(1 + E^((150 - u)/10))) (-1/2 + xD)}, {u,-100,400}],
{xD, -1, 1, Appearance -> "Labeled"}]


As you see I have chosen the same interval for u as you did, and the interval (-1,1) for xD. At xD=-1 it is shown below:

The blue line is the lhp and the yellow one is the rhp of the equation. As we see there is a single solution at this value of the parameter.

At xD=approx=0.616one observs a bifurcation point:

In which a new solution shows up, which then splits into two ones at higher xDvalues:

In the situations involving a bifurcation an analytical solution can be obtained in exceptional cases, and this is not such an exceptional case.

However, the bifurcation point itself, as well as the solution in its close vicinity can be obtained without great difficulties.

Also, a numeric solution of each of these three solutions can be easily tabulated as a function of xD.

However, first of all, it should be decided, which one of these solutions you are interested in.

Yes, it is easy to find the bifurcation point. Its condition is as follows. One needs to calculate a derivative of the equation. In your case equation is:

  eqh1 = u == 126 + 500 (-1 + 2/(1 + E^((150 - u)/10))) (-(1/2) + xD);


And here is the derivative:

eqh2 = D[eqh1[[2]], u] == 1 // Simplify

(*  (50 E^(15 + u/10) (-1 + 2 xD))/(E^15 + E^(u/10))^2 == 1  *)


In the bifurcation point they must be satisfied simultaneously:

sl1 = Solve[{eqh1, eqh2}, {u, xD}, Reals] // N

(* {{u -> 172., xD -> 0.615}}  *)


which gives you the bifurcation point.

Now let us express xD in terms of u:

sl2 = Solve[eqh1, xD][[1, 1, 2]]


and make a table with the elements {xD, u}:

 lst = Table[{sl2 // N, u}, {u, 151, 170, 0.5}]


They are selected such that xD>0.651 and 151<u<170 which garantees that the middle solution is selected.

This solution can be built as a ListPlot. Another way to build it as a ParametricPlot:

    Show[{
ParametricPlot[{sl2, u}, {u, 151, 170}, AspectRatio -> 0.7,
AxesLabel -> {Style["xD", Italic, 16], Style["u", Italic, 16]}],
ListPlot[lst, PlotStyle -> Red]
}]


yielding the following plot:

Thus, the list above represents a tabulated solution that you need. One can further find some simple function that accurately fitts to this solution.

Have fun!

• Much clear now. I'm interested in the middle solution, near u=150. In more common problem parameter xD decreases from 1 to some value, No I see that this value is a bifurcation point xD=0.616. So the solution when xD in [0.616, 1] is what I need. Do you mean, that it's possible to find the bifurcation point value numerically and then linearize the equation in its vicinity? – Artem Zefirov Aug 31 '17 at 12:21
• @Artem Zefirov Have a look at the edits. – Alexei Boulbitch Aug 31 '17 at 14:24