13

2つの変数の関数f（x、y）があり、ゼロと交差する曲線の位置を知る必要があります。ContourPlotはそれを非常に効率的に実行します（つまり、ブルートフォースのきめ細かいスキャンだけでなく、巧妙なマルチグリッド法を使用します）が、プロットを提供します。値のセット{x、y}（特定の解像度）またはこれらの輪郭の位置にアクセスできるようにする補間関数が必要です。

ContourPlotのFullFormからこれを抽出することを考えましたが、これはちょっとしたハックのようです。これを行うためのより良い方法はありますか？

4

2 に答える 2

13

からポイントを抽出することになった場合ContourPlot、これは簡単な方法の1つです。

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

編集

ContourPlot非常に正確な輪郭を生成しないようです。もちろん、これらはプロット用であり、そのためには十分ですが、ポイントは等高線上に正確に配置されていません。

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}

独自の方法で輪郭をトレースすることもできますが、一般的な方法で行うのは大変です。滑らかな輪郭を持つ滑らかに変化する関数に対して機能する概念は次のとおりです。

  1. ある点（pt0）から開始し、の勾配に沿った等高線との交点を見つけますf。

  2. これで、輪郭にポイントができました。一定のステップ（resolution）で輪郭の接線に沿って移動し、ステップ1から繰り返します。

シンボリックに区別できる関数でのみ機能する基本的な実装は次のとおりです。

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]

輪郭検出方法の図

赤い点は「開始点」であり、黒い点は輪郭のトレースです。

編集2

おそらく、同様の手法を使用して、得られるポイントをContourPlotより正確にする方が、より簡単で優れたソリューションです。最初のポイントから開始し、輪郭と交差するまでグラデーションに沿って移動します。

この実装は、シンボリックに区別できない関数でも機能することに注意してください。これが当てはまるかのように関数を定義するだけf[x_?NumericQ, y_?NumericQ] := ...です。

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
于 2011-07-24T09:42:23.887 に答える
8

からポイントを抽出するためのわずかなバリエーションContourPlot（おそらくDavid Parkによる）：

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

または（{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];

編集

ここで説明するように、 MathematicaJournalのPaulAbbottによる記事（間隔内のルートの検索）では、ContourPlotから{x、y}値のリストを取得するための次の2つの代替方法が示されています。

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

上記の例の場合

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

と

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

どこ

ptsXY2 == ptsXY == ptsXY3
于 2011-07-25T11:34:26.800 に答える