Поиск идеального численного алгоритма

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

Для формального анализа устойчивости в вычислительной математике используется простое линейное тестовое уравнение Далквиста: y' = lambda * y, где lambda — комплексное число с отрицательной вещественной частью. Аналитическое решение этого уравнения — затухающая экспонента, стремящаяся к нулю. Логично требовать, чтобы и численное решение, выдаваемое компьютером при любом (даже большом) шаге интегрирования h, также убывало и не уходило в бесконечность. Алгоритмы, обладающие таким свойством для любых lambda в левой полуплоскости, называются А-устойчивыми (абсолютно устойчивыми). Это свойство является Святым Граалем для разработчиков численных методов.

Барьерные теоремы Далквиста: жестокая математическая реальность

В 1956 и 1963 годах шведский математик Жермен Далквист опубликовал серию теорем, которые потрясли мир вычислительной математики. Они установили непреодолимые теоретические пределы для целого класса популярных алгоритмов — линейных многошаговых методов (ЛММ), к которым относятся методы Адамса и формулы дифференцирования назад (BDF).

Первый барьер Далквиста доказал, что порядок точности любого устойчивого явного линейного многошагового метода (k-шагового) не может превышать k при нечетном k, и k+2 при четном. Это означает, что мы не можем бесконечно повышать точность алгоритма, не теряя его устойчивости.

Второй барьер Далквиста оказался еще более суровым. Он математически строго доказал три вещи: 1) Ни один явный линейный многошаговый метод не может быть А-устойчивым. (Поэтому явные методы бесполезны для жестких систем). 2) Порядок точности любого А-устойчивого неявного линейного многошагового метода не может быть больше 2. 3) Среди всех А-устойчивых методов второго порядка наименьшую ошибку усечения дает классический неявный метод трапеций.

L-устойчивость и пути обхода барьеров

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

Решение было найдено в двух направлениях. Первое — это смягчение требований. Было введено понятие «жесткой устойчивости» (stiff stability). Формулы дифференцирования назад (BDF) до 6-го порядка не являются А-устойчивыми (вблизи мнимой оси у них есть небольшие зоны неустойчивости), но они жестко устойчивы и прекрасно работают на практике. Второе направление — выход за рамки линейных многошаговых методов. Неявные методы Рунге-Кутты (такие как методы Радо или Гаусса-Лежандра) не подчиняются теоремам Далквиста. Они могут быть А-устойчивыми и при этом иметь любой, сколь угодно высокий порядок точности (например, метод Гаусса с s стадиями имеет порядок 2s). Кроме того, для борьбы с бесконечно быстрыми процессами (когда действительная часть lambda стремится к минус бесконечности) было введено еще более строгое свойство — L-устойчивость. L-устойчивые алгоритмы (например, неявный метод Эйлера) мгновенно обнуляют высокочастотные компоненты, что делает их идеальными для моделирования ударных волн и резких переключателей в транзисторах.

Обратные задачи и некорректность по Адамару

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

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

Регуляризация: искусство разумного компромисса

Долгое время считалось, что некорректные задачи вообще не имеют физического смысла и не подлежат численному решению. Революция произошла в 1963 году, когда академик Андрей Николаевич Тихонов разработал теорию регуляризации. Идея Тихонова гениальна и философски глубока: если точная математическая задача неустойчива, давайте заменим ее близкой задачей, которая устойчива, и пожертвуем математической точностью ради физического смысла.

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

Выбор параметра регуляризации и L-кривая

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

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

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

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