9
$\begingroup$

This is how the curve looks like:

img = ContourPlot[
                  1/x + 3/4 (((y - 1/Sqrt[3])/x)^2 + 1) Exp[
                  ArcTan[(y - 1/Sqrt[3])/x] - Pi/6] == 0,
           {x, -3, 1}, 
           {y, -(1/5), 4},
       PlotPoints -> 70 ]

enter image description here

@xzczd comes up with a (kinka hacky!) solution, which extract the coordinates forming that curve:

Total[EuclideanDistance @@@ 
  Partition[First@Cases[Normal@img, 
                        Line[a_] :> a, Infinity], 2, 1]
      ]

  (* 9.85614 *)

Is there any better way to do this?

$\endgroup$
5
  • 2
    $\begingroup$ "kinka hacky" is nice as it is more generic? $\endgroup$
    – chris
    Jul 27, 2013 at 17:07
  • $\begingroup$ @chris maybe :D $\endgroup$
    – mmjang
    Jul 27, 2013 at 17:17
  • 2
    $\begingroup$ @chris The fatal defect of this method is, it can't control the Precision. The result is always MachinePrecision, and the influence from PlotPoints is also big: PlotPoints -> 70 gives 9.92371, -> 200 gives 9.85703, -> 300 gives 9.84398, -> 400 gives 9.84211 $\endgroup$
    – xzczd
    Jul 28, 2013 at 6:10
  • $\begingroup$ @xzczd sure; on the other hand its not always trivial to find a parametrization as you did below?? $\endgroup$
    – chris
    Jul 28, 2013 at 7:57
  • $\begingroup$ @chris So I'm still looking forward to a general numeric solution though b.gatessucks had solved it analytically :) $\endgroup$
    – xzczd
    Jul 28, 2013 at 9:44

3 Answers 3

17
$\begingroup$

You can get a parametric representation for your curve :

eqn = 1/x +  3/4 (((y - 1/Sqrt[3])/x)^2 + 1) Exp[ArcTan[(y - 1/Sqrt[3])/x] - Pi/6] ;

aux = First@Solve[(eqn /. {y -> 1/Sqrt[3] + t x}) == 0, x]
(* {x -> -((4 E^(\[Pi]/6 - ArcTan[t]))/(3 (1 + t^2)))} *)

solx = aux[[1, 2]]
(* -((4 E^(\[Pi]/6 - ArcTan[t]))/(3 (1 + t^2))) *)

soly = 1/Sqrt[3] + t x /. aux
(* 1/Sqrt[3] - (4 E^(\[Pi]/6 - ArcTan[t]) t)/(3 (1 + t^2)) *)

Plot[{solx, soly}, {t, -50, 50}, PlotRange -> All]

sol

It looks like you get the curve with t in [-500,500] - you might need to improve on the interval.

ParametricPlot[{solx, soly}, {t, -500, 500}, PlotRange -> All, PlotPoints -> 200]

plot

The last step is just the usual definition of arclength (thanks to @MichaelE2):

NIntegrate[ Sqrt[Simplify[(D[solx, t])^2 + (D[soly, t])^2 ]], {t, -Infinity, Infinity}]
(* 9.83926 *)
$\endgroup$
3
  • 4
    $\begingroup$ Very nice, but I think the limits should be {t, -Infinity, Infinity}, at least on the integral. $\endgroup$
    – Michael E2
    Jul 27, 2013 at 18:35
  • 1
    $\begingroup$ Oh, parametric representation… I should have not given up the effort for getting analytic solution so easily… $\endgroup$
    – xzczd
    Jul 28, 2013 at 6:24
  • $\begingroup$ @MichaelE2 You're right, thanks. $\endgroup$ Jul 28, 2013 at 7:48
4
$\begingroup$

In version 10, you can also do this:

eqn = 1/x + 3/4 (((y - 1/Sqrt[3])/x)^2 + 1) Exp[ArcTan[(y - 1/Sqrt[3])/x] - Pi/6] == 0;
region = ImplicitRegion[eqn, {{x, -3, 1}, {y, -1/5, 4}}];
ArcLength@DiscretizeRegion[region, AccuracyGoal -> 6]
(* 9.83926 *)

But you have to set the AccuracyGoal yourself, and I'm not sure it gives any guarantees on the accuracy of the arc length itself. Sadly applying ArcLength directly to region fails with "Unable to compute the length of region ImplicitRegion[...]."

$\endgroup$
2
$\begingroup$

In the case that one cannot solve the equation for a symbolic parametric representation, then NDSolve can be used to do so numerically. And while we're at it, we may as well integrate the arclength. In the code below, we compute an arc length parametrization, so the parametrization returns to it's starting point when the parameter s equals the total arc length of the loop.

ClearAll[f, x, y];
eqn = 1/x +  3/4 (((y - 1/Sqrt[3])/x)^2 + 1) Exp[ArcTan[(y - 1/Sqrt[3])/x] - Pi/6] == 0;
cplot = ContourPlot[Evaluate @ eqn, {x, -3, 1}, {y, -(1/5), 4}];
f[x_, y_] = eqn /. Equal -> Subtract // Together // Numerator // Simplify;
grad[x_, y_] = D[f[x, y], {{x, y}}];
unitTangent[x_, y_] = #/Sqrt[#.#] &@Cross@grad[x, y];
p0 = NestWhile[                       (* Newton's method to find starting point *)
   With[{g = grad @@ #},                 (* use gradient for derivative *)
     # - (f @@ #) g / g.g                (* Newton's method step *)
     ] &,
   cplot[[1, 1, 1]],                     (* start at a point on the contour plot *)
   Abs[#1 - #2]/Norm[#1] > 1*^-15 &,     (* stopping criterion *)
   2,
   20                                    (* no more than 20 iterations *)
   ];

sol = First@NDSolve[{
     {x'[s], y'[s]} == unitTangent[x[s], y[s]], {x[0], y[0]} == p0,
     WhenEvent[x[s] > p0[[1]], "StopIntegration"]},
    {x, y}, {s, 0, Infinity}];

x["Coordinates"] /. sol // Last // Last
% - NIntegrate[
  Sqrt[Simplify[(D[solx, t])^2 + (D[soly, t])^2]], {t, -Infinity, 
   Infinity}]
(*
  9.83926         - arc length
  8.24538*10^-7   - error (compared to b.gatessucks's result)
*)

Remarks: (1) One might need extra conditions for the stopping condition WhenEvent[x[s] > p0[[1]]... in the case of more complicated curve. (2) The error when starting at p0 = cplot[[1, 1, 1]] is almost 0.001, so it is probably worth improving it in most cases.

$\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.