Извлечение контуров из ContourPlot в Mathematica

У меня есть функция f(x,y) двух переменных, из которых мне нужно знать расположение кривых, на которых она пересекает ноль. ContourPlot делает это очень эффективно (то есть использует умные многосеточные методы, а не просто грубое мелкозернистое сканирование), но просто дает мне график. Я хотел бы иметь набор значений {x,y} (с некоторым заданным разрешением) или, возможно, некоторую интерполирующую функцию, которая позволяет мне получить доступ к местоположению этих контуров.

Думал извлечь это из FullForm ContourPlot, но это похоже на взлом. Есть ли лучший способ сделать это?


contour wolfram-mathematica zero numerical
person Kasper Peeters    schedule 24.07.2011    source источник

Ответы (2)


arrow_upward
12
arrow_downward

Если вы в конечном итоге извлекаете очки из 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
person Szabolcs    schedule 24.07.2011
comment
@Kasper, честно говоря, я не думаю, что есть более простой способ, чем извлечение контуров. Если у вас есть особые требования к точности / местоположению точки и/или если вы можете предоставить нам больше информации о своей функции, мы можем предложить лучшее решение, но оно, скорее всего, будет более сложным. - person Szabolcs; 24.07.2011
comment
Это дает вам контур f[x,y] == 1/2, а не контур f[x,y] == 0, как запрошено, в противном случае я думаю, что это прекрасное решение. - person Sjoerd C. de Vries; 24.07.2011
comment
@Sjoerd, это не имеет значения, я иллюстрировал концепцию. Если хотите, напишите Sin[x] Sin[y] - 1/2 == 0. - person Szabolcs; 24.07.2011
comment
@Szabolcs, да, это, безусловно, легко ;-) Меня немного беспокоила зависимость от внутренних представлений ContourPlot, но это решение настолько простое, что я с удовольствием его использую. - person Kasper Peeters; 24.07.2011
comment
@szabolc. Я понимаю это, но я думаю, что для иллюстративных целей вам лучше использовать f[x,y]==0, чтобы не усложнять излишне. Я бы сказал, что большинство людей увидят RHS как значение уровня, а LHS как спецификацию функции, хотя на самом деле это не обязательно. - person Sjoerd C. de Vries; 25.07.2011

arrow_upward
7
arrow_downward

Небольшой вариант извлечения очков из 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];

Изменить

Как обсуждалось здесь, статья Пола Эбботта в Mathematica Journal (Поиск Корни в интервале) предлагает следующие два альтернативных метода получения списка значений {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[
   Normal@ContourPlot[
     Cos[x] + Cos[y] == 1/2, {x, 0, 4 Pi}, {y, 0, 4 Pi}], 
   Line[{x__}] :> x, Infinity];

где

ptsXY2 == ptsXY == ptsXY3
person tomd    schedule 25.07.2011