Б методы решения жестких систем обыкновенных дифференциальных уравнений

Явление жесткости. Предварительные сведения

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

Рассмотрим две задачи Коши для систем ОДУ^[1]:

с начальными данными и (0) = м₍₎, ^(0) = v₀ (здесь константы я ~ 0(1),? «; 1) и.

Решением первой задачи Коши являются функции а второй —.

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

(u — вектор-столбец, A — матрица с постоянными коэффициентами) существенно различаются. Так, в первом случае Х_{" -е^-1, Л₂~ 0(1); во втором — ~ -1000, Х₂ = -1. В обоих случаях имеем.

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

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

где u_k — концентрации веществ, участвующих в химических реакциях, скорости протекания которых характеризуются коэффициентами cf и ajj. В качестве примера приведем одну из систем химической кинетики, описывающую изменение концентрации трех веществ, участвующих в реакции для случая полного перемешивания^[2].

Пример 6.1.

Обозначим концентрации трех веществ, участвующих в реакциях, через u_t, и₂ и и₃, тогда имеем

Участки решения, характеризующиеся быстрым и медленным его изменением, называются пограничным слоем и квазистационарным режимом соответственно.

Трудности численного решения подобных систем ОДУ, получивших название жестких (определение жесткой системы приведено ниже), связаны с выбором шага интегрирования. Дело в том, что характерные времена исследуемых процессов могут различаться более чем в 10¹² раз. Следовательно, если при численном решении системы.

выбирать шаг из условия, диктуемого теоремами об устойчивости численных методов (для методов типа явных методов Рунге — Кутты соответствующие теоремы приведены в гл. 5),.

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

По этой причине есть следующие альтернативы в выборе подхода к численному решению рассматриваемых задач.

1. Численно решать систему ОДУ с шагом.

т.е. с учетом характерных времен всех процессов, описываемых данной системой.

• 2. Решать систему ОДУ с разными шагами, соответствующими физическим процессам с существенно различными характерными временами. В этом случае необходимо задавать условия перехода к другому шагу интегрирования.
• 3. Пренебречь быстропротекающими процессами и численно рассматривать лишь медленные, проводя интегрирование с шагом, превышающим характерные времена быстрых процессов. В этом случае нам придется конструировать численные методы, позволяющие проводить расчеты с шагом

Дадим теперь формальное определение жесткой системы ОДУ^[3]. Рассмотрим вначале случай автономной системы.

Определение 6.1. Система ОДУ для задачи Коши.

называется жесткой, если спектр матрицы Якоби J = (f'(и)} разделяется на две части:

• 1) жесткий спектр: ReA^(w) < -Л₀, 11шЛ^| < | ReAj, k = 1,…, N_{ — собственные значения матрицы Якоби), А₀ > 0; 7^А₀ «1;
• 2) мягкий спектр: |Х_;-| < X₀, j = 1,…, N₂, Х₀ > 0, при этом Х₀ «с А₀. T_kX₀ — 1.

Отношение А₀/Х₀ называется показателем жесткости системы. В дальнейшем будем полагать Х₀ ~ 0(1).

Определение жесткой системы (ЖС) легко обобщается и на случай неавтономной системы. Естественно считать жесткой такую неавтономную систему, у которой спектр ее матрицы Якоби делится на жесткую и мягкую части при каком-нибудь значении времени t, 0 < t < Т_к. На практике встречаются задачи, являющиеся жесткими в какой-либо части фазового пространства и не являющиеся жесткими в остальном фазовом пространстве. Неавтономная система может менять свою жесткость в зависимости от времени.

Кроме того, как отмечалось выше, жесткость системы зависит и от длины отрезка интегрирования (интервала времени, на котором решается система). Если, как было отмечено в определении, спектр задачи расщепляется, но при этом T_kA_{)< 1, то все неприятности, проявляющиеся при численном решении ЖС ОДУ, несущественны. Систему можно интегрировать любым численным методом с достаточно малым шагом.

Проблему численного решения жестких систем ОДУ рассмотрим на примере модельной линейной системы вида.

Ее точное решение задается формулой.

где константы интегрирования Ср C_k соответствуют жесткой и мягкой частям спектра; Sip со^ — собственные векторы матрицы Якоби, соответствующие собственным значениям Ар X_fr

В этом решении видны две части: первая (жесткая) убывает как на временном интервале [^^(Aq¹)] (пограничный слой), вторая заметно изменяется на интервале [^^(А^¹)] (квазистационарный режим).

Если провести аппроксимацию линейной системы ОДУ с помощью явного метода Эйлера:

то общее решение такой системы разностных уравнений будет.

Второе слагаемое в этом решении аппроксимирует второе слагаемое в точном решении (6.2), а первое быстро растет. Этот нефизический рост приводит к абсурдному результату.

Теперь проведем аппроксимацию линейной системы ОДУ (6.1) с помощью неявного метода Эйлера:

Общее решение такого разностного уравнения.

В этом случае второе слагаемое ведет себя так же, как и точное решение, а первое стремится к нулю как (тЛ₀)^-'⁷, т. е. его поведение качественно совпадает с точным в области пограничного слоя.

В практике численных исследований жестких задач часто оказывается ненужным изучение поведения решения в пограничном слое, и мы можем воспользоваться неявными методами. Но в случае необходимости исследовать этот слой можно с шагом т «с Aq^[4].

Устойчивость методов численного интегрирования жестких систем ОДУ обычно исследуется на примере скалярного уравнения.

Положим, что численный метод, применяемый к решению этого уравнения, может быть записан в виде.

где R (z) называется функцией устойчивости^[4]. О построении функции устойчивости для методов Рунге — Кутты речь пойдет ниже.

Определение 6.2. Численный метод для решения уравнения (6.4) является абсолютно устойчивым, если выполнено условие |R (z)| < 1.

Из определения следует, что |м_я+11 < и_п. Это требование является естественным при Rez < 0, поскольку в таком случае модуль точного решения есть невозрастающая функция.

Множество всех точек z, для которых |i?(z)| < 1, называется областью абсолютной устойчивости.

Определение 6.3. Если область абсолютной устойчивости (|R (z)| < 1) занимает левую полуплоскость комплексной плоскости (Rez < 0, заштрихованная область на рис. 6.1), то метод является А-устойчивым.

Рис. 6.1. Область Л-устойчивости.

Если область абсолютной устойчивости включает в себя угол (в левой полуплоскости комплексной плоскости) с вершиной в нуле и углом полураствора а, то метод называется А (а)-устойчивым. Если метод устойчив на левой полуоси действительной оси, то метод А (0)-устойчив.

Области Л (а) — и Л (0)-устойчивости изображены заштрихованными на рис. 6.2, а и б соответственно.

Рис. 6.2. Области Л (а)-устойчивости и Л (0)-устойчивости.

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

Рис. 6.3. Пример области жесткой устойчивости.

Определение 6.4. Численный метод называется L-устойчивым, если он Л-устойчив и выполняется условие |/?(z)| —" 0 при z —> .

В частности, рассмотренный выше неявный метод Эйлера является L-устойчивым. Решения, полученные L-устойчивыми методами, будут затухающими.

Иногда L-устойчивым методом называют жестко устойчивый (но не Л-устойчивый) метод, для которого также выполняются приведенные выше условия.

• [1] Федоренко Р. П. Указ, соч.; Хайрер Э., Баннер Г. Указ, соч.; Современные численные методы решения обыкновенных дифференциальных уравнений / под ред. Дж. Холла, Дж. Уатта. М.: Мир, 1979; Каханер Д., Моулер К., Нэш С. Численные методы и программноеобеспечение. М.: Мир, 1998.
• [2] Кахсшер Д.} Моулер К., Нэш С. Указ. соч.
• [3] Ризничеико Г. Ю. Указ. соч.
• [4] Самарский А. А., Михайлов А. П. Указ, соч.; Базыкин А. Д. Указ. соч.
• [5] Самарский А. А., Михайлов А. П. Указ, соч.; Базыкин А. Д. Указ. соч.