我有两个变量的函数f(x,y),其中我需要知道曲线在哪个位置穿过零点。 ContourPlot可以非常有效地完成这个任务(即:它使用聪明的多网格方法,而不仅仅是一个强力细粒度扫描),但只是给了我一个情节。我想要有一组值{x,y}(具有某种指定的分辨率),或者可能有一些插值函数,这些函数允许我访问这些轮廓的位置。在Mathematica中从ContourPlot中提取轮廓
曾想过从ContourPlot的FullForm中提取这个,但这似乎有点破解。任何更好的方式来做到这一点?
我有两个变量的函数f(x,y),其中我需要知道曲线在哪个位置穿过零点。 ContourPlot可以非常有效地完成这个任务(即:它使用聪明的多网格方法,而不仅仅是一个强力细粒度扫描),但只是给了我一个情节。我想要有一组值{x,y}(具有某种指定的分辨率),或者可能有一些插值函数,这些函数允许我访问这些轮廓的位置。在Mathematica中从ContourPlot中提取轮廓
曾想过从ContourPlot的FullForm中提取这个,但这似乎有点破解。任何更好的方式来做到这一点?
如果你最终从ContourPlot
提取分,这是一个简单的方法来做到这一点:
points = Cases[
[email protected][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}
我们可以试着拿出我们自己的方法跟踪轮廓,但以一般方式做这件事很麻烦。下面是对具有平滑的轮廓平滑变化的功能工作的概念:
从一些点(pt0
)
开始,找到沿着f
梯度轮廓的交集。
现在我们在轮廓上有一个点。沿着固定的步骤(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 /. [email protected][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]
红点是“起点”,黑点是轮廓线。
EDIT 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[
[email protected][f[x, y] == 0, {x, -3, 3}, {y, -3, 3}],
Line[pts_] -> pts,
Infinity
]
refine[f, #] & /@ points
从ContourPlot
用于提取点A轻微变化(这可能是由于大卫公园):
pts = Cases[
ContourPlot[Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}],
x_GraphicsComplex :> [email protected], Infinity];
或(作为{X,Y}点的列表)
ptsXY = Cases[
Cases[ContourPlot[
Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}],
x_GraphicsComplex :> [email protected], Infinity], {x_, y_}, Infinity];
编辑
正如讨论的here,一个article由Paul雅培在数学杂志(在区间寻踪)给出用于获得{X,Y}由ContourPlot值,包括的一个列表中的以下两种可选方法(!)
ContourPlot[...][[1, 1]]
对于上面的例子
ptsXY2 = ContourPlot[
Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}][[1, 1]];
和
ptsXY3 = Cases[
[email protected][
Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}],
Line[{x__}] :> x, Infinity];
其中
ptsXY2 == ptsXY == ptsXY3
有用,谢谢! –
@Kasper,老实说,我认为没有比提取轮廓更容易的方法。如果您对精度/点位有特殊要求,或者您可以向我们提供有关您的功能的更多信息,那么我们可能会提出更好的解决方案,但这可能会更复杂。 – Szabolcs
这给你的f [x,y] == 1/2轮廓,而不是f [x,y] == 0轮廓的要求,否则我认为这是一个很好的解决方案。 –