Extraer contornos de Contourplot en Mathematica

Question

Tengo una función f (x, y) de dos variables, de las cuales necesito saber la ubicación de las curvas en la que cruza cero. Contourplot lo hace de manera muy eficiente (es decir: utiliza métodos inteligentes de griegas múltiples, no solo un escaneo de grano fino de fuerza bruta) sino que solo me da una parcela. Me gustaría tener un conjunto de valores {x, y} (con alguna resolución especificada) o tal vez alguna función de interpolación que me permita obtener acceso a la ubicación de estos contornos.

He pensado en extraer esto de la forma completa de Contourplot, pero esto parece ser un poco hack. ¿Alguna mejor manera de hacer esto?

Answer 1

Si terminas extrayendo puntos de ContourPlot, esta es una manera fácil de hacerlo:

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 *)

EDITAR

Parece que ContourPlot no produce contornos muy precisos. Por supuesto, están destinados a trazar y son lo suficientemente buenos para eso, pero los puntos no se encuentran con precisión en los contornos:

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}

Podemos tratar de crear nuestro propio método para rastrear el contorno, pero es mucho problema hacerlo de manera general. Aquí hay un concepto que funciona para funciones que varían sin problemas que tienen contornos suaves:

1. Comience desde algún punto (pt0), y encuentre la intersección con el contorno a lo largo del gradiente de f.

2. Ahora tenemos un punto en el contorno. Moverse a lo largo de la tangente del contorno por un paso fijo (resolution), luego repita desde el paso 1.

Aquí hay una implementación básica que solo funciona con funciones que se pueden diferenciar simbólicamente:

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]

Los puntos rojos son los "puntos de partida", mientras que los puntos negros son el rastro del contorno.

Edición 2

Quizás sea una solución más fácil y mejor usar una técnica similar para hacer los puntos que obtenemos ContourPlot más preciso. Comience desde el punto inicial, luego mueva a lo largo del gradiente hasta que se cruzamos el contorno.

Tenga en cuenta que esta implementación también funcionará con funciones que no se pueden diferenciar simbólicamente. Simplemente defina la función como f[x_?NumericQ, y_?NumericQ] := ... si este es el caso.

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

Una ligera variación para extraer puntos de ContourPlot (posiblemente debido a David Park):

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

o (como una lista de puntos {x, y})

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

Editar

Como se discutio aquí, un artículo por Paul Abbott en el Mathematica Journal (Encontrar raíces en un intervalo) proporciona los siguientes dos métodos alternativos para obtener una lista de valores {x, y} de contornplot, incluyendo (!)

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

Para el ejemplo anterior

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

y

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

dónde

ptsXY2 == ptsXY == ptsXY3