contorni estratto da ContourPlot in Mathematica

Question

Ho una funzione f (x, y) di due variabili, di cui devo conoscere la posizione delle curve in cui si annulla. ContourPlot lo fa in modo molto efficiente (che è: esso utilizza metodi multi-grid intelligenti, non solo una forza bruta a grana fine scan), ma appena mi dà una trama. Mi piacerebbe avere un insieme di valori {x, y} (con qualche risoluzione specificata) o forse una qualche funzione di interpolazione che mi permette di ottenere l'accesso alla posizione di questi contorni.

hanno pensato di estrarre questo dal FullForm di ContourPlot, ma questo sembra essere un po 'di hack. Un modo migliore per fare questo?

Answer 1

Se si finisce per l'estrazione di punti da ContourPlot, questo è un modo semplice per farlo:

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

Modifica

Sembra che ContourPlot non produce contorni molto precisi. Sono naturalmente significava per il tracciato e sono abbastanza buono per questo, ma i punti non si trovano esattamente i contorni:

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}

Possiamo cercare di venire con il nostro metodo per tracciare il contorno, ma è un sacco di problemi a farlo in modo generale. Ecco un concetto che le opere per le funzioni senza problemi diversi che hanno contorni morbidi:

1. Inizia da un certo punto (pt0), e trovare l'intersezione con il contorno lungo il gradiente di f.

2. Ora abbiamo un punto sul contorno. Muoversi lungo la tangente del profilo di un passo fisso (resolution), quindi ripetere dal punto 1.

Ecco un'implementazione di base che funziona solo con le funzioni che possono essere differenziati simbolicamente:

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]

I punti rossi sono i "punti di partenza", mentre i punti neri sono la traccia del contorno.

Modifica 2

Forse è una soluzione più facile e meglio usare una tecnica simile per fare i punti che otteniamo da ContourPlot più preciso. Iniziare dal punto iniziale, quindi spostare lungo il gradiente fino a che intersecano il contorno.

Si noti che questa implementazione funziona anche con funzioni che non possono essere differenziati simbolicamente. Basta definire la funzione come f[x_?NumericQ, y_?NumericQ] := ... se questo è il 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 leggera variazione per l'estrazione di punti da ContourPlot (forse a causa di 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 (come un elenco di {x, y} punti)

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

Modifica

Come discusso qui , un articolo da Paul Abbott nel Mathematica Journal (< em> Trovare radici in un intervallo ) ha pronunciato la seguente due metodi alternativi per ottenere un elenco di {x, y} valori da ContourPlot, tra cui (!)

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

Per l'esempio precedente

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

e

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

dove

ptsXY2 == ptsXY == ptsXY3