Решение жестких систем обыкновенных дифференциальных уравнений (ОДУ)
Феномен жесткости: когда классические методы дают сбой
При численном моделировании многих химических, биологических и электронных систем исследователи сталкиваются с загадочным и крайне неприятным вычислительным явлением, получившим название «жесткость» (stiffness). Представьте себе химический реактор, в котором протекают две параллельные реакции: одна завершается за миллисекунды (например, горение или взрыв), а вторая длится часами (медленное окисление или диффузия). Эти процессы описываются системой обыкновенных дифференциальных уравнений (ОДУ).
Математически система называется жесткой, если ее решения содержат компоненты с кардинально различающимися скоростями затухания. Если мы попытаемся решить такую систему классическим явным методом (например, популярным методом Эйлера или методом Рунге-Кутты 4-го порядка), мы столкнемся с катастрофой. Для того чтобы обеспечить вычислительную устойчивость и не дать быстрым компонентам «взорвать» решение, шаг интегрирования должен быть меньше характерного времени самого быстрого процесса. В результате, чтобы промоделировать медленный процесс длительностью в часы, компьютеру придется сделать триллионы крошечных миллисекундных шагов, что займет месяцы машинного времени и приведет к огромному накоплению ошибок округления.
Неявные методы и А-устойчивость
Для преодоления проблемы жесткости классические явные методы категорически не подходят. Спасением становятся неявные (имплицитные) методы интегрирования. В явных методах значение функции на следующем временном шаге вычисляется только на основе уже известных данных с предыдущих шагов. В неявных же методах неизвестное значение входит не только в левую, но и в правую часть уравнения. Это превращает каждый шаг интегрирования в задачу поиска корня нелинейного алгебраического уравнения.
На первый взгляд это кажется усложнением: теперь на каждом шаге нужно запускать метод Ньютона и решать системы линейных уравнений. Однако неявные методы обладают потрясающим свойством — абсолютной устойчивостью (или А-устойчивостью). Для таких методов размер шага интегрирования ограничивается только требованиями к точности описания физического процесса, а не условиями математической устойчивости. Медленный процесс можно интегрировать огромными шагами, полностью игнорируя сверхбыстрые затухающие колебания, не боясь при этом «взрыва» алгоритма. Неявный метод Эйлера является простейшим представителем этого класса.
Формулы дифференцирования назад (BDF) и методы Гира
Простейшие неявные методы обладают низким порядком точности (первым или вторым). Для точного моделирования жестких систем требуются методы более высоких порядков. Американский математик Чарльз Гир в начале 1970-х годов разработал семейство жестко-устойчивых алгоритмов, известных как формулы дифференцирования назад (Backward Differentiation Formulas, BDF).
Вместо того чтобы аппроксимировать подынтегральную функцию правой части, методы BDF строят интерполяционный полином по нескольким последним вычисленным точкам самого решения, а затем дифференцируют этот полином и приравнивают его к значению правой части в новой, еще не вычисленной точке. Это позволяет создавать высокоточные многошаговые неявные схемы (вплоть до 6-го порядка). Современные программные библиотеки (например, LSODE или решатели ode15s в MATLAB) используют адаптивные алгоритмы на основе методов BDF. Они автоматически распознают жесткость системы, на лету меняя как размер шага, так и порядок метода, обеспечивая робастное и эффективное решение сложнейших многомасштабных задач химической кинетики и схемотехники.