Алгоритм/принцип построения графика неявной зависимости f(x,y) = 0?

Question

[limon_spb]

Хабралюди! Подскажите, подкиньте идей!

Нарисовать на canvas, c++, pascal и т.п. график функции y = f(x) — проще простого. Выясняем ОДЗ аргумента и с определенным шагом вычисляем y(x). Ставим точку, и так далее.

А как быть с графиком f(x,y) = 0, если y(x) не выразить? Понимаю, что в общем случае, наверное, можно решить только перебором значений по плоскости. Какие есть частные решения? Нужен что-то вроде brain-storm. Заранее спасибо!

Начну список:

• Перебор значений по x и по y с определенным шагом. Долго, но 100% :-)
  
  
• Нахождение зависимости y(x). Подходит только тогда, когда можно выразить.
  
  
• Переход в полярные координаты, выражение радиуса от угла. Пробегание по углу от 0 до 2pi. Подходит не всегда, но бывает, когда y(x) выразить сложно.
  
  
• Tupper, Jeff. «Reliable Two-Dimensional Graphing Methods for Mathematical Formulae with Two Free Variables» www.dgp.toronto.edu/people/mooncake/papers/SIGGRAP...
  
  Спасибо gribozavr. (не получилось почему-то в вопросе вставить ссылкой)
  
  
• поиск нулей функции градиентным методом (спуском)
  
  Спасибо Lerq
  
  

Answer 1

[Sergey Lerg]

Получается что вам нужно найти пересечение поверхности z = f(x, y) с плоскостью z = 0.
Я знаю сюществуют специальные численные методы для нахождения этих контуров пересечения, нужно хорошенько в сети их поискать.

Comments

[Sergey Lerg]

Если из чего-то серьезного, то вот нашел что-то
citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.953&rep=rep1&type=pdf

Если влоб, то можно воспользоваться градиентным методом. Отображаемая площадь делится на N точек (выше N — больше качество и трудозатратнее). Точки выбираются либо по рандому, либо сеткой с определенным шагом. Далее для каждой точки вычисляется ее градиент. Зная вектор градиента мы знаем примерно в какую сторону нужно перемещать точку, чтобы приблизить ее к z = 0. Величина шага перемещения точки тоже настраивается и зависит от z (ближе — медленнее). Точка перемещается и вся процедура повторяется до тех пор пока либо не достигнем z = 0 с заданной точностью, либо других условий останова (выход за пределы отображаемой области, лимит количества итераций).
Таким образом случайно выбранные точки с течением времени будут перемещаться к контурам пересечения с z = 0 и образовывать эти самые линии.

[limon_spb]

Спасибо за совет.
Возможно я мыслю поверхностно, но мне кажется, что решать 2Д задачку через 3Д — это как-то неординарно по меньшей мере… Даже не знаю, как искать.
В принципе, по первым курсам универа я помню, мы «рассекали» различные поверхности плоскостями, но получались контуры, которые так же неявно заданы. Например x^2+y^2+z^2 = 4 пересекаем плоскостью z = 0, получаем x^2 + y^2 = 4.

[Sergey Lerg]

В конце еще близкие точки можно соединить линиями, если заметен промежуток между ними.

[limon_spb]

Спасибо еще раз. Нужно с этим методом разобраться!

[limon_spb]

Немного поразмыслил…
> Зная вектор градиента мы знаем примерно в какую сторону нужно перемещать точку, чтобы приблизить ее к z = 0

Вектор градиента не факт, что ведет к z=0. Ведь z=0 не обязано быть локальным минимумом f(). Кроме того, f(x,y) при разных x и y может быть как положительным, так и отрицательным при одном и том же градиенте…

[Sergey Lerg]

Локальные минимумы возможны, да. Это нормально и ничего страшного. Их полезно определять и останавливать вычисления в данной точке.
Если f(x,y) положительна, то нам нужен отрицательный градиент (на уменьшение z), а если f(x,y) отрицательна, то нам нужен положительный градиент (на увеличение z). Если знак градиента не тот, что нам нужен, то просто идем в противоположную сторону. По модулю градиенты могут быть равны, но нужно на знаки смотреть. Возможна только одна ситуация, где градиент нам ничего не говорит, это когда он равен нулю.

[limon_spb]

Да, вы правы! Метод безусловно имеет право на существование, не смотря на свои минусы. Наверное неимоверно долго выполняется. Для каждой точки нужно выполнить достаточно много итераций градиентного спуска. Может так случайно получится, что все точки сойдутся лишь к маленькой дуге из всего большого графика. Может так случится, что много точек пропадет в локальных минимумах так и не дойдя до нуля.
Но в целом — очень интересно. Классно было бы визуализировать такое построение и наблюдать, как точки расползаются и принимают форму какой-то кривой!
Скажите, как бы назвать такой метод применения градиентного спуска? Внесу в список :-) Я о таком, к своему стыду, не слышал.

[Sergey Lerg]

У этого метода достаточно широкое поле для оптимизаций. Например, подсмотрел в документе, сначала выбираем K точек (K << N) и смотрим их знаки, может получится так:
— ---++---
--+++--
---++-+
-------+
------++
Если есть точки как положительные, так и отрицательные, то значит пересечение с z=0 на этом участке гарантированно есть. Теперь достаточно выбрать N точек только в области границ плюсов с минусами.
Из-за того, что K достаточно мало, какая-то маленькая область может пройти незамеченной. Для этого для всех K точек вычисляем векторы градиента и сравниваем попарно соседние вектора на коллинеарность и однонаправленность векторным умножением. Если их направление примерно одинаковое, то все нормально, если же оно сильно отличается, то значит между этими двумя точками что-то есть еще интересное и тоже требуется выбрать точек из числа N в этой областе.
После того как N точек выбрано, можно начать градиентный спуск для всех точек, но функцию выбрать не f(x,y), а |f(x,y)| — ее модуль. На пересечение с нулем это никак не повлияет, но вычисления облегчит. В этом случае задача сводится к обыденному поиску минимумов функции, они причем в большинстве своем будут острыми.
Метод можно назвать «поиск нулей функции градиентным методом (спуском)».

[limon_spb]

С модулем — хорошо. А что если брать f^2(x,y)? Тогда поиск нулей функции будет тем же самым почти, что поиск минимумов, где градиент — ноль.

[Sergey Lerg]

Можно и квадрат. Только дополнительная операция умножения.

[bfDeveloper]

Лучше взять модуль, а не квадрат и запустить любой из методов многомерной глобальной оптимизации. Для одномера подходит метод ломаных, который сходится ко всем минимумам. Для многомерного случая он называется метод Пиявского. Из универского курса могу вспомнить метод Хука-Дживса, метод Стронгина (один из лучших, но сложноват), метод деления на 3. Выбирайте любой метод оптимизации, который ищет все минимумы. Думаю, что оптимальнее эту задачу вы не решите.

Answer 2

[Akson87]

Могу предложить немного изменить описанный выше подход с градиентами и 3хмерным описанием проблемы. Возьмите срезы по XZ для каждого Y, который соответствует пикселям на экране (надо помнить, что нам не надо расчитывать результат с точность, которая превышает разрешение монитора). Тогда каждый срез будет представлять из себя стандартную задачу по решению уравнения f(x)=0, тут уж можно какие-нибудь численные итеративные методы стандартные (метод ньютона например или что-то другое) использовать (которые по сути и используют градиенты в каком-то смысле). Проблемой с такими срезами будет то, что будут плохо отображаться горизонтальные линии. Можно попробовать последовательно делать с фиксированными Х и У, тогда должно получиться получше. Но все это мысли вслух так сказать, как бы я подходил к проблеме.

Answer 3

[gribozavr]

Наиболее хороший с моей точки зрения метод
Tupper, Jeff. «Reliable Two-Dimensional Graphing Methods for Mathematical Formulae with Two Free Variables» www.dgp.toronto.edu/people/mooncake/papers/SIGGRAPH2001_Tupper.pdf

Comments

[limon_spb]

Отличная статья! Спасибо!

[Mrrl]

Возможно, я чего-то в алгоритме не понял. Я попробовал его применить, чтобы построить график уравнения
2*x*x-3*x*y-2*y*y-x+2*y=0 на квадрате [0,1]x[0,1] с разрешением 1024x1024. При этом глубину разбиения ограничил 30 шагами (т.е. интервальному алгоритму разрешалось спускаться только до квадратов со стороной 2^(-30) ). Программа тут же заявила, что 30 шагов ей мало (остаются красные пиксели), а всего (даже для неполного результата) потребовалось проанализировать 258 миллионов квадратиков. Правда, нарисовать итоговую картинку я еще не пытался.
А у Вас этот алгоритм заработал?

[gribozavr]

Я не реализовывал, только статью читал. Вот результат GraphEq от авторов: www.habrastorage.com/images/aaa.png Много красного (точек, где нет уверенности).

[Mrrl]

После небольшой правки с первым уравнением программа справилась (перебрав 8.6 миллионов квадратов, размером до 2^(-21)). Но вот построить график 9*x*x-6*x+9*y*y-1.25=0 ей не под силу: окружность касается прямой y=1/2 в точке, не являющейся двоично-рациональной. Доказать, что в смежном квадратике (где y>=1/2)есть точка кривой, интервальный метод не сможет :(
Конечно, если результат «много красного» допустим и кажется более ценным, чем «пропущены части графика, которые меньше ячейки сетки», то алгоритм имеет смысл применять — он дает хорошую оптимизацию на огромных (по числу ячеек мелкой сетки) площадях. Но надо хорошо подготовить функцию f к вычислению в интервальной арифметике, там разный порядок вычислений может дать очень разные по качеству результаты.

[gribozavr]

По-моему лучше пользователю показать где именно точные результаты, а где неопределённость, чем выдать результат, а корректности которого судить никак нельзя.

[Mrrl]

Зависит от цели и ресурсов. Научно правильное решение легко может оказаться стрельбой из пушки по воробьям. И кроме ожидаемого воробья выдать пользователю еще и стаю убитых комаров, которые тому совсем не нужны. Если программа будет 15 секунд гарантированно строить все компоненты графика, а потом еще 5 секунд выделять и выбраывать компоненты размером меньше порогового, то это будет несколько хуже, чем если она сза 0.1 секунду построит основную компоненту каким-нибудь эвристическим методом.

Answer 4

[Mrrl]

Можно сделать так. Если известна хотя бы одна пара соседних точек растра, на которых функция имеет разные знаки (например, f(x,y)>0, f(x+1,y)<0), то линия проходит межу ними. Можно взять один из квадратов, ограничивающих этот отрезок (например, (x,y)-(x+1,y)-(x+1,y+1)-(x,y+1) ), и, проверив знаки функции f(x+1,y+1) и f(x+1,y), узнать, какую вторую сторону квадрата пересечет график. Процесс продолжать, пока график не дойдет до границы области или не замкнется.
Для поиска стартовых точек можно взять значения в точках (Nx,Ny) при каком-нибудь N, и уже на этом растре искать соседние точки с разными знаками: график наверняка пройдет между ними. Маленькие компоненты при этом могут потеряться, но для их поиска пришлось бы использовать очень серьезные методы — например, решать систему f(x,y)=0, df(x,y)/dx=0 (ее решения дадут экстремальные точки всех компонент). Можно для этого применить метод Ньютона (тоже стартуя из всех точек (Nx,Ny)).

Comments

[Mrrl]

А насчет того, что перебирать все (x,y) долго, возникает вопрос — а какие требования? F(x,y) долго вычислять? Или разрешение экрана очень высокое? Или нужно какое-нибудь особое быстродействие (перерисовка в реальном времени)? Потому что в остальных случаях вычислить 10 миллионов значений функций не должно быть проблемой.

[limon_spb]

Действительно, для каких-то задач перебор всех х, у подойдет, поэтому этот метод я не вычеркиваю совсем. Но в моем случае — нужна, как вы говорите, перерисовка в реальном времени. Более того — нужно уметь подбирать такую f(x,y), которая бы по методу наименьших квадратов проходила ближе всего к заранее заданному массиву точек.

[Mrrl]

Тогда всесовсем просто. Для заданных точек найдется две, для которых f будет разных знаков, где-то между ними — отрезок, на котором f меняет знак (ищется за логарифм, делением пополам). А дальше — обход компоненты, как я написал. Общее время работы — линейное от длины кривой.

[gribozavr]

И где вы будете брать «точки»? В середине, в левом нижнем/левом верхнем/… углу? И тогда пропустите все части графика, которые меньше ячейки сетки? Прочитайте статью из моего ответа, всё совсем не просто, но вполне реализуемо и пугать необходимостью решать такие системы как вы говорите не нужно.

[Mrrl]

Так написали же — «такую f(x,y), которая бы по методу наименьших квадратов проходила ближе всего к заранее заданному массиву точек». Точки я возьму из «заданного массива». При такой формулировке очевидно, что надо искать компоненту, которая проходит где-то между точками. Конечно, если пространство функций слишком замысловато, то и самая близкая кривая может оказаться «пеной» из многих компонент, но это будет означать, что пространство выбрано не совсем удачно.
А статью прочитаю. Как можно справиться с цилиндрической декомпозицией без решения систем, это уже интересно :)

Answer 5

[bfDeveloper]

Есть ещё одно интересное соображение по трёхмерному решению, навеянное методом Хука-Дживса.
Пусть мы нашли одну точку нашего графика (например, рассмотрев сечение x = const и применив в нём метод Пиявского). Назовём эту точку опорной. Если предположить, что мы ищем одну непрерывную кривую, то можно утверждать, что рядом с найденной точкой есть другие, принадлежащие этой кривой. Возьмём некоторый шаг h и посчитаем значение функции F(x, y) в нескольких точках на расстоянии h от первой. Выберем точку со значением максимально близким к 0 за следующую опорную и повторим действия. Таким образом мы как бы пройдёмся по дну оврага (если брать модуль F) или просто по склону, сохраняя высоту.
Если хорошо подумать над тем как искать следующую точку, то можно сообразить как от начальной точки пройти сразу по двум направлениям (она же не обязательно будет на границе, и описанный выше метод найдёт кривую только по одну сторону от точки).
Получив массив точек, мы можем применить метод градиентного спуска для каждой из них для получения более точного результата. Так как точки уже близки к кривой, практически исключается возможность скатывания в локальные минимумы.
Вообще в методах глобальной оптимизации принято сначала локализовывать точки минимума, а потом запускать методы локальной оптимизации в найденных областях. Здесь это тоже неплохо работает.

Answer 6

[limon_spb]

Спасибо всем за ответы! Жаль, как решение, можно только один пост поставить!