Фундамент итерационной математики

В самых первых статьях нашего цикла мы рассматривали метод дихотомии и метод Ньютона для поиска корней нелинейных уравнений вида f(x) = 0. Метод Ньютона быстр, но требует вычисления производных и часто расходится при плохом выборе начального приближения. Существует еще один, невероятно элегантный и концептуально прозрачный класс алгоритмов — метод простой итерации (или метод последовательных приближений). Его глубокая математическая суть базируется не на геометрии касательных, а на фундаментальных понятиях топологии и функционального анализа.

Чтобы применить этот метод, исходное уравнение f(x) = 0 необходимо алгебраически преобразовать к эквивалентному виду: x = phi(x). Это можно сделать бесконечным числом способов (например, просто прибавив x к обеим частям, или выразив x из сложного логарифма). Теперь наша задача поиска корня превратилась в задачу поиска неподвижной точки отображения phi (fixed point). Неподвижная точка — это такое значение x, которое при подстановке в функцию phi выдает на выходе в точности то же самое значение x. Вычислительный алгоритм становится тривиальным: мы берем начальное случайное приближение x0, вычисляем x1 = phi(x0), затем x2 = phi(x1) и так далее до бесконечности.

Сжимающие отображения и теорема Банаха

Главный вопрос вычислительной математики: при каких условиях эта бесконечная последовательность x_n действительно сойдется к истинному корню, а не улетит в бесконечность и не зациклится? Ответ дает знаменитая теорема Банаха о неподвижной точке (принцип сжимающих отображений), доказанная польским математиком Стефаном Банахом в 1922 году.

Теорема гласит: если функция phi(x) является сжимающим отображением на замкнутом интервале (то есть расстояние между любыми двумя точками после применения функции строго уменьшается с коэффициентом сжатия q < 1), то на этом интервале существует ровно одна неподвижная точка, и метод простой итерации гарантированно сойдется к ней из абсолютно любого начального приближения! С точки зрения дифференциального исчисления, условие сжатия означает, что модуль производной функции phi(x) в окрестности корня должен быть строго меньше единицы. Чем ближе производная к нулю, тем меньше коэффициент сжатия q, и тем феноменально быстрее алгоритм будет сходиться. Если же производная по модулю больше единицы (отображение растягивающее), алгоритм неминуемо разойдется, как бы близко к корню мы ни стартовали. Искусство применения метода простой итерации заключается именно в таком хитром алгебраическом преобразовании f(x) = 0 в x = phi(x), чтобы искусственно сделать производную phi(x) максимально близкой к нулю.

Коллективный интеллект на службе математики

В области глобальной стохастической оптимизации генетические алгоритмы (ГА) и дифференциальная эволюция (DE), основанные на теории естественного отбора, десятилетиями удерживали лидерство. Однако в конце 1990-х годов возникла мощная альтернативная парадигма эвристического поиска, вдохновленная не эволюцией видов, а социальным поведением животных. Как стая птиц находит лучшее место для кормежки в огромном лесу? Как муравьи без единого лидера строят оптимальные кратчайшие маршруты от муравейника к источнику пищи? Оказалось, что математическая формализация этих природных механизмов дает потрясающе быстрые численные алгоритмы оптимизации сложных недифференцируемых функций.

Этот класс алгоритмов получил название Swarm Intelligence (Роевой интеллект). В отличие от генетических алгоритмов, где особи (решения) безжалостно уничтожаются и заменяются мутантами (дарвиновская философия), в алгоритмах роевого интеллекта элементы популяции не умирают. Они сохраняют свою жизнь на протяжении всего процесса вычислений, перемещаясь в пространстве параметров и обмениваясь социальной информацией со своими соседями. Эта кооперативная природа алгоритмов делает их сходимость более плавной и управляемой.

Оптимизация роем частиц (Particle Swarm Optimization, PSO)

Алгоритм PSO был разработан в 1995 году Кеннеди и Эберхартом. Пространство поиска заполняется «частицами» (математическими векторами возможных решений). Каждая частица имеет свою координату и вектор скорости. На каждом шаге частица оценивает значение целевой функции в своей точке. Направление ее следующего шага (новая скорость) вычисляется как векторная сумма трех составляющих.

Первая составляющая — инерция: частица склонна двигаться туда, куда летела на предыдущем шаге. Вторая составляющая — личная (когнитивная) память: частица всегда помнит ту точку пространства, где она лично находила наилучший результат (наименьший минимум), и ее притягивает к этой точке. Третья составляющая — социальный опыт: частицу притягивает к той точке, которая является наилучшей среди всей стаи (глобальный лидер). Коэффициенты доверия к себе и к социуму являются настроечными параметрами. В результате этого простого математического взаимодействия частиц возникает невероятно сложный и целенаправленный танец роя, который быстро стягивается в самые глубокие минимумы сложных ландшафтов (функций Розенброка, Растригина и др.). PSO превосходно подходит для настройки весов в системах управления и проектирования СВЧ-антенн.

Муравьиные алгоритмы (Ant Colony Optimization, ACO)

Если PSO отлично работает в непрерывных пространствах, то для сложнейших задач дискретной комбинаторной оптимизации (например, задача коммивояжера, разводка печатных плат, маршрутизация пакетов в интернете) Марко Дориго в 1992 году создал Муравьиные алгоритмы (ACO). Идея заимствована из механизма стигмергии: непрямой коммуникации муравьев через отложение химических следов (феромонов).

Алгоритм запускает в граф (сеть дорог) популяцию искусственных муравьев. Каждый муравей строит свой маршрут вероятностно. Выбор следующего узла зависит от двух факторов: эвристической видимости (короткие ребра предпочтительнее длинных) и концентрации феромона на этом ребре. Построив маршрут, муравей откладывает на нем виртуальный феромон: чем короче (лучше) получился итоговый маршрут, тем больше феромона остается на его гранях. В то же время феромоны на всех гранях графа подвергаются математическому испарению на каждой итерации. Это предотвращает преждевременную сходимость к локальному минимуму (забывание плохих путей). Благодаря положительной обратной связи, ребра, принадлежащие глобально оптимальному маршруту, быстро накапливают огромную концентрацию феромона, и весь алгоритмический муравейник дружно сходится к идеальному решению (найденному кратчайшему пути).

Математика таяния льда и кристаллизации металлов

Процессы с фазовыми переходами первого рода (плавление льда, затвердевание отливок в металлургии, рост монокристаллов из расплава, замораживание грунтов при строительстве) представляют собой сложнейший вызов для вычислительной математики. Математически эти процессы описываются классической задачей Стефана (сформулированной Йозефом Стефаном в 1889 году при изучении замерзания полярных льдов). Основная проблема задачи Стефана — наличие четкой, геометрически тонкой границы раздела фаз (фронта плавления или кристаллизации), которая физически перемещается в пространстве с течением времени.

На этой подвижной границе происходит выделение или поглощение скрытой теплоты фазового перехода. Математически на фронте ставится жесткое условие теплового баланса (условие Стефана): скорость движения границы прямо пропорциональна разности тепловых потоков, приходящих из твердой и жидкой фаз. Трудность заключается в том, что положение этой границы заранее неизвестно — оно является частью искомого решения (задача со свободной границей). Решать такие дифференциальные уравнения классическими методами сеток крайне неудобно, так как фронт постоянно пересекает узлы сетки, вызывая сингулярности в расчетах.

Методы явного отслеживания фронта (Front Tracking)

Исторически первым подходом были методы с отслеживанием фронта. В них расчетная область разделяется на две отдельные зоны: твердую и жидкую, в каждой из которых решается свое линейное уравнение теплопроводности Фурье. На границе между ними сетка динамически перестраивается (или деформируется) на каждом временном шаге так, чтобы узлы сетки всегда строго лежали на фронте фазового перехода.

Этот подход (например, с использованием ALE сеток) дает высочайшую точность определения формы кристалла и распределения температур. Однако он терпит сокрушительное фиаско, если топология фронта начинает меняться. Если в процессе литья в форме образуется несколько независимых кусков льда, которые затем сливаются (или наоборот, одна жидкая ванна распадается на части), алгоритм отслеживания сетки запутывается в узлах, порождает математические ошибки и часто останавливает программу. Для сложной 3D-металлургии этот подход оказался слишком хрупким.

Метод энтальпии: сквозной счет без выделения границы

Чтобы обойти топологические кошмары, вычислительная теплофизика разработала гениально простой и абсолютно надежный класс алгоритмов — методы сквозного счета, самым популярным из которых является Метод энтальпии (Enthalpy method). В этом подходе мы отказываемся от явного отслеживания резкой границы между льдом и водой. Вместо этого вводится понятие эффективной теплоемкости и полной энтальпии (теплосодержания) материала.

Математически скрытая теплота плавления не выделяется как сингулярный источник на границе, а искусственно «размазывается» (сглаживается) по небольшому температурному интервалу (например, между 0 и 1 градусом Цельсия). В этом искусственном интервале (называемом двухфазной зоной или mushy zone) эффективная теплоемкость материала становится просто гигантской. Теперь мы решаем одно единое, нелинейное уравнение энергии для всей расчетной области (и для твердой, и для жидкой фазы одновременно), вообще не задумываясь о том, где находится фронт! Позиция фронта легко восстанавливается после расчетов просто путем поиска изотермы, соответствующей температуре плавления. Метод энтальпии легко справляется с любыми топологическими изменениями (слиянием и распадом зон кристаллизации), легко программируется на фиксированных декартовых сетках и сегодня является стандартом во всех литейных CAE-модулях (например, в ProCAST или MAGMASOFT).