5
$\begingroup$

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. *)
$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Oct 9 '16 at 15:18
  • $\begingroup$ See How can I plot implicit equations? $\endgroup$ – user9660 Oct 9 '16 at 15:22
  • 1
    $\begingroup$ @Louis This is different because the OP wants, I think, Plot[e[t1[t2], t2], {t2, a, b}], where t1[t2] is defined by f[t1[t2], t2] == 0. This seems closer to a duplicate: Define a function with variables linked implicitly $\endgroup$ – Michael E2 Oct 9 '16 at 15:30
  • $\begingroup$ @MichaelE2, yes, agree. $\endgroup$ – user9660 Oct 9 '16 at 15:33
  • $\begingroup$ The equation f[t1, t2] == 0 is discussed here. $\endgroup$ – Michael E2 Oct 9 '16 at 19:32
2
$\begingroup$

(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"}]

enter image description here

and

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

enter image description here

So

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

enter image description here

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]

enter image description here

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"}]

enter image description here


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]

enter image description here

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

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

enter image description here


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]

enter image description here

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}]

enter image description here

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

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

enter image description here

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]

enter image description here

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

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

enter image description here

| improve this answer | |
$\endgroup$
2
$\begingroup$

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]

Mathematica graphics


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]
 ]

Mathematica graphics

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}]]
 ]

Mathematica graphics

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}]
 ]

Mathematica graphics

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

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.