我有兩個變量的函數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輪廓的要求,否則我認爲這是一個很好的解決方案。 –