Численное решение дифференциальных уравнений с запаздывающим аргументом

Когда настоящее зависит от далекого прошлого

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

Такие процессы моделируются дифференциальными уравнениями с запаздывающим аргументом (ДДУЗА, или Delay Differential Equations, DDEs). В простейшем виде уравнение выглядит как y'(t) = f(t, y(t), y(t-tau)), где tau — это время запаздывания. Эта небольшая добавка (t-tau) радикально меняет всю математическую природу задачи, превращая конечномерное пространство состояний ОДУ в бесконечномерное. Теперь для запуска процесса интегрирования нам недостаточно знать начальное положение системы в одной точке (t=0). Нам необходимо задать начальную функцию (предысторию) на всем интервале от -tau до 0!

Распространение разрывов и сглаживание

Первая вычислительная трудность при решении ДДУЗА — это феномен распространения разрывов производных. Если начальная функция на отрезке [-tau, 0] не склеена идеально гладко с самим уравнением в точке t=0 (что бывает почти всегда), то в точке t=0 возникает разрыв первой производной решения. Поскольку значение из точки t=0 попадет в правую часть уравнения через время tau, этот разрыв породит излом уже второй производной в точке t=tau. Затем появится разрыв третьей производной в точке t=2*tau, и так далее.

К счастью, дифференциальный оператор обладает эффектом интегрирования: с каждым шагом на величину запаздывания порядок разрыва повышается, и решение становится все более гладким. Тем не менее, для высокоточных численных методов (например, Рунге-Кутты 4-го порядка) эти начальные изломы являются катастрофой — метод теряет свой порядок точности. Современный решатель ДДУЗА должен уметь алгоритмически выявлять эти точки разрывов (breaking points) и принудительно ставить узлы сетки интегрирования точно в эти точки, чтобы не интегрировать через излом.

Плотный вывод (Dense Output) и интерполяция

Вторая и главная проблема: алгоритму нужно знать значение функции y(t-tau) в прошлом. Если алгоритм интегрирует с переменным адаптивным шагом (что обязательно для эффективных расчетов), точка (t-tau) практически никогда не совпадет с теми узлами сетки, которые алгоритм уже вычислил и сохранил в памяти. Точка запаздывания попадет куда-то между прошлыми узлами.

Чтобы найти значение между узлами, алгоритму необходим механизм непрерывного (плотного) вывода (Dense Output). Классические дискретные методы Рунге-Кутты здесь не подходят. Математикам пришлось разрабатывать непрерывные методы Рунге-Кутты, которые на каждом шаге выдают не только число в конце интервала, но и специальный интерполяционный полином (сплайн Эрмита), который локально описывает функцию внутри шага с сохранением заявленного порядка точности. Обращаясь к истории полиномов, решатель вычисляет точное значение запаздывающего аргумента. Современные программы (например, dde23 в MATLAB) изящно решают системы с постоянным, переменным и даже распределенным (интегральным) запаздыванием, позволяя моделировать хаотические аттракторы Маккея-Гласса в лазерной физике.