1
$\begingroup$

I want to get equal value contour of a function. But I want the curve to be precise.

A naive way would be using ContourPlot, then extract the contour data. To do this, we need a utility to extract contours like below

ClearAll[extractContourAssoc];
extractContourAssoc[contourPlot_] := Module[{},
   tooltipData = 
    Cases[Normal@contourPlot, Tooltip[x__] -> {x}, Infinity];
   contourValueList = tooltipData[[;; , -1]];
   contourDataList = 
    Table[Cases[i, Line[pts__] -> pts, Infinity], {i, tooltipData}];
   AssociationThread[contourValueList -> contourDataList]];

Then for example

func[x_, y_] := x^2 + y^2;
cplot = ContourPlot[func[x, y] == 1, {x, 0, 3}, {y, 0, 3}];

to extract the contour in cplot , we use

pts = First@First@extractContourAssoc@cplot;

However, the precision is limited, for example, it we substitute back the coordinate of curves.

func @@ # - 1 & /@ pts
(*{-0.0025047, -0.00131712, -0.00137604, -0.00225737, -0.00143394, \
-0.00116087, -0.00113379, -0.00280045, -0.00107967, -0.00161928, \
-0.00143129, -0.000948526, -0.00279327, -0.00140771, -0.000650118, \
-0.000687417, -0.000736281, -0.000528992, -0.00106009, -0.000248064, \
-0.00138826, -0.00156147, -0.000895443, -0.00156147, -0.000757219, \
-0.00138826, -0.000248064, -0.00106009, -0.00126354, -0.000736281, \
-0.00190364, -0.00140771, -0.00034149, -0.000606891, -0.000948526, \
-0.00143129, -0.00161928, -0.00107967, -0.00280045, -0.000860253, \
-0.00113379, -0.00116087, -0.00143394, -0.00225737, -0.00137604, \
-0.00131712, -0.000643397}*)

If the curve is precise enough, then the numbers should all be zero.

Is there a way to get more precise contours? Though there is MaxRecursion option, but it seems that there is no precision control in ContourPlot.


update

Sorry I found I made a mistake about extractContourAssoc. Line[pts__] -> pts should be Line[pts__] :> pts

ClearAll[extractContourAssoc];
extractContourAssoc[contourPlot_] := Module[{},
   tooltipData = 
    Cases[Normal@contourPlot, Tooltip[x__] -> {x}, Infinity];
   contourValueList = tooltipData[[;; , -1]];
   contourDataList = 
    Table[Cases[i, Line[pts__] :> pts, Infinity], {i, tooltipData}];
   AssociationThread[contourValueList -> contourDataList]];

And MaxRecursion actually works. For example, Set MaxRecursion->10 will have precision roughly about 5x10^-6. So the problem is solved.

$\endgroup$
  • $\begingroup$ you can take your contourplot results and use as initial points for FindRoot $\endgroup$ – george2079 Mar 4 '18 at 13:53
1
$\begingroup$

Why not use an NDSolve based approach instead? For instance:

eqn = x[t]^2 + y[t]^2 == 1;
{xsol, ysol, diff} = Quiet @ NDSolveValue[
    {D[eqn, t], x'[t]^2 + y'[t]^2 == 1, x[0]==1, y[0]==0},
    {x, y, x[t]^2+y[t]^2-1}, 
    {t, 0, 2 Pi}
];

ParametricPlot[{xsol[t], ysol[t]}, {t, 0, 2Pi}]

enter image description here

I included the difference you were interested in as an output of the NDSolveValue call:

Plot[diff, {t, 0, 2Pi}]

enter image description here

Finally, if you extract the Line object from the parametric plot, you will see that its errors are also less than $10^{-5}$ as well.

$\endgroup$
  • $\begingroup$ Hi, @Carl Woll. It seems that your method needs to know at least one precise point at first, am I right? $\endgroup$ – matheorem Feb 2 '18 at 9:54
1
$\begingroup$

A typical way is to polish with FindRoot. Use the contour for one equation and the normal line for the second.

extractContourAssoc::noimpr = "No improvement possible in contour tagged ``";
extractContourAssoc[contourPlot_, {f_, x_, y_}] := 
  Module[{grad = D[f, {{x, y}}], contour},
   tooltipData = Cases[Normal@contourPlot, Tooltip[t__] :> {t}, Infinity];
   contourValueList = tooltipData[[;; , -1]];
   contourDataList = Table[
     Cases[
      i,
      Line[pts__] :> Block[{},
        (* check the Tooltip tag to see how to use it *)
        If[Head[i[[-1]]] === Equal,   (* tag is an equation with held variables *)
         contour = i[[-1]] /. Equal -> Subtract // ReleaseHold,
         If[NumericQ[i[[-1]]],        (* tag is numerical value *)
          contour = f - i[[-1]],
          Message[extractContourAssoc::noimpr, i[[-1]]];  
          contour = $Failed
          ]];
        (* improve the position if possible *)
        If[contour =!= $Failed,
         {x, y} /. 
            FindRoot[{contour, Cross[grad].({x, y} - #)}, Thread@{{x, y}, #}] & /@ pts,
         pts]
        ],
      Infinity],
     {i, tooltipData}];
   AssociationThread[contourValueList -> contourDataList]];

Example:

func[x_, y_] := x^2 + y^2;
cplot = ContourPlot[func[x, y] == 1, {x, 0, 3}, {y, 0, 3}];

The one-argument form of extractContourAssoc is the OP's function:

pts = First@First@extractContourAssoc@cplot;  (* OP's extractContourAssoc *)
func @@ # - 1 & /@ pts // RealExponent // Histogram

Mathematica graphics

The two-argument form is the one above that attempts to improve the position of the points with FindRoot:

pts = First@First@extractContourAssoc[cplot, {func[x, x], x, y}]; (* improved *)
func @@ # - 1 & /@ pts // RealExponent // Histogram[#, {1}] &

Mathematica graphics

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