Solve_ODE (972307), страница 3
Текст из файла (страница 3)
Для дифференциальных уравнений химической кинетики развивались специфические численные методы, использовавшие особенности формы этих уравнений. Например, учитывался тот факт, что форма решений хорошо передается линейной комбинацией экспоненциальных функций (см. выше пример системы ОДУ с постоянными коэффициентами (11)). Предлагались также подходы, основанные на временной замене дифференциальных уравнений с большими алгебраическими уравнениями подобно тому, как это делается в приближении квазистационарных концентраций. Некоторые из развитых подходов оказались достаточно успешными, хотя область их применимости, как правило, ограничена частными задачами.
Только на рубеже 70‑х годов появился действительно эффективный универсальный метод интегрирования жестких систем ОДУ произвольной структуры, предложенный Гиром [5-7]. Он стал, по существу, единственным методом, применяемым в современных программах для математического моделирования кинетики химических процессов [3, 4, 8].
Метод Гира разработан на основе теоретических исследований проблемы устойчивости численных решений дифференциальных уравнений, проводившихся в 60‑х годах. Было сформулировано понятие жесткой устойчивости, которая является более слабой по сравнению с A‑устойчивостью, но тем не менее позволяет интегрировать жесткие системы ОДУ, ограничивая шаг только условием приемлемости локальной ошибки. Оказалось, что для многошаговых схем, обладающих жесткой устойчивостью, нет ограничений на порядок метода (в отличие от A‑устойчивых методов), причем действительно удалось построить такие схемы 3‑го, 4‑го и более высоких порядков.
Метод Гира представляет собой вариант метода прогноза и коррекции переменного порядка (от 1 до 6) с автоматическим выбором шага. Переменный порядок фактически означает, что имеется не одна, а целый набор формул интегрирования (по паре для каждого порядка). Прогноз проводится по явным формулам, которые получаются усечением ряда Тейлора (2) до членов соответствующего порядка. Необходимые значения высших производных последовательно оцениваются через разделенные разности первой, второй и т. д. производных в предыдущих точках (т. наз. схема Нордсика). Для коррекции применяются неявные формулы до 6‑го порядка включительно, обладающие свойством жесткой устойчивости (по мере роста порядка область устойчивости несколько сужается). Коррекция выполняется с помощью нескольких ньютоновских итераций (в отличие от однократного применения простой итерации в других методах прогноза и коррекции, таких, как методы Милна или Адамса). Это снижает требования к точности прогноза и обеспечивает дополнительную устойчивость процесса интегрирования.
Эффективное интегрирование жестких уравнений обеспечивается специальной стратегией одновременного управления величиной шага и порядком метода: в каждой точке выбирается оптимальный порядок, который позволяет использовать максимальный шаг при условиях сохранения устойчивости и обеспечения заданного уровня локальной погрешности решения. Переменный порядок делает метод Гира самостартующим, поскольку процесс интегрирования всегда начинается с одноточечной схемы первого порядка.
Высокую эффективность метода Гира в задачах химической кинетики можно проиллюстрировать на примере модели радиолиза водного раствора метанола [4], упоминавшейся в предыдущем разделе. Для интегрирования до t=4·105 с (предельно достижимый интервал в случае использования обычного метода Рунге-Кутта) потребовалось сделать 166 шагов, а для достижения точки t=100 с — 806 шагов; при этом время счета на ЭВМ БЭСМ‑6 (1 Мфлоп/с) составило около 15 с. Заметим, что хотя метод Рунге-Кутта сохранял устойчивость при шаге до 109 с, полное совпадение решений с методом Гира было получено только при шаге 1011 с.
Литература
1. G. Dahlquist, BIT, 1963, v. 3, p. 27.
2. Э. Камке, «Справочник по обыкновенным дифференциальным уравнениям». М.: 1971.
3. R. G. Gelinas, J. Comput. Phys., 1973, v. 11, p. 455.
4. А. В. Абраменков, И. А. Абраменкова, Химия высоких энергий, 1979, т. 13, с. 557 (полный текст статьи депонирован в ВИНИТИ за № 1943-79, 76 с.)
5. C. W. Gear, Information Processing 68. Amsterdam, 1969, p. 187.
6. C. W. Gear, IEEE Trans., 1971, CT 18, p. 89.
7. C. W. Gear, Comm. ACM, 1971, v. 14, p. 176.
8. Л. С. Полак, М. Я. Гольденберг, А. А. Левицкий. «Вычислительные методы в химической кинетике». М.:, 1984.