Transcendental equation inside another equation when 3D plotting

Question

I have been trying to solve this plotting problem myself, and I have received a lot of very useful advice from here on Stack Exchange - but I am now at a complete loss as to how to achieve the plot I want. I am trying to do a 3D plot of the equation z = -t*Log[2*Cosh[(2y + x)/t]] where y is defined by the transcendental equation y = Tanh[(2y + x)/t]. Specifically the relation between the two equations is that y = - ∂z/∂x.

In particular, y = Tanh[(2y + x)/t] contains a supercritical bifurcation - hence, the equation has three solutions for t < 2, and I have to pick one solution depending on the sign of x. However, I have figured out that a simple way to get around this problem, is to use that y = Tanh[(2y + x)/t] is symmetric around x = 0, and then simply just doing the 3D plot in two halves up to x = +/- 0.001.

Now, I have succeeded in plotting y = Tanh[(2y + x)/t] on my own with no problems - however, when I try to plot z = -t*Log[2*Cosh[(2y + x)/t]] as best as I can, I get a host of error messages and a beautiful surface plot that looks very close to correct - yet, but I can tell from cuts through the surface that the plot simply isn't correct, because the slope of z just isn't fulfilling y = - ∂z/∂x.

I am very much a noob at Mathematica, so I have done the best I could by bodging and tinkering - but I am completely stomped as to whether I am missing something mathematically or whether it is Mathematica problem, because I am not plugging y correctly into z. Could anybody possibly show me the correct way of doing this 3D plot of z, because I think I just may be in over my head? <:)

Addition:

I thought I'd better post the code I have used:

g[x_, t_] := y /. FindRoot[y == Tanh[(2 y + x)/t], {y, -5}];
p1 = ContourPlot3D[ z == - t * Log[2*Cosh[(2 g[x, t] + x)/t]] , 
{t, 0.001, 3}, {x, -1.5, -0.001}, {z, -4.0, -1.0}, 
MeshFunctions -> {Function[{t, x, z}, 
Evaluate[  D[z + t * Log[2*Cosh[(2 g[x, t] + x)/t]], z]]    ]}, 
Mesh -> {{0}}, MeshShading -> Lighter /@ ColorData[1] /@ {1, 2}, 
AxesLabel -> {t, x, z}, Lighting -> "Neutral"];

k[x_, t_] := y /. FindRoot[y == Tanh[(2 y + x)/t], {y, 5}];
p2 = ContourPlot3D[ z == - t * Log[2*Cosh[(2 k[x, t] + x)/t]] , 
{t, 0.001, 3}, {x, 0.001, 1.5}, {z, -4.0, -1.0}, 
MeshFunctions -> {Function[{t, x, z}, 
Evaluate[  D[z + t * Log[2*Cosh[(2 k[x, t] + x)/t]], z]]    ]}, 
Mesh -> {{0}}, MeshShading -> Lighter /@ ColorData[1] /@ {1, 2}, 
AxesLabel -> {t, x, z}, Lighting -> "Neutral"];

Show[p1, p2, PlotRange -> {{0.0, 3.0}, {-1.5, 1.5}, {-4.0, -1.0}}, 
ViewPoint -> {0.9*Pi, -0.45*Pi, 1.05*Pi} ]

As I wrote above, this is the product of quite a bit of tinkering and bodging because I am very new to this, so I am not completely sure what everything is doing - this just seems to do the trick. The main errors I am getting are: FindRoot::nlnum: "The function value {-5.-1.\ Tanh[(-10.+x)/t]} is not a list of numbers with dimensions {1} at {y} = {-5.}." and ReplaceAll::reps: "{FindRoot[y==Tanh[(2\y+x)/t],{y,-5}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing." - but I am not sure what to make of them.

Answer 1

I think the technique I used in this answer could be applied to your question as well.

paraRegion = ContourPlot3D[
  y == Tanh[(2 y + x)/t], (* the constraint on x, y, t *)
  {x, -2, 2},
  {y, -1, 1},
  {t, 0, 2},
  AxesLabel -> (Style[#, 20, Bold] & /@ {x, y, t})]

Cases[paraRegion,
   GraphicsComplex[pts_, others_,
     opts1___, VertexNormals -> vn_, opts2___] :>
    GraphicsComplex[
     Function[{x, y, t},
       {x, y, -t*Log[2*Cosh[(2 y + x)/t]]} (* the surface in xyz space *)
       ] @@@ pts,
     others, opts1, opts2], ∞][[1]] //
 Graphics3D[#, Axes -> True, 
   AxesLabel -> (Style[#, 20, Bold] & /@ {x, y, z}), 
   PlotRange -> All] &

Edit:

To illustrate the surface indeed fulfills $y=-\partial_x z$ (and actually a ruled surface), we can add options as following to the ContourPlot3D for paraRegion to emphasize the feature:

ContourStyle -> None,
MeshFunctions -> {#1 &, #2 &},
Mesh -> {10, 20},
MeshStyle -> {Gray, Directive[Darker[Blue], Thick]}

which give

Edit 2:

According to OP's request in the comment, to plot the surface in $txz$ space, there is no essential difference. We can just change the parametric equations in $xyz$ space in the line I marked with (* the surface in xyz space *) to the equations in $txz$ space, that is

{t, x, -t*Log[2*Cosh[(2 y + x)/t]]}

The key point is, once we have a self-contained feasible set $\left\{\boldsymbol{p}\mid\boldsymbol{p}=(t,x,y,z)\land F(t,x,y,z)=0\right\}$ (e.g. the points on the surface of paraRegion along the relation z == f[x, y, t]), it's up to us in which space we'd like to represent it - just strip out the extra dimension and it's done.

Answer 2

One can adapt (simplify somewhat) this answer to this case.

It will be convenient to compute a NormalsFunction and to parametrize the plot in terms of y as defined by zero set of F = t ArcTanh[y] - (2*y + x) (which is computed by g and k in the OP from the equivalent equation y == Tanh[(2 y + x)/t]).

Surface normal function

From F == 0, we can compute $\partial y/\partial t$ (dydt) and $\partial y/\partial x$ (dydx). These can be plugged into the gradient (normal) of the surface z == -t*Log[2*Cosh[(2 y[t, x] + x)/t]], the normal can be multiplied by a scalar, both to clear the denominator and to get it to point "up" (with respect to z). Then we can apply a couple of identities and simplify.

F = t ArcTanh[y] - (2*y + x);
dydt = -D[F, t]/D[F, y];
dydx = -D[F, x]/D[F, y];

2 t (-2 + t + 2 y^2)*
     D[z + t*Log[2*Cosh[(2 y[t, x] + x)/t]], {{t, x, z}}] /.
        {Derivative[1, 0][y][t, x] -> dydt, Derivative[0, 1][y][t, x] -> dydx, 
              Cosh[u_] :> 1/Sqrt[1 - Tanh[u]^2]} /. y[t, x] -> y /. 
          Tanh[(x + 2 y)/t] -> y // PowerExpand // Expand // FullSimplify

(* {-2 y (x + 2 y) (-2 + t + 2 y^2) +
       t (4 y (-1 + y^2) ArcTanh[y] + t Log[4] + (-1 + y^2) Log[16] -
          (-2 + t + 2 y^2) Log[1 - y^2]), 
    2 t^2 y,
    2 t (-2 + t + 2 y^2)} *)

The formula from this output may be used to define a surface normal function for plotting:

normalsFn = 
  Compile[{t, x, y},
    {-2 y (x + 2 y) (-2 + t + 2 y^2) +
         t (4 y (-1 + y^2) ArcTanh[y] + t Log[4] + (-1 + y^2) Log[16] -
            (-2 + t + 2 y^2) Log[1 - y^2]), 
     2 t^2 y,
     2 t (-2 + t + 2 y^2)} ];

Parametrization

There is a bifurcation where the y partial of F is zero. We can solve F == 0 for x in terms of y and t; it happens that x is a linear function of t, which makes the problem quite tractable. We can use this solution to find a parametrization of t, x, such that the bifurcation happens when u == 0, and the end points u == -1 and u == 1 correspond to t == 0 and t = 3 respectively.

The functions tP, xP, zP parametrize the t, x, z coordinates of the surface in terms of y and a parameter u. We use some utility functions in defining them:

• xOFty: returns x for a given t, y;
• tOFxy: returns t for a given x, y;
• xBrOFy: returns the x coordinate on the bifurcation curve for a given y.

We give different definitions of tP and xP depending on the side of the bifurcation curve. (See also the linked answer.) We define the parametrization in terms of the domain {{tMin, tMax}, {xMin, xMax}}, for the purposes of rescaling, which domain is taken from the OP.

xOFty[t_, y_] := Evaluate[x /. Simplify@First@Solve[F == 0, x]];
tOFxy[x_, y_] := Evaluate[t /. Simplify@First@Solve[F == 0, t]];
xBrOFy[y_] := Evaluate[x /. Simplify@First@Solve[F == 0 /. First@Solve[D[F, y] == 0, t], x]];

Clear[xP, tP, zP];
xP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; u >= 0 && y >= 0 := 
  Rescale[u, {0, 1}, {xBrOFy[y], Min[x2, xOFty[t2, y]]}];
xP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; u >= 0 && y < 0 := 
  Rescale[u, {0, 1}, {xBrOFy[y], Max[x1, xOFty[t2, y]]}];
xP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; u < 0 := 
  Rescale[u, {-1, 0}, {-2 y, xBrOFy[y]}];

tP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; u >= 0 && y == 0 := 
  Rescale[u, {0, 1}, {2, t2}];
tP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; u < 0 && y == 0 := 
  Rescale[u, {-1, 0}, {t1, 2}];
tP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; -1 < y < 1 := 
  tOFxy[xP[u, y, {{t1, t2}, {x1, x2}}], y];

zP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] := -tP[u, y, {{t1, t2}, {x1, x2}}] *
   Log[2*Cosh[(2 y + xP[u, y, {{t1, t2}, {x1, x2}}]) / tP[u, y, {{t1, t2}, {x1, x2}}]]];

{tMin, tMax} = {0, 3}; {xMin, xMax} = {-2, 2};
domain = {{tMin, tMax}, {xMin, xMax}};
param[u_?NumericQ, y_?NumericQ] /; -1 < y < 1 :=
   {tP[u, y, domain], xP[u, y, domain], zP[u, y, domain]};

Plots

As before (in my other answer), we need more plot points near y == -1 and y == 1, so we use Sin[(Pi/2) v] to interpolate between -1, +1. (We use it for the sake of convenience.)

The whole surface. To get good result near y == -1 and y == 1, we increase WorkingPrecision a little.

ParametricPlot3D[param[u, Sin[(Pi/2) v]],
 {u, -1, 1}, {v, -1, 1},
 MeshFunctions -> {#3 &}, NormalsFunction -> (normalsFn[#1, #2, Sin[(Pi/2) #5]] &), 
 AxesLabel -> {t, x, z}, ImageSize -> 300, WorkingPrecision -> 20]

The part indicated in the OP. We can restrict the plot with a slight redefinition of xP (keeping the other definitions) and limiting the domain of the parameter u; however, over this domain, we do not need to increase WorkingPrecision.

xP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; u >= 0 && y >= 0 := 
  Rescale[u, {0, 1}, {(*xBrOFy[y]*)0, Min[x2, xOFty[t2, y]]}];
xP[u_, y_, {{t1_, t2_}, {x1_, x2_}}] /; u >= 0 && y < 0 := 
  Rescale[u, {0, 1}, {(*xBrOFy[y]*)0, Max[x1, xOFty[t2, y]]}];

ParametricPlot3D[param[u, Sin[(Pi/2) v]],
 {u, 0, 1}, {v, -1, 1},
 MeshFunctions -> {#3 &}, NormalsFunction -> (normalsFn[#1, #2, Sin[(Pi/2) #5]] &), 
 AxesLabel -> {t, x, z}, ImageSize -> 300]