Что такое жесткие системы уравнений

Жесткие системы ОДУ

До сих пор мы имели дело с «хорошими» уравнениями, которые надежно решались численными методами (типа Рунге-Кутты). Однако имеется класс так называемых жестких (stiff) систем ОДУ, для которых стандартные методы практически неприменимы, поскольку их решение требует исключительно малого значения шага численного метода. Для этих систем разработаны специальные алгоритмы

Что такое жесткие ОДУ?
Строгого общепринятого математического определения жестких ОДУ нет. В рамках этой книги будем считать, что жесткие системы — это те уравнения, решение которых получить намного проще с помощью определенных неявных методов, чем с помощью явных методов (типа тех, что были рассмотрены в предыдущих разделах). Примерно такое определение было предложено в 1950-х годах классиками в этой области Кертиссом и Хиршфельдером. Начнем разговор о жестких ОДУ с примера нежесткого уравнения, решение которого показано на рис.1. На том же графике показано решение подобного ОДУ с другим коэффициентом при правой части, равным не -10, а -30. Решение обоих уравнений не составило труда, и численный метод Рунге-Кутты дал правильный результат.

Рис. 1. Решение нежестких ОДУ методом Рунге-Кутты

На рис. 2 показано решение того же ОДУ с коэффициентом -50. Вас, несомненно, должен насторожить результат, выданный методом Рунге-Кутты. Характерная «разболтка» решения говорит о неустойчивости алгоритма. Первое, что можно сделать, — увеличить количество шагов в методе Рунге-Кутты. Можно убедиться, что при step>20 разболтка пропадает, и решение становится похожим на графики, показанные на рис. 1.

Рис. 2. Неверное решение более жесткого ОДУ методом Рунге-Кутты

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

Замечание Некоторые ученые замечают, что в последние годы методы Рунге-Кутты стали уступать свое главенствующее положение среди алгоритмов решения ОДУ методам, способным решать жесткие задачи.

Исторически, интерес к жестким системам возник в середине XX века при изучении уравнений химической кинетики с одновременным присутствием очень медленно и очень быстро протекающих химических реакций. Тогда неожиданно оказалось, что, как считалось, исключительно надежные методы Рунге-Кутты стали давать сбой при расчете этих задач. Рассмотрим классическую модель взаимодействия трех веществ (Робертсон, 1966), которая как нельзя лучше передает смысл понятия жесткости ОДУ.
Пусть вещество «0» медленно превращается в «1»: «0»—>»1″ (со скоростью 0.1), вещество «1» при каталитическом воздействии самого себя превращается очень быстро в вещество «2»: «1»+»1″—>»2″+»1″ (10 3 ). И, наконец подобным образом (но со средней скоростью) реагируют вещества «2» и «1»: «1»+»2″—>»0″+»2″ (10 2 ). Система ОДУ, описывающая динамику концентрации реагентов, приведена ниже.

Бросается в глаза сильно различающийся порядок коэффициентов при разных слагаемых. Именно степень этого различия чаще всего и определяет жесткость системы ОДУ. В качестве соответствующей характеристики выбирают матрицу Якоби (якобиан) векторной функции F(t,y), т.е. функциональную матрицу, составленную из производных F(t,y). Чем вырожденнее матрица Якоби, тем жестче система уравнений. В приведенном примере определитель якобиана и вовсе равен нулю при любых значениях y₀, y₁ и y₂:

Для приведенного примера стандартным методом Рунге-Кутты все-таки удается найти решение. Однако, для этого требуется очень большое число шагов, M=20000, что делает расчеты очень медленными. При меньшем числе шагов численному алгоритму не удается найти решение. Решение данной системы ОДУ приведено в разделе Модель химической кинетики главы Динамические системы.

В принципе, можно было бы снизить жесткость системы «вручную», применяя масштабирование. Для этого нужно искусственно уменьшить искомую функцию y₁, к примеру, в тысячу раз, разделив все слагаемые в системе ОДУ, содержащие y_1, на 1000. После масштабирования для решения полученной системы методом Рунге-Кутты будет достаточно взять всего M=20 шагов.

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

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

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

Рассмотрим в качестве примера две задачи Коши для систем обыкновенных дифференциальных уравнений (ОДУ) [9.1], [9.2]:

с начальными данными u(0) = u₀, v(0) = v₀ ; здесь ; и линейную систему с постоянными коэффициентами

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

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

( u — вектор — столбец, A — матрица с постоянными коэффициентами) существенно различаются. Так, в первом случае ; во втором: В обоих случаях имеем:

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

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

Пример 1. Обозначим концентрации трех веществ, участвующих в реакции, через u₁ , u₂ и u₃ , тогда

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

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

выбирать шаг из условия

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

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

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

Что такое жесткие системы уравнений

Глава 5. Решение дифференциальных уравнений

5.6 Жесткие системы дифференциальных уравнений

Система дифференциальных уравнений, записанная в матричном виде, Y = A X , считается жесткой, если матрица коэффициентов почти вырожденная. В этом случае решение, возвращаемое функцией rkfixed , может быть неустойчивым. При решении жесткой системы необходимо использовать одну из трех функций, специально разработанных для решения жестких систем дифференциальных уравнений. К ранее существовавшим функциям Stiffb и Stiffr , реализующим методы Булирша – Штера или Розенброка соответственно, в MathCAD 2001 i добавилась функция Radau , реализующая метод RADAU 5. рассмотрим вначале функции Stiffb и Stiffr .

Вид матрицы, возвращаемой этими функциями, совпадает с возвращаемым функцией rkfixed . Однако функции Stiffb и Stiffr требуют дополнительно задания якобиана системы уравнений.

Обращение к этим функциям:

Здесь у – вектор начальных условий размерности m , где m – порядок ОДУ или число уравнений в системе ОДУ; х₁ и х₂ – начало и конец интервала интегрирования, на котором ищется решение системы ОДУ; начальные условия, заданные вектором у, – это значение решения системы в точке х₁; n – число точек (не считая начальной), в которой ищется решение; D ( x , y ) – m – мерный вектор, который содержит первые производные неизвестных функций; J ( x , y ) – функция, которая возвращает матрицу размером . Первый столбец ее содержит производные F ( t , y ) по у. Остальные столбцы представляют собой матрицу Якоби системы ОДУ.

Пример решения жесткой системы ОДУ приведен на рис. 5.14.

Рис. 5. 14 Пример решения жесткой системы ОДУ

Следует заметить, что увеличение числа шагов в 1000 раз позволяет решить приведенную систему с помощью функции rkfixed ( рис. 5.15), но это оттого, что заданная система не слишком жесткая. Если увеличить разницу в коэффициентах уравнений, то функция rkfixed потребует гигантского числа шагов и в принципе не сможет решить систему.

Рис. 5. 15 Решение жесткой системы ОДУ функцией rkfixed

Новая функция Radau , введенная в MathCAD 2001 i , имеет такой же список аргументов, как и функция Rkadapt и Bulstoer , а именно: Radau (у, x ₁ , x ₂ , D ).

Функция Radau предназначена для решения систем жестких ОДУ, как и функции Stiffr и Stiffb . Преимуществом функции Radau перед Stiffr и Stiffb является то, что она требует указания якобиана в качестве параметра функции.

Сравнение результатов расчета жесткой системы ОДУ с помощью функций Radau и Stiffr приведено на рис. 5 . 16 . Численные результаты, естественно, совпадают, но использование функции Radau проще.