Решение жесткой системы ОДУ химической кинетики
Листинг 11.14. Решение жесткой системы ОДУ химической кинетики
Расчеты показывают, что для получения того же результата (см. рис. 11.18) оказалось достаточно в тысячу раз меньшего количества шагов численного алгоритма, чем для стандартного метода Рунге-Кутты! Примерно во столько же раз требуется меньше компьютерного времени на проведение расчетов. Стоит ли говорить, что, если Вы имеете дело с жесткими (в той или иной степени) системами, применение описанных специальных алгоритмов просто необходимо.
Важно заметить, что до сих пор мы имели дело с примером не очень жесткой системы. Попробуйте вместо скоростей упомянутых химических реакций (см. разд. 11.5.1), 0.1, 103 и 102 взять другие числа, например 0.05, 10 и 107, соответственно. Заметим, что такое соотношение скоростей часто встречается в прикладных задачах химической кинетики и определяет куда более жесткую систему ОДУ. Ее уже никак не удается решить стандартными методами, поскольку число шагов численного метода должно быть просто гигантским. А между тем, алгоритмы для жестких ОДУ справляются с этой задачей с легкостью (рис. 11.19), причем практически при тех же значениях шага, что были взяты в листинге 11.14. Обратите внимание, что порядки величины решений для концентраций различных веществ на рис. 11.19 различаются еще сильнее, чем в предыдущем (менее жестком) примере.
Это еще раз доказывает, что одна и та же система ОДУ с различными коэффициентами может быть жесткой в разной степени В частности, приведенный выше пример генератора Ван дер Поля с параметром ц=5000 — это уже пример жесткой задачи.
В заключение приведем соответствующие встроенные функции, которые применяются для решения жестких систем ОДУ не на всем интервале, а только в одной заданной точке t1.
- radau(y0, t0, t1,acc,F,k,s) —метод RADAUS
- stiffb(y0,t0,t1,acc,F,j,k,s) — метод Булирша-Штера
- stiffr(y0,t0,t1,acc,F,j,k,s) —метод Розенброка
Имена этих функций пишутся с маленькой буквы, а их действие и набор параметров полностью аналогичны рассмотренным нами ранее для функций, относящихся к решению в заданной точке нежестких систем (см. разд. 11.3.2). Отличие заключается в специфике применяемого алгоритма и необходимости задания матричной функции якобиана j(t,y).
