contours Extrait de ContourPlot en Mathematica

Question

I ont une fonction f (x, y) de deux variables, dont je besoin de connaître l'emplacement des courbes au cours de laquelle il passe par zéro. ContourPlot fait cela très efficace (qui est: il utilise des méthodes intelligentes multi-grilles, pas seulement un balayage à grain fin force brute), mais me donne juste un terrain. Je voudrais avoir un ensemble de valeurs {x, y} (avec une résolution spécifiée) ou peut-être une fonction interpoler qui me permet d'avoir accès à l'emplacement de ces contours.

ont pensé à extraire ce de la FullForm de ContourPlot, mais cela semble être un peu d'un hack. Une meilleure façon de le faire?

Answer 1

If you end up extracting points from ContourPlot, this is one easy way to do it:

points = Cases[
  Normal@ContourPlot[Sin[x] Sin[y] == 1/2, {x, -3, 3}, {y, -3, 3}],
  Line[pts_] -> pts,
  Infinity
]

Join @@ points (* if you don't want disjoint components to be separate *)

EDIT

It appears that ContourPlot does not produce very precise contours. They're of course meant for plotting and are good enough for that, but the points don't lie precisely on the contours:

In[78]:= Take[Join @@ points /. {x_, y_} -> Sin[x] Sin[y] - 1/2, 10]

Out[78]= {0.000163608, 0.0000781187, 0.000522698, 0.000516078, 
0.000282781, 0.000659909, 0.000626086, 0.0000917416, 0.000470424, 
0.0000545409}

We can try to come up with our own method to trace the contour, but it's a lot of trouble to do it in a general way. Here's a concept that works for smoothly varying functions that have smooth contours:

1. Start from some point (pt0), and find the intersection with the contour along the gradient of f.

2. Now we have a point on the contour. Move along the tangent of the contour by a fixed step (resolution), then repeat from step 1.

Here's a basic implementation that only works with functions that can be differentiated symbolically:

rot90[{x_, y_}] := {y, -x}

step[f_, pt : {x_, y_}, pt0 : {x0_, y0_}, resolution_] :=
 Module[
  {grad, grad0, t, contourPoint},
  grad = D[f, {pt}];
  grad0 = grad /. Thread[pt -> pt0];
  contourPoint = 
    grad0 t + pt0 /. First@FindRoot[f /. Thread[pt -> grad0 t + pt0], {t, 0}];
  Sow[contourPoint];
  grad = grad /. Thread[pt -> contourPoint];
  contourPoint + rot90[grad] resolution
 ]

result = Reap[
   NestList[step[Sin[x] Sin[y] - 1/2, {x, y}, #, .5] &, {1, 1}, 20]
];

ListPlot[{result[[1]], result[[-1, 1]]}, PlotStyle -> {Red, Black}, 
 Joined -> True, AspectRatio -> Automatic, PlotMarkers -> Automatic]

The red points are the "starting points", while the black points are the trace of the contour.

EDIT 2

Perhaps it's an easier and better solution to use a similar technique to make the points that we get from ContourPlot more precise. Start from the initial point, then move along the gradient until we intersect the contour.

Note that this implementation will also work with functions that can't be differentiated symbolically. Just define the function as f[x_?NumericQ, y_?NumericQ] := ... if this is the case.

f[x_, y_] := Sin[x] Sin[y] - 1/2

refine[f_, pt0 : {x_, y_}] := 
  Module[{grad, t}, 
    grad = N[{Derivative[1, 0][f][x, y], Derivative[0, 1][f][x, y]}]; 
    pt0 + grad*t /. FindRoot[f @@ (pt0 + grad*t), {t, 0}]
  ]

points = Join @@ Cases[
   Normal@ContourPlot[f[x, y] == 0, {x, -3, 3}, {y, -3, 3}],
   Line[pts_] -> pts,
   Infinity
   ]

refine[f, #] & /@ points

Answer 2

A slight variation for extracting points from ContourPlot (possibly due to David Park):

pts = Cases[
   ContourPlot[Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}], 
   x_GraphicsComplex :> First@x, Infinity];

or (as a list of {x,y} points)

ptsXY = Cases[
   Cases[ContourPlot[
     Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}], 
    x_GraphicsComplex :> First@x, Infinity], {x_, y_}, Infinity];

Edit

As discussed here, an article by Paul Abbott in the Mathematica Journal (Finding Roots in an Interval) gives the following two alternative methods for obtaining a list of {x,y} values from ContourPlot, including (!)

ContourPlot[...][[1, 1]]

For the above example

ptsXY2 = ContourPlot[
    Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}][[1, 1]];

and

ptsXY3 = Cases[
   Normal@ContourPlot[
     Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}], 
   Line[{x__}] :> x, Infinity];

where

ptsXY2 == ptsXY == ptsXY3