Высшая математика

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

Альтернатива Галеркину: сила точечного совпадения

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

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

Базис B-сплайнов: локальность и высокая гладкость

Теоретически в качестве пробной функции можно взять один длинный глобальный полином (многочлен). Но, как мы помним из феномена Рунге, глобальные многочлены высоких степеней сильно осциллируют и ведут себя нестабильно. Идеальным базисом для метода коллокации оказались B-сплайны (базисные сплайны), разработанные Исааком Шёнбергом и Карлом де Бором.

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

Оптимальные узлы: корни полиномов Чебышева и Гаусса

Самым критическим моментом в методе сплайн-коллокации является выбор расположения самих точек коллокации. Если расставить эти точки равномерно, метод может потерять устойчивость или сходиться очень медленно. Математический анализ, проведенный Карлом де Бором и Ричардом Сварцем в 1973 году, доказал потрясающую вещь: точность метода кардинально возрастает (явление суперсходимости), если в качестве точек коллокации внутри каждого сеточного интервала использовать корни ортогональных полиномов (узлы Гаусса-Лежандра или Чебышева).

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

Подробнее

Численное интегрирование функций с сингулярностями: метод Tanh-Sinh

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

Классические методы численного интегрирования (квадратуры Ньютона-Котеса, формулы Симпсона или Гаусса) блестяще справляются с гладкими функциями, которые заданы на конечном отрезке. Однако в физических задачах математикам постоянно приходится сталкиваться с интегралами, которые имеют математические особенности — сингулярности. Например, при расчетах потенциала заряженной нити или в задачах механики разрушения (напряжения на кончике трещины) подынтегральная функция может уходить в бесконечность на краях отрезка или иметь вертикальные асимптоты. Классический алгоритм, столкнувшись с бесконечностью (делением на ноль), выдаст ошибку Not-a-Number (NaN) или будет сходиться катастрофически медленно.

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

Магия нелинейных замен переменных

Второй, гораздо более мощный алгоритмический подход основан на нелинейных заменах координат. Идея заключается в том, чтобы так исказить ось интегрирования, чтобы «растянуть» проблемные участки (около сингулярности) и свести бесконечности к нулю. Одной из первых успешных замен такого рода стала подстановка Игараси-Мори-Такахаси (IMT).

Но настоящим прорывом и шедевром японской вычислительной школы стало создание в 1974 году квадратурной формулы гиперболического тангенса-синуса (Tanh-Sinh quadrature, или метод двойной экспоненциальной интеграции). Этот метод был разработан Хидетаси Такахаси и Масатаке Мори. Идея заключается во введении сложнейшей подстановки: x = tanh(pi/2 * sinh(t)). Эта функция обладает потрясающим топологическим свойством: она отображает бесконечный интервал (от минус до плюс бесконечности) строго в конечный отрезок [-1, 1].

Двойное экспоненциальное затухание

При выполнении этой замены дифференциал интегрирования (якобиан) порождает функцию, которая убывает к нулю с феноменальной скоростью — по закону двойной экспоненты: exp(-C * exp(|t|)). Эта скорость убывания настолько сокрушительна, что она способна мгновенно подавить (помножить на ноль) абсолютно любую математическую сингулярность исходной функции на краях отрезка!

После такой замены переменной мы получаем новый, гладкий и колоколообразный интеграл на бесконечной прямой, для вычисления которого применяется самое простое правило трапеций (да, базовая формула трапеций на бесконечном интервале для гладких убывающих функций работает точнее сложных методов Гаусса). Алгоритм Tanh-Sinh позволяет с феноменальной легкостью брать интегралы с любыми концевыми сингулярностями, логарифмическими разрывами или бесконечными производными, достигая невероятной точности в тысячи знаков после запятой в системах произвольной разрядности. Сегодня этот алгоритм встроен в ядра таких мощных математических систем, как MPMath в Python и Boost.Math в C++.

Подробнее

Алгоритмы на графах для разреженных матриц: Катхилл-Макки и минимальная степень

Скрытая угроза прямых решателей: проблема заполнения

Мы неоднократно обращали внимание на то, что огромные системы уравнений (СЛАУ), возникающие в инженерных расчетах, являются сверхразреженными: в них присутствуют миллионы нулей и лишь малая доля значащих чисел. Для их хранения созданы эффективные форматы вроде CSR. Однако если мы попытаемся решить такую систему точным, прямым методом (LU-разложение или разложение Холецкого для симметричных матриц), мы столкнемся с катастрофическим явлением, которое в вычислительной алгебре называется «заполнением» (fill-in).

В процессе прямого исключения переменных (вспомните метод Гаусса), вычитая одну строку из другой, алгоритм неизбежно складывает ненулевые элементы. Если в исходной строке стоял спасительный ноль, то в процессе прямого хода на его месте может появиться значащее число. К концу алгоритма матрица, которая изначально была почти пустой, рискует превратиться в гигантскую, абсолютно плотную структуру. Миллионы нулей исчезнут, алгоритм сожрет всю оперативную память кластера (Out of Memory) и процесс остановится. Как же предотвратить это катастрофическое разрастание данных?

Матрица как топологический граф

Гениальное решение пришло на стыке алгебры и дискретной математики. Любую симметричную матрицу можно представить в виде неориентированного графа. Каждая строка (уравнение) — это узел графа (вершина). Если матричный элемент на пересечении строки A и столбца B не равен нулю, значит между узлом A и узлом B в графе существует физическое ребро связи. Заполнение матрицы в процессе LU-разложения в терминах теории графов означает добавление новых, фиктивных ребер между узлами.

Порядок, в котором мы перечисляем уравнения (нумеруем переменные x1, x2, x3), не меняет физического решения системы, но он критически, на порядки, влияет на уровень заполнения матрицы нулями в памяти! Следовательно, перед тем как запустить вычислительно тяжелое LU-разложение, нам нужно запустить быстрый алгоритм переупорядочивания (эвристику на графе), который найдет оптимальную нумерацию строк и столбцов матрицы для минимизации будущего заполнения.

Ширина ленты и алгоритмы Reverse Cuthill-McKee (RCM)

Первым классом таких алгоритмов стали методы уменьшения ширины ленты (Bandwidth reduction). Исторически самым известным из них является алгоритм Катхилла-Макки (Cuthill-McKee), а точнее его улучшенная версия — обратный алгоритм Катхилла-Макки (RCM). Алгоритм начинает поиск с периферийного узла графа (узла с минимальным количеством связей) и послойно, методом поиска в ширину (BFS), перенумеровывает все связанные с ним узлы. В результате все ненулевые элементы матрицы «сбегаются» и плотно прижимаются к главной диагонали, образуя узкую ленту. Вне этой ленты остаются чистые нули, которые гарантированно никогда не заполнятся при разложении.

Для более сложных и неструктурированных 3D-сеток ленточные алгоритмы уступают место семейству методов минимальной степени (Minimum Degree Algorithm — AMD). Алгоритм AMD работает на опережение: на каждом шаге он симулирует процесс исключения Гаусса и всегда выбирает для исключения ту переменную, узел которой в графе в данный момент имеет наименьшее число соседей (минимальную степень). Это локально жадная эвристика, но она минимизирует создание новых ребер на каждом микрошаге. Сегодня мощнейшие программные комплексы на графах (такие как METIS или Scotch) выполняют глубокое иерархическое разбиение графа матрицы на кластеры, что является обязательным «нулевым» шагом перед любым суперкомпьютерным матричным расчетом.

Подробнее

Метод спектральных элементов: объединение гибкости МКЭ и точности рядов

Объединение гибкости и невероятной точности

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

В 1984 году Энтони Патера предложил элегантный математический мост между этими двумя мирами, создав Метод спектральных элементов (Spectral Element Method, SEM). Идея метода заключается в том, чтобы взять лучшее от обоих подходов. Расчетная область со сложной геометрией разбивается на крупные макроэлементы (гексаэдры в 3D или четырехугольники в 2D), подобно тому, как это делается в МКЭ. Эта операция обеспечивает необходимую геометрическую гибкость алгоритма и возможность распараллеливания.

Внутри элемента: полиномы высоких порядков

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

Узлы интерполяции внутри макроэлемента располагаются не равномерно, а сгущаются к границам элемента (так называемые узлы Гаусса-Лобатто-Лежандра). Такое неравномерное распределение математически необходимо для предотвращения катастрофических осцилляций на краях (феномена Рунге) при использовании полиномов 10-й, 15-й или даже 20-й степени. Матрицы масс и жесткости, возникающие при дискретизации по методу Галеркина на таких узлах, вычисляются с использованием спектрально-точных квадратурных формул Гаусса. Благодаря этому массовые матрицы получаются строго диагональными, что позволяет обращать их мгновенно, делая метод невероятно эффективным для нестационарных задач.

Применение в гидродинамике и глобальной сейсмологии

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

Первая область — это прямое численное моделирование (DNS) турбулентных течений несжимаемой жидкости при низких и средних числах Рейнольдса. Программный комплекс Nek5000, разработанный Полом Фишером в Аргоннской национальной лаборатории, использует метод спектральных элементов для расчета гидродинамики ядерных реакторов на крупнейших мировых суперкомпьютерах. Вторая область — это глобальная сейсмология. Распространение сейсмических волн от землетрясений сквозь мантию Земли идеально описывается упругими волновыми уравнениями. Комплекс SPECFEM3D использует метод спектральных элементов для создания высокоточных томографических моделей строения нашей планеты, так как алгоритм позволяет избежать численной дисперсии (искусственного искажения скорости волн), от которой сильно страдают обычные конечно-разностные схемы.

Подробнее

Уравнения Гамильтона-Якоби-Беллмана в численном оптимальном управлении

Поиск идеальной траектории в пространстве состояний

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

Существуют два главных математических пути решения таких задач. Первый подход — это принцип максимума Понтрягина, который приводит к сложнейшим краевым задачам для обыкновенных дифференциальных уравнений (ОДУ). Второй, не менее мощный подход — метод динамического программирования, разработанный американским математиком Ричардом Беллманом в 1950-х годах. Принцип оптимальности Беллмана философски гласит: «Оптимальное поведение обладает тем свойством, что каковы бы ни были начальное состояние и решение в начальный момент, последующие решения должны составлять оптимальное поведение относительно состояния, получающегося в результате первого решения». В непрерывном времени этот принцип приводит к знаменитому уравнению Гамильтона-Якоби-Беллмана (HJB).

Сложность УЧП и вязкостные решения

Уравнение HJB — это нелинейное дифференциальное уравнение в частных производных (УЧП) первого или второго порядка (если в системе есть стохастический шум). Его неизвестной является функция ценности (Value function) — минимальная ожидаемая стоимость достижения цели из любого текущего состояния. Решив уравнение HJB, мы получаем закон управления в виде обратной связи (closed-loop control), что критически важно для робастности реальных систем управления, подверженных помехам.

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

Численные схемы: метод Годунова и Upwind

Численное решение УЧП Гамильтона-Якоби-Беллмана является экстремально сложной задачей вычислительной математики. Простые конечно-разностные схемы (например, метод центральных разностей) здесь абсолютно не работают — они генерируют дикие нефизичные осцилляции в местах изломов функции ценности.

Для стабильного решения необходимо применять специализированные односторонние разностные схемы (Upwind schemes) первого порядка. Идея, заимствованная из вычислительной газовой динамики (схемы Годунова и Лакса-Фридрихса), заключается в том, что разностная производная должна браться строго против направления потока информации. Кроме того, уравнение HJB требует решения задачи минимизации (вычисления Гамильтониана) на каждом узле сетки, что делает алгоритм крайне ресурсоемким. Но главная преграда — это «проклятие размерности» Беллмана. Размерность сетки для УЧП равна количеству переменных состояния системы. Для 6-мерного дрона-квадрокоптера построение полной сетки невозможно ни на одном суперкомпьютере мира. В последние годы для пробивания этого барьера начали активно применяться методы машинного обучения, где уравнения HJB решаются с использованием глубоких нейронных сетей (Physics-Informed Neural Networks), открывая новую эру в робототехнике и теории управления.

Подробнее

Методы решения интегральных уравнений Вольтерры: системы с эффектом затухающей памяти

Когда прошлое непрерывно влияет на будущее

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

Такие физические процессы с «исторической памятью» адекватно описываются интегральными уравнениями Вольтерры. Главной отличительной чертой уравнения Вольтерры (в отличие от уравнений Фредгольма) является то, что верхний предел интеграла не является фиксированным числом, а представляет собой саму независимую переменную (чаще всего время t). То есть мы интегрируем искомую функцию от момента зарождения системы (t=0) до текущего момента времени. Искомая функция может стоять только под знаком интеграла (уравнения первого рода) или присутствовать и вне его (уравнения второго рода).

Численные квадратуры и сведение к рекуррентной системе

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

Идея метода квадратур интуитивно понятна. Мы вводим на оси времени сетку дискретных узлов. Затем мы заменяем непрерывный интеграл из уравнения Вольтерры его конечно-разностным приближением (например, по формуле трапеций или методу Симпсона). Благодаря тому, что верхний предел интеграла переменный, мы получаем не единую гигантскую матрицу (как у Фредгольма), а удобную треугольную систему линейных алгебраических уравнений (СЛАУ).

Эту систему не нужно решать сложными матричными методами. Она решается последовательно, шаг за шагом (методом прямой подстановки). Чтобы вычислить значение функции на n-ном шаге, алгоритму достаточно подставить в формулу уже известные, вычисленные значения функции на всех предыдущих шагах. Этот процесс очень напоминает численное интегрирование ОДУ, но с существенной оговоркой: объем вычислений растет катастрофически. На каждом новом временном шаге алгоритм должен заново пересчитать интеграл по всей предыстории.

Методы коллокации и аналоги Рунге-Кутты

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

Для уравнений Вольтерры второго рода математиками были разработаны специальные адаптации методов Рунге-Кутты (V-R-K методы). Они обеспечивают высокий порядок сходимости, но анализ их численной устойчивости невероятно сложен. В отличие от ОДУ, устойчивость схемы для интегрального уравнения зависит не только от шага по времени, но и от поведения ядра интегрального оператора. Ошибки округления могут накапливаться в интегральной сумме и приводить к экспоненциальной раскачке алгоритма, что требует применения жестко-устойчивых неявных квадратурных формул дифференцирования назад (BDF).

Подробнее

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

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

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

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

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

В 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 независимых коротких выстрелов, стартующих одновременно из разных точек. Поскольку каждый интервал короткий, быстрорастущие решения просто не успевают раздуться до машинной бесконечности, и интегрирование каждого кусочка проходит успешно и стабильно.

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

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

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

Подробнее