Метод прямых (Method of Lines) для дифференциальных уравнений в частных производных

Устранение пространственных переменных

При математическом моделировании теплопроводности, химической кинетики в реакторах или процессов диффузии инженеры сталкиваются с эволюционными (нестационарными) уравнениями в частных производных (УЧП) параболического и гиперболического типов. Традиционный подход к их численному решению подразумевает одновременную дискретизацию как по пространственным координатам, так и по времени, что приводит к созданию двумерных разностных схем (например, явной схемы Эйлера или неявной схемы Кранка-Николсон). Однако существует альтернативная, весьма мощная математическая стратегия, известная как метод прямых (Method of Lines, MOL).

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

Дискретизация пространства и появление системы ОДУ

Рассмотрим, например, одномерное уравнение теплопроводности. Мы заменяем вторую пространственную производную в правой части уравнения классическим конечно-разностным приближением (по трем соседним узлам). В результате для каждого внутреннего узла сетки мы получаем отдельное обыкновенное дифференциальное уравнение (ОДУ) первого порядка по времени. Если мы разбили стержень на 1000 узлов, наше одно УЧП чудесным образом превратилось в жестко связанную систему из 1000 ОДУ.

Огромное преимущество такого подхода заключается в том, что теория и программные комплексы для решения систем ОДУ проработаны в вычислительной математике до совершенства. Нам больше не нужно вручную программировать сложные шаги по времени. Мы просто отдаем эту гигантскую систему ОДУ в руки высококлассного, «умного» интегратора (такого как LSODE, VODE или функция ode15s в MATLAB). Этот интегратор возьмет на себя самую сложную работу: он будет автоматически оценивать локальную ошибку, адаптивно менять размер шага по времени (сокращая его при резких скачках температуры и увеличивая в спокойные периоды) и автоматически менять порядок аппроксимации для поддержания идеального баланса между скоростью и точностью.

Проблема жесткости и алгоритмы Гира

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

Если мы попытаемся скормить такую систему обычному явному решателю ОДУ (например, классическому адаптивному методу Рунге-Кутты-Фельберга), он потерпит сокрушительное фиаско. Из-за жесткости алгоритм будет вынужден уменьшить временной шаг до микроскопических долей секунды, чтобы не допустить взрыва вычислительной неустойчивости, и симуляция одной секунды физического процесса займет месяцы машинного времени. Поэтому для успешного применения метода прямых критически важно использовать специализированные неявные решатели ОДУ, основанные на формулах дифференцирования назад (BDF, или методы Гира). Они обладают абсолютной (L-устойчивостью), что позволяет им бесстрашно шагать по оси времени гигантскими адаптивными шагами, игнорируя сверхкороткие численные флуктуации и блестяще решая задачи математической физики.