Extraiga 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 a las que cruza cero. ContourPlot lo hace de manera muy eficiente (es decir, utiliza métodos inteligentes de múltiples cuadrículas, no solo un escaneo de grano fino de fuerza bruta), pero solo me da un argumento. Me gustaría tener un conjunto de valores {x, y} (con alguna resolución especificada) o quizás alguna función de interpolación que me permita acceder a la ubicación de estos contornos.

Hemos pensado en extraer esto del FullForm de ContourPlot pero parece ser un truco. ¿Alguna mejor manera de hacer esto?

Answer 1

Si al final de la extracción de 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. Son, por supuesto, significaba para el trazado y son lo suficientemente buenos para eso, pero los puntos no se encuentran precisamente 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 llegar a nuestro propio método para trazar el contorno, pero es un montón de problemas hacerlo de una manera general. Aquí hay un concepto que funciona para problemas diferentes funciones que tienen contornos lisos:

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

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

He aquí una implementación básica que sólo funciona con las 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 la traza del contorno.

EDIT 2

Tal vez sea una solución más fácil y mejor usar una técnica similar para hacer los puntos que obtenemos de ContourPlot más precisa. Comience desde el punto inicial, luego avance a lo largo del gradiente hasta que crucemos 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 de los puntos de extracción 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 {x, y} puntos)

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 discutió here, una article por Paul Abbott en los Mathematica Journal (Roots encontrando en un intervalo) da los siguientes dos métodos alternativos para la obtención de una lista de {x, y} valores de ContourPlot, 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]; 

donde

ptsXY2 == ptsXY == ptsXY3