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

Метод Галеркина: проекционный подход к решению краевых задач

Переход от дифференциальных уравнений к вариационным задачами

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

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

Пробные и тестовые функции: слабая формулировка

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

Затем исходное дифференциальное уравнение умножается на весовую (или тестовую) функцию и интегрируется по всей расчетной области. В классическом методе Бубнова-Галеркина в качестве тестовых функций выбираются те же самые функции, что и в базисе пробных функций. Приравнивая полученный интеграл от невязки к нулю, мы гарантируем, что ошибка нашего приближения не содержит проекций на выбранный базис. Интегрирование по частям (формула Грина) позволяет перебросить часть производных с неизвестной функции на гладкую тестовую функцию. Это приводит к так называемой «слабой формулировке» задачи: требования к гладкости искомого решения радикально снижаются. Уравнение второго порядка теперь можно решать, используя базисные функции, имеющие лишь первые производные.

Сведение к СЛАУ и сходимость

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

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

Подробнее

Квази-Монте-Карло: последовательности с низким расхождением и Соболь

Проблемы псевдослучайных чисел в методе Монте-Карло

Классический метод Монте-Карло (МК), как мы выяснили ранее, является незаменимым инструментом для вычисления интегралов в пространствах высокой размерности (от 5 до 100 и более измерений), где традиционные сеточные методы бессильны из-за «проклятия размерности». Однако стандартный МК имеет серьезный недостаток — очень медленную сходимость, пропорциональную O(1/sqrt(N)). Чтобы увеличить точность результата всего на один знак (в 10 раз), нужно сгенерировать в 100 раз больше случайных точек.

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

Последовательности с низким расхождением

Чтобы решить эту проблему, математики разработали метод Квази-Монте-Карло (КМК). Ключевое отличие заключается в том, что КМК полностью отказывается от случайности! Вместо псевдослучайных генераторов используются жестко детерминированные алгоритмы, которые генерируют так называемые последовательности с низким расхождением (Low-discrepancy sequences). Эти последовательности заполняют многомерный гиперкуб максимально равномерно, избегая образования как кластеров, так и пустых дыр. Точки расставляются так, чтобы каждое новое значение ложилось точно в самый большой незаполненный пробел.

С точки зрения математики, расхождение (discrepancy) — это мера того, насколько плотность точек в произвольном параллелепипеде отклоняется от идеальной равномерной плотности. Для случайных чисел расхождение велико. Для сеток оно мало, но их невозможно построить в высоких размерностях. Последовательности КМК (такие как последовательности Холтона, Фора, Нидеррайтера) дают идеальный компромисс. Самой популярной и эффективной в индустрии является ЛП-тау последовательность Соболя, разработанная выдающимся советским математиком Ильей Мееровичем Соболем в 1967 году.

Применение КМК на финансовых рынках и в графике

Замена генератора случайных чисел на последовательность Соболя превращает классический Монте-Карло в Квази-Монте-Карло, при этом скорость сходимости меняется радикально: она приближается к O(1/N). Это означает, что для достижения той же точности требуется не в 100 раз больше вычислений, а всего лишь в 10 раз. Это колоссальный скачок в эффективности.

Метод КМК произвел настоящую революцию в инвестиционных банках. При ценообразовании сложных финансовых производных (например, азиатских или ипотечных опционов), где выплаты зависят от траектории цены актива за каждый из 365 дней в году, возникает интеграл 365-й размерности. Применение КМК позволило оценивать такие портфели за секунды вместо часов. Также последовательности Соболя стали индустриальным стандартом в компьютерной графике и рендеринге (алгоритмы Path Tracing). Фотореалистичное освещение в современных голливудских фильмах (Pixar, Disney) рассчитывается путем трассировки лучей, где направления лучей выбираются именно на основе последовательностей КМК, обеспечивая гладкую картинку без цифрового шума при минимальном количестве лучей на пиксель.

Подробнее

Численное конформное отображение: геометрия на службе аэродинамики

Магия комплексного анализа в геометрии

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

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

Преобразование Жуковского и профили крыла

В начале 20-го века Николай Егорович Жуковский использовал простейшую функцию комплексного переменного w = z + 1/z для создания первого теоретически рассчитанного профиля крыла. Эта аналитическая формула (преобразование Жуковского) позволяет отобразить окружность из плоскости z на каплевидный профиль крыла в плоскости w. Рассчитав поле скоростей и давлений вокруг цилиндра (что делается легко), можно перенести эти данные обратно на профиль крыла и получить подъемную силу.

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

Интеграл Кристоффеля-Шварца для многоугольников

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

Для поиска этих параметров используются специализированные численные методы многомерной оптимизации или алгоритмы продолжения по параметру. Существуют мощные программные пакеты (например, SC Toolbox для MATLAB, разработанный Дрисколлом и Трефетеном), которые автоматически рассчитывают эти параметры и строят конформные сетки. Численные конформные отображения активно применяются не только в расчете подъемной силы, но и в проектировании СВЧ-волноводов, где нужно точно рассчитать емкость и индуктивность кабелей сложного сечения.

Подробнее

Методы Рунге-Кутты-Нистрёма для дифференциальных уравнений второго порядка

Специфика систем второго порядка в механике

Фундаментальный закон классической механики (второй закон Ньютона) утверждает, что сила равна массе, умноженной на ускорение (которое является второй производной от координаты по времени). Поэтому подавляющее большинство задач кинематики, баллистики, теории колебаний и небесной механики изначально формулируются в виде обыкновенных дифференциальных уравнений (ОДУ) второго порядка: y'' = f(t, y, y').

Стандартный подход к численному решению таких уравнений, которому учат в базовых курсах, заключается в их понижении до систем уравнений первого порядка. Мы просто вводим новую переменную (скорость v = y') и получаем систему: y' = v и v' = f(t, y, v). Эту систему затем решают классическим методом Рунге-Кутты (например, RK4). Однако этот универсальный подход далеко не всегда является оптимальным. Во-первых, он искусственно увеличивает размерность задачи в два раза. Во-вторых, он не учитывает внутреннюю структуру исходного уравнения второго порядка, что может приводить к потере точности при вычислении координат.

Специализированные методы: алгоритмы Нистрёма

В 1925 году финский математик Эверт Йоханнес Нистрём разработал специализированное семейство методов Рунге-Кутты, предназначенное для прямого интегрирования дифференциальных уравнений второго порядка без предварительного сведения их к системам первого порядка. Эти методы получили название методов Рунге-Кутты-Нистрёма (RKN).

Особую ценность представляют алгоритмы RKN для специального класса уравнений вида y'' = f(t, y), в которых правая часть (сила) не зависит от первой производной (скорости). Такие задачи описывают консервативные механические системы (без трения и диссипации энергии), например, движение планет вокруг Солнца или колебания идеального маятника. Применяя методы Нистрёма к таким системам, можно добиться феноменальной эффективности: они требуют значительно меньше вычислений правой части (вызовов функции силы) на каждом шаге интегрирования по сравнению со стандартным методом Рунге-Кутты того же порядка точности, сохраняя при этом идеальную гладкость траектории.

Симплектические интеграторы и сохранение энергии

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

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

Подробнее

Быстрое преобразование мультиполей (FMM): прорыв в задачах N тел

Гравитация, кулоновские силы и проблема квадратичной сложности

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

Если мы попытаемся вычислить силы «в лоб» прямым суммированием, нам придется рассчитать взаимодействие каждой частицы с каждой другой. Это означает, что вычислительная сложность алгоритма растет пропорционально квадрату числа частиц (O(N^2)). Для миллиона частиц потребуется триллион операций на каждом шаге интегрирования по времени. Это абсолютный вычислительный тупик: даже самые мощные суперкомпьютеры не способны моделировать динамику плазмы или формирование галактик в течение приемлемого времени, используя прямые методы расчета.

Идея кластеризации и мультипольное разложение

В 1987 году Владимир Рохлин и Лесли Грингард опубликовали статью, описывающую Быстрое преобразование мультиполей (Fast Multipole Method, FMM). Этот алгоритм был признан журналом Computing in Science & Engineering одним из 10 главных алгоритмов 20 века. FMM радикально снижает сложность задачи N тел с O(N^2) до O(N log N), а в некоторых реализациях — и до строгого O(N).

Идея алгоритма базируется на простом физическом наблюдении: если вы находитесь на Земле, вам не нужно вычислять гравитационное притяжение от каждого отдельного камня на Луне. Вы можете рассматривать Луну как единый точечный объект с суммарной массой, расположенный в ее центре масс. FMM применяет этот принцип иерархически. Пространство разбивается на кубическое дерево (octree в 3D или quadtree в 2D). Частицы, находящиеся близко друг к другу (в соседних ячейках), взаимодействуют напрямую. Но для частиц в далеких кластерах их суммарное гравитационное или электрическое поле заменяется компактным математическим рядом — мультипольным разложением.

Математика сферических гармоник

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

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

Подробнее

Многосеточные методы (Multigrid): масштабируемое решение сверхбольших СЛАУ

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

При решении эллиптических дифференциальных уравнений в частных производных (например, уравнения Пуассона для электростатического потенциала или давления в гидродинамике) возникают огромные системы линейных алгебраических уравнений. Мы уже упоминали методы подпространств Крылова, которые отлично справляются с этой задачей. Однако существует класс алгоритмов, который математически превосходит их по масштабируемости, достигая теоретического предела вычислительной эффективности — линейной сложности O(N). Это многосеточные методы (Multigrid methods, MG).

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

Переход между масштабами: сужение и продолжение

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

Процесс начинается на мелкой сетке. Метод Гаусса-Зейделя делает несколько итераций, быстро сглаживая ошибку (высокочастотный шум). Оставшаяся плавная ошибка (невязка) переносится на более грубую сетку с помощью математического оператора сужения (restriction). На грубой сетке эта ошибка кажется «острой», и сглаживатель снова эффективно ее давит. Этот процесс погружения продолжается вплоть до самой грубой сетки (состоящей из десятка узлов), где СЛАУ решается точно прямым методом практически мгновенно. Затем найденная поправка возвращается обратно на мелкие сетки через оператор интерполяции или продолжения (prolongation), корректируя решение на каждом уровне.

V-циклы, W-циклы и алгебраический многосеточный метод (AMG)

Последовательность спусков на грубые сетки и подъемов на мелкие формирует так называемые циклы. Самый популярный из них — V-цикл (вниз и вверх). Существуют также W-циклы, которые проводят больше времени на грубых уровнях для лучшей очистки от низкочастотных ошибок, и F-циклы (Full Multigrid), которые вообще начинают работу с самой грубой сетки, постепенно переходя к точным вычислениям на мелких.

Классический геометрический многосеточный метод требует явного знания геометрии расчетной области для построения грубых сеток, что трудно реализовать для деталей со сложной криволинейной границей. Поэтому в современных промышленных пакетах (ANSYS, COMSOL) применяется Алгебраический Многосеточный Метод (AMG). Он строит иерархию сеток и операторы перехода не на основе геометрических координат узлов, а исключительно на основе значений коэффициентов в самой матрице СЛАУ. AMG работает как мощный черный ящик (или предобуславливатель для метода сопряженных градиентов), обеспечивая сверхбыстрое решение систем с миллиардами уравнений на суперкомпьютерах.

Подробнее

Численное обращение преобразования Лапласа: борьба с некорректными задачами

Могущественный инструмент операционного исчисления

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

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

Некорректность по Адамару и взрыв ошибок

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

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

Алгоритм Гавера-Стефеста и методы Фурье

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

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

Подробнее

Интервальный анализ: вычисления с гарантированной математической точностью

Катастрофическая потеря значимости чисел с плавающей запятой

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

Для управления полетами космических аппаратов, расчета доз облучения в медицинской лучевой терапии и проверки безопасности атомных реакторов такой риск неприемлем. Нам нужен способ не просто получить приближенный ответ, но и получить строгие математические гарантии того, что истинный ответ гарантированно находится внутри определенного диапазона. Эту революционную идею в 1960-х годах воплотил в жизнь американский математик Рамон Мур, создав новую ветвь математики — интервальный анализ.

Замена чисел на множества

В интервальном анализе базовым объектом оперирования является не конкретное число (скаляр), а замкнутый интервал (отрезок) на числовой оси: [a, b]. В этом отрезке a — это строгая нижняя граница, а b — строгая верхняя. Любые исходные данные с погрешностью (например, длина детали 10 мм с погрешностью 0.1 мм) сразу представляются как интервал [9.9, 10.1].

Далее переопределяются все арифметические операции (сложение, вычитание, умножение, деление) и базовые функции (синус, экспонента). Когда компьютер складывает два интервала, он складывает их нижние и верхние границы с особым аппаратным округлением (нижняя граница всегда округляется в меньшую сторону, верхняя — в большую). Это гарантирует, что истинный ответ никогда не выскользнет за пределы вычисленного результирующего интервала. Иными словами, интервальные вычисления предоставляют стопроцентную, математически строгую оценку глобальной погрешности алгоритма. Если в результате расчета орбиты спутника мы получили интервал расстояния до Земли [400 км, 405 км], мы можем быть абсолютно уверены, что спутник не упадет.

Проблема зависимости (эффект переоценки)

Главным препятствием на пути повсеместного внедрения интервальной арифметики является так называемая «проблема зависимости». Если в классической математике x - x всегда равно 0, то в интервальной арифметике, если мы возьмем интервал X = [1, 2] и вычтем его сам из себя, мы получим интервал [-1, 1], а не ноль. Алгоритм «забывает», что интервалы физически связаны (зависят от одной переменной), и рассматривает их как два независимых множества, принимающих наихудшие возможные значения.

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

Подробнее

Численные методы в финансовой математике: опционы и модель Блэка-Шоулза

От законов теплопроводности к ценообразованию на Уолл-Стрит

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

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

Биномиальные деревья и сеточные методы

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

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

Подробнее

Решение жестких краевых задач: метод ортогональной прогонки Годунова

Вычислительная катастрофа при интегрировании

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

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

Спасение через ортогонализацию

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

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

Подробнее