Main menu

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

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

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

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

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

Чтобы обуздать экспоненциальный взрыв и восстановить вычислительную стабильность, математиками Келлером, Осборном и Деуфлхардом был разработан метод многократной стрельбы (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), значит наша квадратичная модель отлично описывает физику! Алгоритм с радостью принимает шаг и на следующей итерации оптимистично увеличивает радиус доверия. Если же реальное значение функции оказалось значительно хуже предсказанного (парабола «обманула» расчеты из-за сверхвысокой нелинейности), шаг немедленно отвергается, точка не меняется, а радиус доверительной области принудительно и жестко сжимается. Этот класс методов обладает потрясающей робастностью сходимости.

Подробнее

Численная оптимизация топологии: как математика создает бионический дизайн

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

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

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

Метод SIMP: плотность материала как непрерывная переменная

Самым распространенным и вычислительно эффективным алгоритмом топологической оптимизации в мировой индустрии является метод SIMP (Solid Isotropic Material with Penalization). Алгоритм начинает работу с того, что разбивает всю доступную расчетную область на мелкие кубические конечные элементы (воксели). Изначально каждому такому элементу программно присваивается математическая «плотность» от 0 (полная воздушная пустота) до 1 (цельный монолитный металл). Проблема заключается в том, что в реальном физическом мире материал детали не может быть полупрозрачным, с плотностью 0.5. Нам на выходе нужен строгий дискретный ответ: либо там есть сталь, либо там воздух.

Чтобы заставить математический алгоритм стремиться к четким черным и белым границам, метод SIMP вводит штрафование (Penalization) промежуточных серых плотностей. Жесткость каждого элемента привязывается к его расчетной плотности через нелинейную степенную функцию (обычно плотность возводится в куб). Таким образом, полупустой элемент со значением 0.5 будет иметь всего 12.5% от максимальной жесткости, но «весить» как половина полного. Оптимизатору (обычно используется метод MMA - Method of Moving Asymptotes) становится математически крайне невыгодно и энергетически дорого оставлять такие серые зоны. Он агрессивно выталкивает плотность каждого элемента к строгим краям: либо в 0, либо в 1. На каждом шаге итерации алгоритм решает гигантскую СЛАУ упругости, вычисляет градиенты функции чувствительности и плавно обновляет плотности элементов. С развитием промышленной 3D-печати этот метод стал абсолютным стандартом в аэрокосмической отрасли.

Подробнее

Метод граничных элементов (МГЭ): снижение размерности задачи на единицу

Снижение размерности: от объема к поверхности

Метод конечных элементов (МКЭ) и метод конечных разностей (МКР) безусловно доминируют в современных инженерных расчетах. Однако оба эти метода являются полностью объемными: они требуют построения густой расчетной сетки во всей 3D-области математического моделирования. Если вам нужно рассчитать обтекание подводной лодки в открытом океане или напряжение вокруг микротрещины в огромном стальном блоке, объемная сетка съест колоссальное количество памяти суперкомпьютера. В таких специфических случаях на помощь приходит элегантный математический аппарат — Метод граничных элементов (МГЭ, или Boundary Element Method, BEM).

Основная геометрическая идея МГЭ заключается в снижении размерности задачи ровно на одну физическую единицу. Для массивного трехмерного тела сетка строится только на его внешней двумерной поверхности (оболочке). Для плоского двумерного сечения крыла самолета сетка представляет собой просто одномерную кривую, очерчивающую внешний контур. Это радикально, на несколько порядков, уменьшает количество узлов и неизвестных переменных, позволяя решать задачи высочайшей сложности на обычных ПК. Но возникает парадокс: как можно математически рассчитать то, что происходит внутри объема, зная только числовые данные на его поверхности?

Интегральные уравнения и фундаментальные решения

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

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

Подробнее

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

От классических полиномов к глубокому обучению

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

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

Физически информированные нейросети (PINNs)

Настоящей революцией в численном решении дифференциальных уравнений в частных производных (УЧП) стало появление физически информированных нейронных сетей (Physics-Informed Neural Networks, PINNs). Идея PINN заключается во внедрении физических законов (самих дифференциальных уравнений) непосредственно в функцию потерь (Loss function) нейронной сети во время ее математического обучения.

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

Подробнее
Subscribe to this RSS feed

Соц. сети