Аналитически рассчитать всё ты не в силах, И подгонять решения не можешь ты, но сможешь сделать результат красивым, осуществив все вычислимые мечты

Ранее в одном из постов уже было сказано, что полиномы степени 5 и выше в общем виде неразрешимы в радикалах, а значит и множество всех точек полиномиальной поверхности F(x,y,z) = 0 пятой степени и выше аналитически вычислить невозможно. И тут методы приближённых вычислений должны сказать; «Не смогли смириться с поражением? И куда вас это привело? Снова ко мне!». Представим, что нам необходимо найти точки, где полиномиальная функция обращается в нуль. Можно взять диапазон построения [a, b] и каждую ось разбить на N отрезков, для каждой точки найти F(xi, yj, zk) и проверить меняет ли функция знак при переходе от исходной точки к соседним. Это самое первое, что приходит в голову, но подобный метод медленный, сложный и требует солидных вычислительных мощностей при незначительном увеличении точности. При этом, нельзя не заметить, что мы вычисляем значение функции по всему объёму, а в результате должны получить точки на двумерной поверхности, а значит бОльшая часть вычислений функции уходит «в молоко». Тогда, для экономии времени следует вначале построить грубую сетку, а между узлами, где произошла смена знака увеличить концентрацию узлов.

Так, очень странно, вроде всё рассказал, а текст продолжается. А, точно! Надо рассказать об анализе поведения самой полиномиальной функции для построения качественной адаптивной сетки. Так как эта задача, требующая кОмплексного решения и творческого подхода к проблеме, то будем, есть этого слона постепенно. На уровне грубой сетки можно определить количество компонент связности. Для подсчёта компонент связности используют алгоритмы поиска букетов графа (например, с помощью методов "обхода в ширину" или "союз-находка". Названия странные потому что использован дословный перевод). Обычно рассматривают 6-, 18- или 26-связность в 3D (по граням, рёбрам или вершинам). Рассмотрим подробнее выше упомянутые методы «обхода в ширину» и «союз-находка» для полиномиальной поверхности F = x^7 + y^7 + z^7 - 4(xy^2z^4 + zx^2y^4 + yz^2x^4) + 3(xy^3z^3 + x^3yz^3 + x^3y^3z)

1) Обход в ширину (BFS, Breadth-First Search). Как же я люблю BFS. Вот он, слева-направо: Метод из теории графов, который помогает определить количество компонент связности. В этом алгоритме выбирается начальная вершина. Далее в память заносятся все её соседи, потом соседи соседей и т.д. до тех пор, пока все связанные вершины не будут посещены. А потом берутся вершины, которые не были посещены и цикл повторяется заново, и так до тех пор, пока все вершины графа не будут перечислены во всех связанных компонентах.

2) Второй метод «Союз-находка» (Union-Find, DSU) похож на BFS, выполненный в обратном порядке. Изначально каждую вершину мы считаем отдельной компонентой. Далее, для каждой точки проверяем: не происходит ли смена знака при переходе на соседнюю вершину. После прохождения цикла по всем узлам, и расчёта отметаются компоненты связности, имеющие 1 вершину.

Это неплохой метод, позволяющий увидеть общую структуру алгебраической поверхности, но если мы выберем такой шаг сетки, что случайно пропустим отдельную компоненту связности между узлами, то обнаружить её этим методом мы не сможем. Как видно на примере, показанном на рисунке выше, по приближённым расчётам поверхность имеет 5 компонент связности. Синим цветом показан граф, полученный методом BFS, красным цветом DSU. Как видно по рисунку, эти методы дают одинаковый результат, что говорит о том, что грубой сетки достаточно, чтобы получить общий вид структуры алгебраической поверхности вне зависимости от выбранного метода. Мы видим, что «кубики», в которых содержатся отдельные компоненты связности лежат на одной прямой (стоит сказать, что уравнение этой прямой x = y = z). Это может быть как отдельные компоненты связности, так и сингулярность, так и артефакт приближённых расчётов. Без более тонкой сетки и различных методов анализа отдельных компонент связности в этих участках однозначный вывод сделать нельзя. К слову, при использовании приближённых методов и грубых сеток также сложно различить такие структуры, как каспы и точки ветвления от артефактов.

Продолжаем уходить вглубь адаптивных сеток. Как было сказано выше: при смене знака между узлами, мы исследуем поведение поверхности более детально. Между узлами грубой сетки более точную сетку можно локально ориентировать вдоль градиента функции F(x,y,z), что даст больше возможности для обнаружения каких-либо структур алгебраической поверхности. Но не только градиентами едины. Стоит описать принципы вариационных методов адаптивных сеток. Общая схема решения этим методом представлена на фотографии, приложенной к посту, но если подробнее:

Как было сказано выше алгебраическая поверхность это двумерная поверхность, а значит, множество точек, принадлежащие ей можно описать следующим образом: ?(u, v) = (x(u, v), y(u, v), z(u, v)). На практике этот факт реализуется следующим образом: предположим, что мы уже построили грубую сетку и получили местоположение различных компонент связности и в каждом элементе объёма, в котором должна находится поверхность, мы строим трёхмерную сетку, разбивая её таким образом, что чем больше проекция градиента ?F(x,y,z) на ось смещения, тем ближе между собой располагаются узлы сетки. Выше представлен рисунок, показывающий участок алгебраической поверхности, построенный на равномерной сетке, адаптивной по направлению градиента функции и третья сетка построена на основе величины частной производной по переменным (x,y,z).

Следующим шагом в применении вариационного метода задаём метрику для выбранного уравнения поверхности. Для приближённого вычисления точек алгебраической поверхности и при невозможности аналитически определить параметризацию x(u, v), y(u, v), z(u, v), вблизи алгебраической поверхности метрику можно задать следующим образом: M = exp(-F^2/?^2)*(E + ?(?F??F)/(?F^2 + ?)), где ? определяет ширину слоя вокруг нулевой поверхности, ? — степень анизотропии, ? - некий произвольный положительный параметр чуть больше нуля, E – единичная матрица 2*2 (ну или символ Кронекера, если вам больше нравится язык тензоров, а не матриц). На самом деле, несмотря на обилие параметров в формуле, большую их часть мы можем выбрать произвольно, что называется в пределах нормы. Если у нас есть метрический тензор и преобразования координат, то не грех и построить матрицу Якоби J для локального преобразования координат (в некоторых случаях для приближённой сетки во главу вычислений ставят Гессиан). Далее, для каждой точки уточнённой сетки мы ищем те, где минимизируется функционал, зависящий от метрического тензора и Якобиана. Если некоторые детали поверхности после построения уточнённой сетки всё ещё сложно различимы, то повторяем весь выше описанный алгоритм уточнения сетки. Но в нашем разговоре о вариационном методе мы упустили хедлайнера нашей вечеринки: сам функционал, который необходимо минимизировать. И тут есть разные подходы, как учитывающие особенности и инварианты заданной алгебраической поверхности, так и более универсальные. Затронем некоторые из них.

1) Функционал Дирихле. В выделенном объёме мы выбираем точки минимизирующие функционал следующего вида:

E = ??(||?x||^2 +||?y||^2 + ||?z||^2)*?(u, v)*du*dv. Где ? – весовая функция. Обычно, в качестве весовой функции служит величина ?(u, v) = 1 + ?*||??(x(u,v),y(u, v),z(u,v))||, ? – положительный параметр, регулирующий степень адаптации. И снова мы имеем дело с параметром, который можно выбрать в пределах разумного и не прогадать.

2) Функционал Уинслоу. Да, мышь из «Котопса» тоже умеет в этом вашем матанализе. Но если серьёзно: Функционал Уинслоу минимизирует искажения и обеспечивает равномерное распределение узлов, особенно в областях с большими градиентами функции. В основе лежит идея минимизации энергии, аналогичной энергии диффузии.

I(u,v) = ?? (||?u||^2 + ||?w||^2)/?(x, y, z)dx*dy*dz, где ? – область интегрирования, то есть сама алгебраическая поверхность или в случае приближённых вычислений её участок, на котором осуществляется повышение точности грубой сетки, ну и весовая функция может быть той же самой, что и в предыдущем примере (с поправкой на то, что зависеть она будет от других переменных). Функционал Уинслоу минимизируется по всем возможным отображениям сетки, что приводит к системе эллиптических уравнений (уравнения Уинслоу). ?*( ?(x, y, z)^-1*?u) = 0, ?*( ?(x, y, z)^-1*?v) = 0. Решение этой системы определяет оптимальное распределение узлов.

3) Функционал Хуанга. Применять этот функционал для адаптивных сеток сложнее, но зато его точность выше, чем у Уинслоу и Дирихле. ?? dx*dy*dz* p^-1*tr((J^-1*M*J(T))^p) – log(det(J^-1*M*J(T)))), где J(T) обозначается транспонированная матрица Якоби, p – параметр, который обычно берут равным 1 или 2. Этот функционал минимизируется по всем возможным отображениям, что приводит к оптимальному распределению узлов с учётом геометрии задачи. Он Совмещает критерии равномерности и выравнивания. Функционал Хуанга также обеспечивает, чтобы элементы сетки были как можно ближе к равномерным (по объёму в метрике M) и как можно более изотропными (без сильных искажений формы). При правильной реализации функционал Хуанга не допускает вырождения ячеек, а также позволяет учитывать произвольные метрические тензоры, что важно для анизотропных и адаптивных задач. Этот метод больше подходит для приближённого решения дифференциальных уравнений с граничными условиями в области с нестандартной геометрией, но в вопросе выбора функционала никто нам не запрещает поджигать спички лазером.

На графике выше представлена алгебраическая поверхность для нашего уравнения F, полученная с помощью вариационного метода адаптивных сеток с использованием функционала Дирихле. И на данном графике так же видны 5 отдельных компонент связности. Построим градиент функции F в областях малых компонент связности. В файле, прикреплённом к данному посту показано изображение проекции векторного поля на плоскость, в которой лежат компоненты связности. На данном рисунке видно, что для компонент 1 и 2 (зелёная и синяя точки) векторы пресекаются и наслаиваются друг на друга, а для компонент 3 и 4 (жёлтая и фиолетовая точки соответственно) вдоль прямой наблюдается пробор. Можно предположить, что в области компонент 1 и 2 наблюдается сток векторного поля, а в компонентах 3 и 4 – источник. Такие точки, где градиент меняет направление и величину привлекают наше внимание потому что это может влиять на топологию и локальную структуру поверхности. В некоторых областях топологии (например, в теории Морса) такие точки и вовсе играют одну из ключевых ролей в анализе заданной алгебраической поверхности. Но следует ответить на вопрос: что это за маленькие компоненты связности? И почему они наблюдаются и при повышении точности сетки? При этом, стоит отметить, что если взять поверхность не F(x,y,z) = 0 а F(x,y, z) = d, где d – некая величина, находящаяся в окрестности нуля, то мы увидим уже не отдельные компоненты связности, а тонкую цилиндрическую поверхность вдоль линии, где раньше располагались отдельные компоненты связности (в зависимости от знака d линия будет направлена или от нуля в положительном направлении вектора (1, 1, 1), как это показано на одном из графиков, или антипараллельно ему). В данном случае это означает, что нулевой уровень функции "разорван" — поверхность не соединена, а разбита на отдельные компоненты, возможно, из-за особенностей самого уравнения или его параметров. Ведь не стоит упускать из виду то, что алгебраические поверхности могут иметь сложную топологию: нулевой уровень может быть "разорван" на изолированные друг от друга части. Здесь ещё следует учесть, что связность уровня функции может резко меняться при малом изменении параметра, особенно вблизи особых точек или при наличии "тонких горлышек".

Подведём итоги: Методы конечных элементов - это мощной инструмент для получения графического представления изучаемого объекта. А тема приближённых решений настолько обширна, что все вопросы и аспекты этой области и упомянуть трудно. Но при всём их многообразии и широте области применения топология и алгебраическая геометрия даёт нам более точные и детальные представления о алгебраических многообразиях, которые невозможно получить сколь угодно точной адаптивной сеткой и т.п. Можно сказать, что приближённые решения служат приятным дополнениям к точным математическим изыскания, но никогда их не заменят. Обобщая эти рассуждения на более широкий класс задач, я бы сказал так: аналитическое и приближённое решение это своего рода Инь и Янь между физиками и математиками, где в результате их борьбы и объединения рождается объективная истина.

Источник: vk.com

Комментарии: