13

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

曾想过从 ContourPlot 的 FullForm 中提取它，但这似乎有点 hack。有没有更好的方法来做到这一点？

4

2 回答 2

13

如果您最终从中提取点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
于 2011-07-24T09:42:23.887 回答
8

ContourPlot从（可能是由于大卫公园）提取点的轻微变化：

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 杂志上的一篇文章（Finding Roots in an Interval）给出了以下两种从 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
于 2011-07-25T11:34:26.800 回答