从Mathematica中的CONTOURPOT提取轮廓

Question

我有两个变量的函数f（x，y），我需要知道它交叉零的曲线的位置。 ContourPlot非常有效地做到这一点（也就是说：它使用巧妙的多机方法，而不仅仅是蛮力细粒扫描），而只是给了我一个情节。我想拥有一组值{x，y}（带有指定的分辨率）或一些插值功能，使我能够访问这些轮廓的位置。

已经想到了从Contourplot的完整形式中提取此内容，但这似乎有点黑客。有什么更好的方法吗？

Answer 1

如果您最终从中提取点 ContourPlot, ，这是一种简单的方法：

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

Answer 2

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

编辑

如所讨论的 这里, ， 一个 文章 保罗·雅培（Paul Abbott） Mathematica杂志 (间隔寻找根部）提供以下两种替代方法，用于从Contourplot获得{x，y}值列表，包括（！）

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