Plot of implicit function

Question

I have an equation F(t1,t2) = 0 (1). Please, could you tell me, how I can make a plot of function E = g(t1(t2),t2), where t1(t2) is the solution of equation (1) (I don't know its explicit form)

f[t1_, t2_] := -Sin[t1] - (-Cos[t1] + Cos[t2])/(t2 - t1); (* - this is Equation (1) *)
Eee[t1_, t2_] :=                                (* - this is the function, *)
 (-Sin[t2] - (-Cos[t1] + Cos[t2])/(t2 - t1))^2/2; (* which plot from t2 I want to make. *)

Answer 1

(I'm not very satisfied with the method of obtaining the result; I haven't found yet another way to accomplish the task.)

See also a closely related question.

Eee[t1_, t2_] := (-Sin[t2] - (-Cos[t1] + Cos[t2])/(t2 - t1))^2/2
f[t1_, t2_] := -Sin[t1] - (-Cos[t1] + Cos[t2])/(t2 - t1)

For some insight:

ContourPlot[{Eee[t1, t2]}, {t1, 0, 2 π}, {t2, 0.001, 2 π}, FrameLabel -> {"t1", "t2"}]

and

ContourPlot[{f[t1, t2]}, {t1, 0, 2 π}, {t2, 0.001, 2 π}, 
 FrameLabel -> {"t1", "t2"}, PlotLegends -> Automatic]

So

plot = ContourPlot[{f[t1, t2] == 0}, {t1, 0, 2 π}, {t2, 0.001, 2 π}, FrameLabel -> {"t1", "t2"}]

The line $y=x$ is asymptotic, as f[t,t]=1/0 and Limit[f[t1, t2], t1 -> t2] is 0.

One can extract the points directly from the plot:

points = Cases[
   Cases[Normal@plot, List[___], Infinity], {_?NumericQ, _?NumericQ}, 
   Infinity];

Verify:

ListPlot[points, AspectRatio -> 1, Frame -> True]

So one can generate a list of pairs {t2, Eee}:

Epoints = Eee[#[[1]], #[[2]]] & /@ points;
data = Transpose@{points[[All, 2]], Epoints};

ListPlot[data, Frame -> True, FrameLabel -> {"t2", "Eee"}]

---

Or without the line $y=x$ (which happens to correspond to Eee[t2]==0):

points2 = DeleteCases[points, {x_, y_} /; Abs[x - y] < 0.01];

ListPlot[points2, AspectRatio -> 1, Frame -> True]

Epoints2 = Eee[#[[1]], #[[2]]] & /@ points2;
data2 = Transpose@{points2[[All, 2]], Epoints2};

ListPlot[data2, Frame -> True, FrameLabel -> {"t2", "Eee"}]

---

One needs to examine a wider range of t1 and t2, as the contour plot of f[t1, t2]==0 seems to display different branches, which is reflected in the plot of Eee vs. t2:

plot = ContourPlot[{f[t1, t2] == 0}, {t1, 0, π}, {t2, -8 π, 8 π}, 
  FrameLabel -> {"t1", "t2"}, PlotPoints -> 100]

points = Cases[
   Cases[Normal@plot, List[___], Infinity], {_?NumericQ, _?NumericQ}, 
   Infinity];
points2 = DeleteCases[points, {x_, y_} /; Abs[x - y] < 0.05];

ListPlot[points2, AspectRatio -> 1, Frame -> True, 
  PlotRange -> {{-0.1, π}, All}]

Epoints2 = Eee[#[[1]], #[[2]]] & /@ points2;
data2 = Transpose@{points2[[All, 2]], Epoints2};

plot12 = ListPlot[Sort @ data2, Frame -> True, 
  FrameLabel -> {"t2", "Eee"}, Joined -> True]

The straight line comes from some spurious point from the contour plot.

---

Or:

plot = ContourPlot[{f[t1, t2] == 0}, {t1, Pi/2, 3/2 Pi}, {t2, -8 Pi, 
   8 Pi}, FrameLabel -> {"t1", "t2"}, PlotPoints -> 100]

Epoints = Eee[#[[1]], #[[2]]] & /@ points;
data = Transpose@{points[[All, 2]], Epoints};

ListPlot[data, Frame -> True, FrameLabel -> {"t2", "Eee"}]

Answer 2

The more or less exact solution t1B below is nice, but this numeric one is ten times as fast. The heuristic starting point for FindRoot was done by eye & manual iteration. This one is more practical, certainly for plotting. When used for the graph at the bottom of this answer, the plot takes less than a second to compute.

ClearAll[t1A];
t1A[t2_?NumericQ] := \[FormalT] /.    (* base solution *)
   FindRoot[
    f[\[FormalT], t2],
    {\[FormalT], 
     Pi/2 - ArcTan[t2/2] - (ArcTan[t2/2]/(8 + Abs[t2]) - 0.8 Sinc[t2^2/(1 + Abs[t2])])}];
t1A[t_?NumericQ, n_Integer] :=        (* translated to branch n -- see original answer *)
  t1A[t - π n] + π n;

t1B[3., 4] // AbsoluteTiming
t1A[3., 4] // AbsoluteTiming
(*
  {0.00724, 15.7072}
  {0.000531, 15.7072}
*)

---

Rahul's method from Define a function with variables linked implicitly:

ParametricPlot[{t2, Eee[t1, t2]}, {t1, 0, Pi}, {t2, -4 Pi, 4 Pi},
 MeshFunctions -> {-Cos[#3] + Cos[#4] + Sin[#3] (-#3 + #4) &}, 
 Mesh -> {{0}}, 
 MeshStyle -> 
  Directive[RGBColor[0.368417, 0.506779, 0.709798], 
   AbsoluteThickness[1.6`], Opacity[1]], PlotPoints -> 150, 
 MaxRecursion -> 1, PlotStyle -> None, BoundaryStyle -> None, 
 Exclusions -> None, AspectRatio -> 0.6]

---

Original answer

A solution to the OP's equation f[t1, t2] == 0, which happens to be equivalent to the equation in How to find roots of equation, can be found by adapting my answer there. My answer gave a parameterization of part of the solution set. This may be adapted to parametrize a complete branch, and by periodicity, this can be extended to all branches. InverseFunction may be used to construct t1 as a function of t2 from the parametrization.

res = θ -> -(t/2) + ArcTan[2 - t Cot[t/2], t];  (* from my answer *)

(* check res *)
FullSimplify[f[θ + t, θ] /. res]
(*  0  *)

ClearAll[param];
(* parametrizes {t2, t1} *)
param[t_, n_] = Pi*n - Pi*Quotient[t, 2 Pi] + {θ, θ + t} /. res;

Constructing a relatively efficient t1 requires that the value of InverseFunction[] be computed only once per evaluation. (The code can be generated from param[#1, #2].) The argument n determines which branch is computed.

ClearAll[t1B];
t1B[t_?NumericQ, n_Integer] := 
  ArcTan[2 - Cot[#/2] #, #] - π Quotient[#, 2 π] + #/2 + π*n &[
   InverseFunction[ArcTan[2 - Cot[#/2] #, #] - π Quotient[#, 2 π] - #/2 + π*n &][t]];

We need to check that InverseFunction computes computes the "right" branch, or at least stays on one branch:

cplot = ContourPlot[        (* a plot of several branches *)
   f[t1, t2] == 0,
   {t2, -8 π, 8 π}, {t1, 0.00001 - 8 π, 8 π},
   FrameLabel -> Automatic, 
   GridLines -> {Range[Pi/2 - 8 Pi, 8 Pi, Pi], Range[Pi/2 - 8 Pi, 8 Pi, Pi]}];

Show[                       (* comparison *)
 cplot,
 Plot[t1B[t2, 0], {t2, -8. Pi, 8. Pi}, PlotStyle -> Red]
 ]

Here are a bunch of them, which completely mask the contour plot.

Show[
 cplot,
 ListLinePlot[Table[{t2, t1B[t2, n]}, {n, -9, 7}, {t2, -8 Pi, 8 Pi, 0.2}]]
 ]

On to the OP's problem of plotting Eee[t1B[t2, 0], t2]. It's somewhat faster to not evaluate Eee until t1B[t2, 0] has been computed, because InverseFunction takes a small but appreciable amount of time to evaluate. It's much faster to evaluate it on approximate reals, since with exact inputs to InverseFunction[], Mathematica can sometimes take a very long time to determine the value.

Block[{y},
 y[t_?NumericQ] := Eee[t1B[N@t, 0], t];
 Plot[y[t2], {t2, -20, 20}]
 ]

(The plot takes almost 13 seconds on my laptop.)