Численные методы

Аппроксимация Паде: рациональные функции против рядов Тейлора

Ограничения полиномиального разложения Тейлора

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

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

Могущество дробно-рациональных функций Паде

В 1890 году французский математик Анри Паде (ученик знаменитого Эрмита) формализовал метод, который кардинально превосходит полиномы Тейлора по качеству аппроксимации. Аппроксимация Паде предлагает искать приближение функции не в виде одного многочлена, а в виде рациональной функции — отношения двух полиномов (числителя степени M и знаменателя степени N).

Идея построения аппроксимации Паде поразительно элегантна. Мы приравниваем исходную функцию (или ее известный ряд Тейлора) к искомой дроби. Затем умножаем обе части на знаменатель дроби. Приравнивая коэффициенты при одинаковых степенях независимой переменной слева и справа, мы получаем систему линейных алгебраических уравнений. Решив эту простую СЛАУ, мы находим искомые коэффициенты полиномов в числителе и знаменателе. Обозначается такая аппроксимация обычно как [M/N].

Преодоление полюсов и аналитическое продолжение

Наличие полинома в знаменателе наделяет аппроксимацию Паде математической «суперспособностью». Если исходная физическая функция имеет полюс (то есть обращается в бесконечность, как, например, функция тангенса при 90 градусах), ряд Тейлора ломается. Аппроксимация Паде, напротив, блестяще моделирует этот физический разрыв: знаменатель дроби просто сам обращается в ноль в этой точке, идеально повторяя топологию исходной функции!

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

Подробнее

Дробное исчисление: численные методы для производных нецелого порядка

Когда порядок производной становится дробным числом

В классическом математическом анализе, заложенном Ньютоном и Лейбницем, мы привыкли оперировать производными только целого порядка: первая производная (скорость), вторая производная (ускорение), третья производная (рывок) и так далее. Однако еще в 1695 году маркиз де Лопиталь задал Лейбницу в письме интригующий вопрос: «А что будет, если порядок производной n сделать равным 1/2?». Лейбниц ответил, что это приведет к удивительным парадоксам, из которых однажды будут извлечены великие следствия. Так зародилась совершенно невероятная ветвь математики — дробное исчисление (Fractional Calculus).

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

Операторы Римана-Лиувилля и Капуто

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

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

Методы Грюнвальда-Летникова и проблема вычислительной памяти

Поскольку аналитически решить дробные дифференциальные уравнения удается лишь в тривиальных случаях, возникает острая необходимость в эффективных численных методах. Самым интуитивным подходом к дискретизации является схема Грюнвальда-Летникова. Она напрямую обобщает классические конечные разности. Если первая производная вычисляется по двум точкам, вторая — по трем, то для вычисления дробной производной по Грюнвальду-Летникову требуется бесконечный биномиальный ряд.

В этом кроется главная вычислительная сложность («проклятие») дробного исчисления. Классическая производная локальна: чтобы найти скорость автомобиля прямо сейчас, нам нужно знать его положение только за последнюю миллисекунду. Дробная же производная глобальна. Это математический оператор с эффектом «затухающей памяти». Чтобы вычислить состояние дробной системы на n-ном временном шаге, компьютеру необходимо учесть всю историю состояния системы, начиная с самого первого шага t=0! Это приводит к тому, что время вычислений растет пропорционально квадрату количества шагов (O(N^2)), а требования к оперативной памяти растут линейно. Для длительного интегрирования таких систем математики вынуждены разрабатывать хитрые методы «короткой памяти» (усечения хвоста) и применять алгоритмы быстрых дискретных сверток на базе БПФ.

Подробнее

Автоматическое дифференцирование: скрытый двигатель алгоритмов машинного обучения

Почему не подходят символьные и численные методы?

Для оптимизации сложных функций (например, при обучении глубоких нейронных сетей с миллиардами параметров) нам необходимо быстро и сверхточно вычислять градиенты. Классическая вычислительная математика исторически предлагала два подхода, но оба оказались непригодны для задач современных масштабов. Первый подход — символьное дифференцирование (которым занимаются системы компьютерной алгебры типа Mathematica). Программа аналитически применяет жесткие правила дифференцирования к исходной формуле. К сожалению, для сложных алгоритмов с циклами ветвления этот метод приводит к феномену «экспоненциального разбухания выражений», когда итоговая формула производной занимает гигабайты памяти. Второй подход — численное дифференцирование конечными разностями. Как мы обсуждали ранее, этот метод страдает от катастрофической потери точности из-за неизбежных ошибок округления и требует N+1 вызовов функции для функции от N переменных, что при миллионах переменных замедлит расчеты до бесконечности.

Третий, поистине революционный путь — это Автоматическое Дифференцирование (АД). Это не символьная аналитика и не численное приближение. АД основывается на том простом факте, что любая, даже самая сложная компьютерная программа в конечном итоге состоит из последовательности примитивных арифметических операций (сложение, умножение) и базовых функций (синус, экспонента). АД применяет цепное правило дифференцирования (chain rule) к этим элементарным шагам прямо в процессе выполнения программы, сохраняя идеальную машинную точность и не создавая раздутых математических формул.

Прямой проход и алгебра дуальных чисел

Автоматическое дифференцирование имеет два основных режима работы. Прямой режим (Forward mode) концептуально основан на использовании так называемых дуальных чисел — специальной алгебраической структуры вида a + b*epsilon, где число epsilon обладает уникальным свойством: его квадрат строго равен нулю (по аналогии с комплексной мнимой единицей, квадрат которой равен -1). Применяя арифметику дуальных чисел к коду программы, мы одновременно за один единственный вычислительный проход получаем и итоговое значение функции, и ее точную направленную производную. Прямой режим невероятно эффективен и быстр, если мы анализируем функцию с малым числом входов и огромным числом выходов.

Обратный проход (Reverse Mode) и алгоритм Backpropagation

Однако в машинном обучении ситуация строго обратная: функция потерь (Loss function) — это всего лишь одно единственное скалярное число (один выход), которое сложнейшим образом зависит от миллионов весовых коэффициентов (входов). Для решения таких задач прямой режим потребовал бы миллионов повторных проходов. Здесь на сцену выходит обратный режим АД (Reverse mode), который в индустрии искусственного интеллекта больше известен под термином «алгоритм обратного распространения ошибки» (Backpropagation).

Обратный режим АД работает в два математических этапа. Сначала выполняется обычный прямой проход, во время которого алгоритм строит в памяти гигантский направленный вычислительный граф и бережно запоминает промежуточные значения абсолютно всех переменных (что требует колоссального объема оперативной видеопамяти VRAM). Затем алгоритм начинает двигаться по этому графу в обратном направлении, от выхода к входам, передавая сигналы чувствительности (градиенты). Благодаря гениальному применению цепного правила с конца, обратный режим позволяет вычислить точный градиент функции сразу по всем миллионам параметров всего за один обратный проход графа! Без этого потрясающего вычислительного алгоритма современный бум генеративных нейронных сетей был бы абсолютно невозможен.

Подробнее

Жесткие краевые задачи: метод многократной пристрелки по параметру

Когда одиночный выстрел уходит в бесконечность

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

Суть проблемы заключается в жесточайшей чувствительности дифференциальной системы к начальным данным. Если в задаче присутствуют быстрорастущие экспоненциальные компоненты, то даже микроскопическая ошибка (в миллиардных долях) при выборе начального угла «орудия» приведет к тому, что к середине расчетного интервала численное решение улетит в математическую бесконечность из-за переполнения регистров памяти компьютера. Алгоритм просто не сможет дойти до правого края, чтобы оценить ошибку промаха. Невозможно прицелиться, если любое малейшее дрожание руки приводит к отклонению снаряда на миллионы километров.

Метод многократной (параллельной) стрельбы

Чтобы обуздать экспоненциальный взрыв и восстановить вычислительную стабильность, математиками Келлером, Осборном и Деуфлхардом был разработан метод многократной стрельбы (Multiple Shooting Method). Идея этого метода заключается в том, чтобы не пытаться «прострелить» весь длинный и сложный интервал одним выстрелом.

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

Сшивка решений и матрица Якоби

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

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

Подробнее

Алгоритмы вычисления Быстрых Преобразований Фурье (БПФ) в многомерных пространствах

Алгоритм, изменивший цифровую эпоху

Преобразование Фурье — это математическая основа обработки любых сигналов. Оно позволяет раскладывать сложные волновые формы на набор простых синусоид. Дискретное преобразование Фурье (ДПФ) выполняет эту задачу для цифровых (сэмплированных) данных. Однако прямое, «школьное» вычисление ДПФ по формуле имеет квадратичную алгоритмическую сложность O(N^2). Для аудиофайла длиной в одну секунду (44 100 отсчетов) потребуется около двух миллиардов операций умножения. Если бы мы пользовались этим прямым методом, современная потоковая передача видео, мобильная связь 4G/5G и магнитно-резонансная томография (МРТ) были бы просто невозможны из-за вычислительных ограничений.

Проблема была решена в 1965 году, когда Джеймс Кули и Джон Тьюки опубликовали алгоритм Быстрого Преобразования Фурье (БПФ, или FFT). (Справедливости ради, похожие идеи высказывал еще Карл Фридрих Гаусс в 1805 году, но они опередили свое время). Алгоритм БПФ радикально снижает количество операций с O(N^2) до O(N log N). Для массива в миллион точек ускорение составляет почти 50 000 раз! Журнал Computing in Science & Engineering заслуженно включил БПФ в десятку величайших алгоритмов XX века.

Разделяй и властвуй: принцип бабочки

Самая популярная версия алгоритма Кули-Тьюки базируется на принципе «разделяй и властвуй» (divide and conquer) по основанию 2 (Radix-2). Алгоритм требует, чтобы количество точек N в массиве было степенью двойки (например, 256, 512, 1024). Если это не так, массив дополняется нулями (zero-padding).

Алгоритм рекурсивно разбивает задачу вычисления одного ДПФ размера N на вычисление двух ДПФ размера N/2 (одно для четных индексов массива, другое — для нечетных). Затем каждый из этих подмассивов снова делится пополам, и так далее, пока размер массива не достигнет 1. Затем результаты начинают объединяться (синтезироваться) обратно. На этапе объединения используется базовая вычислительная структура, которая из-за графа информационных потоков получила название «операция бабочки» (Butterfly operation). Она комбинирует два комплексных числа, используя умножение на специальные комплексные экспоненты (поворачивающие множители). Перед выполнением БПФ элементы массива переставляются в специальном порядке, называемом битовой инверсией (Bit-reversal permutation), что позволяет выполнять все вычисления прямо на месте (in-place), не выделяя дополнительную оперативную память.

Многомерное БПФ: от звука к изображениям

Обычное одномерное БПФ применяется для обработки звука или одномерных радиосигналов. Но как анализировать цифровые фотографии (поиск краев, сжатие JPEG, фильтрация шума) или решать трехмерные дифференциальные уравнения Пуассона на сетке? Для этого используется многомерное БПФ.

Счастье для инженеров заключается в том, что многомерное дискретное преобразование Фурье математически сепарабельно (разделимо). Это означает, что для вычисления двумерного БПФ изображения (матрицы пикселей) не нужно писать новый сложный алгоритм. Достаточно сначала применить обычное одномерное БПФ к каждой строке матрицы независимо. Получив промежуточную матрицу комплексных чисел, мы затем применяем то же самое одномерное БПФ к каждому ее столбцу. Этот строчный-столбцовый подход (Row-Column algorithm) легко масштабируется на 3D и 4D пространства, идеально распараллеливается на тысячи ядер графических процессоров (с помощью библиотек вроде cuFFT от NVIDIA) и является бьющимся сердцем современной компьютерной томографии и молекулярной химии.

Подробнее

Вычисление определителей и обращение матриц: численные аспекты

Разрыв между теоретической алгеброй и вычислительной реальностью

В курсе высшей алгебры студенты учатся решать системы линейных алгебраических уравнений (СЛАУ) вида Ax = b с помощью правила Крамера (через определители) или путем явного вычисления обратной матрицы (x = A^-1 * b). Эти методы красивы в теории и отлично подходят для ручного решения простейших систем 3x3 на листке бумаги. Однако в реальной вычислительной математике при написании программного обеспечения инженеры избегают этих методов как огня. Почему же аналитические фавориты терпят крах при встрече с кремниевыми процессорами?

Все дело в вычислительной сложности. Вычисление определителя по классической формуле разложения по строке (формула Лейбница) требует факториального количества операций (O(N!)). Для матрицы размером всего 20x20 вычисление определителя таким способом потребует 20! операций. Если суперкомпьютер будет выполнять миллиард операций в секунду, ему потребуется несколько сотен тысяч лет, чтобы решить эту крошечную задачу. Использование правила Крамера для СЛАУ требует вычисления N+1 таких определителей, что делает этот метод абсолютно бесполезным для практических расчетов в физике и инженерии.

Почему обращать матрицы явно — плохая идея?

Вторая распространенная ошибка новичков в программировании — явное вычисление обратной матрицы A^-1 для решения системы уравнений. Во-первых, вычисление обратной матрицы методом Гаусса-Жордана требует в 3 раза больше арифметических операций, чем решение самой системы методом исключения Гаусса (или с помощью LU-разложения). Мы тратим драгоценное время процессора впустую.

Во-вторых, обращение разреженных матриц приводит к катастрофическому заполнению памяти. Как мы помним из предыдущих статей, в методе конечных элементов матрица A на 99% состоит из нулей. Однако ее обратная матрица A^-1 почти всегда является полностью плотной (все элементы не равны нулю). Если исходная матрица из миллиона уравнений в формате CSR занимала несколько десятков мегабайт, то обратная матрица займет терабайты оперативной памяти. Именно поэтому золотое правило вычислительной линейной алгебры гласит: «Никогда не вычисляйте обратную матрицу явно, если вам нужно только решить систему уравнений». Вычисляйте LU-разложение или используйте итерационные методы.

Обусловленность и устойчивость обратной матрицы

Еще одной колоссальной проблемой является накопление ошибок округления. Каждая математическая матрица характеризуется числом обусловленности (Condition Number) — грубо говоря, это мера того, насколько матрица близка к вырожденной (к матрице с нулевым определителем).

Если число обусловленности велико (матрица плохо обусловлена), то микроскопические ошибки в исходных данных или ошибки округления с плавающей запятой внутри процессора при вычислении обратной матрицы будут многократно усилены (иногда в миллионы раз). Полученная компьютерная «обратная матрица» при умножении на исходную даст результат, крайне далекий от единичной матрицы. Анализ собственных и сингулярных чисел (SVD), а также использование алгоритмов предобуславливания (Preconditioning) стали обязательными инструментами для стабилизации процесса решения СЛАУ и борьбы с физической некорректностью инженерных моделей.

Подробнее

Численное моделирование в электродинамике: метод FDTD (сетка Йи)

Динамика электромагнитных полей в реальном времени

Система уравнений Максвелла является фундаментом всей современной классической электродинамики, описывая генерацию и распространение электромагнитных волн. Для расчета сложных антенн смартфонов, микроволновых печей, оптоволоконных кабелей и радарных систем-невидимок инженерам необходимо решать эти уравнения с высочайшей точностью. Одним из самых мощных и интуитивно понятных численных алгоритмов для этой задачи является Метод конечных разностей во временной области (Finite-Difference Time-Domain, FDTD).

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

Магия смещенной сетки Кейна Йи (Yee Grid)

Гениальность метода FDTD, предложенного Кейном Йи в 1966 году, заключается в особой топологии пространственной дискретизации. Если в обычном методе сеток все неизвестные величины (например, температура и давление) вычисляются в одних и тех же узлах, то в сетке Йи компоненты электрического (E) и магнитного (H) полей разнесены в пространстве.

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

Проблема открытого пространства: поглощающие условия (PML)

При моделировании излучающей антенны электромагнитная волна должна уходить в бесконечность. Однако память компьютера конечна, и мы вынуждены обрезать расчетную область математическими границами. Если просто оборвать сетку, волна отразится от этой границы, как от металлического зеркала, и разрушит вычисления внутри области.

Для решения этой фундаментальной проблемы Жан-Пьер Беренжер в 1994 году изобрел Идеально Согласованные Слои (Perfectly Matched Layer, PML). Это специальный математический фиктивный материал, который располагается по краям расчетной сетки. Волна любой поляризации, падающая на PML под любым углом, проникает в него без отражения на границе раздела, а внутри слоя ее амплитуда экспоненциально затухает до нуля. Изобретение PML сделало метод FDTD индустриальным стандартом в вычислительной электродинамике, оптике и фотонике.

Подробнее

Искусственный интеллект и численные методы: нейронные операторы и суррогатное моделирование

Сдвиг парадигмы в вычислительной математике

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

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

Нейронные операторы: обучение независимости от сетки

Фундаментальным недостатком классических нейронных сетей (например, сверточных CNN) в задачах физического моделирования является их жесткая привязка к дискретной сетке. Если сеть обучена на сетке 100x100 узлов, она не сможет работать на сетке 200x200 без полного переобучения. Настоящим прорывом стало создание нейронных операторов, самым известным из которых является Фурье-нейронный оператор (Fourier Neural Operator, FNO).

Нейронные операторы учатся отображать одно бесконечномерное функциональное пространство в другое. FNO использует теорему о свертке и выполняет основную часть вычислений не в физическом пространстве, а в спектральном пространстве Фурье. Сеть отбрасывает высокочастотный шум и обучается на низкочастотных гармониках, определяющих основную физику процесса. Результатом является алгоритм, который абсолютно независим от расчетной сетки (mesh-free). Обучив FNO на грубой сетке, инженер может запросить у сети предсказание на сверхдетальной мелкой сетке (Zero-Shot Super-Resolution), и алгоритм мгновенно выдаст физически корректный результат. Эта технология уже применяется для сверхбыстрого предсказания свойств многофазных жидкостей в пористых средах (добыча нефти) и прогнозирования погоды.

Подробнее

Решение уравнений Навье-Стокса: алгоритмы SIMPLE и PISO для несжимаемых жидкостей

Вычислительный парадокс давления в несжимаемых потоках

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

В сжимаемой сверхзвуковой газовой динамике давление жестко связано с текущей плотностью и температурой через уравнение состояния, а сама плотность вычисляется из эволюционного уравнения неразрывности (сохранения массы). Для воды или медленного воздуха уравнение неразрывности превращается в жесткое геометрическое кинематическое условие: математическая дивергенция вектора скорости должна быть строго равна абсолютному нулю в каждой физической точке пространства (сколько массы втекло в контрольный объем, ровно столько же и вытекло). Эволюционного дифференциального уравнения по времени для вычисления давления просто не существует в природе! Давление здесь выступает абстрактным математическим множителем Лагранжа — мгновенной фиктивной силой, которая невидимо распределяется по всему гидродинамическому объему так, чтобы итоговое поле скоростей всегда и везде удовлетворяло условию строгой бездивергентности.

Связка давления и скорости: алгоритм предиктор-корректор SIMPLE

Для блестящего решения этой нетривиальной проблемы связи давления и скорости (pressure-velocity coupling) в 1972 году физики Сполдинг и Патанкар разработали легендарный алгоритм SIMPLE (Semi-Implicit Method for Pressure Linked Equations). Идея этого численного алгоритма концептуально основана на классическом методе предиктор-корректор.

На самом первом этапе вычислений (расчетный шаг предиктора) алгоритм берет какое-то начальное, грубо угаданное поле давления. Используя это физически неточное давление, решатель интегрирует уравнения импульса и находит промежуточное (рабочее) поле скоростей. Естественно, это сырое поле скоростей будет «неправильным» — оно грубо нарушит физический закон сохранения массы (дивергенция не будет равна нулю, жидкость будет «возникать» из ниоткуда). На втором этапе (шаг математического корректора) алгоритм составляет специальное уравнение для поправки давления (Уравнение Пуассона). Суть его в том, что вычисленная поправка давления должна создать такие микроскопические дополнительные векторы скорости, которые в сумме с промежуточными скоростями строго и математически точно обнулят дивергенцию! Вычислив эту сложную поправку, алгоритм корректирует и поле давлений, и поле скоростей. Поскольку система гидродинамики сильно нелинейна, этот вычислительный цикл приходится повторять десятки итераций на каждом микрошаге по времени. Дальнейшим мощным математическим развитием этой концепции стал алгоритм PISO, обладающий дополнительным вторым этапом коррекции и ставший мировым стандартом в открытых CFD-кодах.

Подробнее

Методы доверительных областей (Trust Region) в нелинейной оптимизации

Мощная альтернатива классическому линейному поиску

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

Методы доверительных областей (Trust Region Methods) используют совершенно противоположную вычислительную философию. Вместо того чтобы сначала угадать направление, а потом мучительно искать длину геометрического шага, этот мощный алгоритм сначала определяет максимальную допустимую длину шага (сферический радиус доверия), а уже строго внутри этой очерченной пространственной области ищет оптимальное направление и точку минимума.

Построение суррогатной квадратичной модели и адаптивный радиус

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

Затем численный алгоритм ищет строгий локальный минимум именно этой простой суррогатной квадратичной модели, но с очень жестким математическим ограничением — итоговое решение ни при каких условиях не должно выходить за пределы вычисленного радиуса доверительной гиперсферы. Найдя перспективную пробную точку, алгоритм проверяет ее на реальной (исходной сложнейшей) целевой функции. Вычисляется важный коэффициент согласования: точное отношение реального падения функции к теоретически предсказанному падению на параболе. Если предсказание совпало с суровой реальностью (коэффициент оказался близок к 1), значит наша квадратичная модель отлично описывает физику! Алгоритм с радостью принимает шаг и на следующей итерации оптимистично увеличивает радиус доверия. Если же реальное значение функции оказалось значительно хуже предсказанного (парабола «обманула» расчеты из-за сверхвысокой нелинейности), шаг немедленно отвергается, точка не меняется, а радиус доверительной области принудительно и жестко сжимается. Этот класс методов обладает потрясающей робастностью сходимости.

Подробнее