Глава 9.4.3. Пример: химическая кинетика

Рассмотрим классическую модель химической кинетики (Робертсон, 1966), которая как нельзя лучше передает смысл понятия жесткости ОДУ.

Попытка решения стандартными методами

Рассмотрим составную схему химического взаимодействия трех веществ. Пусть вещество "О" медленно превращается в "1": "0"->"1" (со скоростью 0.1), вещество "1" при каталитическом воздействии самого себя превращается очень быстро в вещество "2": "1"+"1"->"2"+"1" (10³). И, наконец, подобным образом (но со средней скоростью) реагируют вещества "2" и "1": "1"+"2"->"0" + "2" (10²). Система ОДУ, описывающая динамику концентрации реагентов, с попыткой решения методом Рунге—Кутты, приведена, в символической форме в первой строке листинга 9.10.

Листинг 9.10. Жесткая система ОДУ химической кинетики

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

Листинг 9.11. Якобиан рассматриваемой системы ОДУ химической кинетики

Для примера, приведенного в листинге 9.10, стандартным методом Рунге— Кутты все-таки удается найти решение (оно показано на рис. 9.15). Однако для этого требуется очень большое число шагов, M=20000, что делает расчеты очень медленными. При меньшем числе шагов численному алгоритму не удается найти решение. В процессе работы алгоритма оно расходится, и Mathcad вместо результата выдает ошибку о превышении предельно большого числа.

Еще один факт, на который стоит обратить внимание, — это различие в порядке величины получающегося решения. Как видно из рис. 9.15, концентрация первого реагента уг существенно (в тысячи раз) превышает концентрацию остальных. Это свойство также очень характерно для жестких систем.

ПРИМЕЧАНИЕ

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

Рис. 9.15. Решение жесткой системы ОДУ химической кинетики методом Рунге—Кутты (продолжение листинга 9.10)

Применение специфических алгоритмов

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

Листинг 9.12. Решение жесткой системы ОДУ химической кинетики

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

Важно заметить, что до сих пор мы имели дело с примером не очень жесткой системы. Попробуйте вместо скоростей упомянутых химических реакций 0.1, 10³ и 10² взять другие числа, например, 0.05, 10⁴ и 10⁷ соответственно. Заметим, что такое соотношение скоростей часто встречается в прикладных задачах химической кинетики и определяет куда более жесткую систему ОДУ. Ее уже никак не удается решить стандартными методами, поскольку число шагов численного метода должно быть просто гигантским. А между тем алгоритмы для жестких ОДУ справляются с этой задачей с легкостью (рис. 9.16), причем практически при тех же значениях шага, что были взяты в листинге 9.12. Обратите внимание, что порядки величины решений для концентраций различных веществ на рис. 9.16 различаются еще сильнее, чем в предыдущем (менее жестком) примере.

Рис. 9.16. Решение более жесткой системы ОДУ химической кинетики методом Розенброка

Для решения очень жестких систем особенно подходит функция Radau, которая соблазнительна еще и тем, что избавляет от необходимости предварительного расчета якобиана. В случае очень жесткой задачи химической кинетики результат, выдаваемый функцией Radau, практически тот же, что и в случае использования алгоритма Розенброка (рис. 9.16). Любопытно, что решение менее жесткой системы (из предыдущего листинга) при помощи функции Radau удается получить только для первой половины интервала, примерно до t1=20. Для больших интервалов она дает сбой, выводя вместо решения сообщение об ошибке.

Содержание раздела