Когда в уравнения вмешивается случайность

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

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

Метод Эйлера-Маруямы: база стохастического моделирования

Самым простым и исторически первым численным методом для СДУ является схема Эйлера-Маруямы (названная в честь Леонарда Эйлера и японского математика Гисиро Маруямы). Этот метод является естественным стохастическим обобщением классического метода Эйлера для обычных ОДУ.

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

Метод Мильштейна и сильная vs слабая сходимость

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

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

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

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

Допустим, у нас есть простой алгоритм для вычисления какой-либо величины (производной, определенного интеграла или решения ОДУ), который имеет низкий порядок точности. Пусть его ошибка пропорциональна квадрату шага сетки O(h^2). Если мы посчитаем значение с шагом h, мы получим грубый ответ. Если мы посчитаем значение с вдвое меньшим шагом (h/2), мы получим более точный ответ (ошибка уменьшится в 4 раза). Идея Ричардсона заключается в том, чтобы взять эти два приближенных результата и алгебраически скомбинировать их так, чтобы главный член ошибки O(h^2) взаимно уничтожился (вычелся). В результате мы получаем новую оценку, погрешность которой пропорциональна уже O(h^4) — то есть алгоритм второго порядка мгновенно превращается в алгоритм четвертого порядка практически бесплатно!

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

Наиболее известное и триумфальное применение экстраполяции Ричардсона было предложено Вернером Ромбергом в 1955 году для задачи численного интегрирования. Алгоритм получил название метод (или процесс) Ромберга. Базой для этого метода служит самая простая и грубая квадратурная формула — составной метод трапеций.

Формула Эйлера-Маклорена доказывает, что ошибка метода трапеций раскладывается в бесконечный ряд только по четным степеням шага h (h^2, h^4, h^6 и т.д.). Алгоритм Ромберга использует этот факт на полную мощность. Вычисления организуются в виде треугольной таблицы. В первый столбец записываются результаты вычисления интеграла методом трапеций с последовательным делением шага пополам (1 разбиение, 2, 4, 8, 16...). Это самая ресурсоемкая часть алгоритма. Затем начинается магия экстраполяции.

Построение треугольной таблицы и феноменальная сходимость

Второй столбец таблицы формируется путем применения экстраполяции Ричардсона к соседним элементам первого столбца. Удивительно, но элементы второго столбца математически строго совпадают с результатами, которые дал бы метод Симпсона (точность O(h^4)). Но мы на этом не останавливаемся! Мы применяем ту же формулу экстраполяции (с измененными коэффициентами) к элементам второго столбца, формируя третий столбец, ошибка которого имеет порядок O(h^6) (метод Буля). Процесс продолжается до тех пор, пока мы не дойдем до вершины треугольника.

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

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

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

Для формального анализа устойчивости в вычислительной математике используется простое линейное тестовое уравнение Далквиста: 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-устойчивые алгоритмы (например, неявный метод Эйлера) мгновенно обнуляют высокочастотные компоненты, что делает их идеальными для моделирования ударных волн и резких переключателей в транзисторах.