Абсолютная устойчивость и барьеры Далквиста в решении дифференциальных уравнений

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

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

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