Алгоритм построения приближенных формул для решения прямой кинетической задачи на основе численного решения Текст научной статьи по специальности «Математика»

УДК 541.127 : 519.62 : 517.928

АЛГОРИТМ ПОСТРОЕНИЯ ПРИБЛИЖЕННЫХ ФОРМУЛ ДЛЯ РЕШЕНИЯ ПРЯМОЙ КИНЕТИЧЕСКОЙ ЗАДАЧИ НА ОСНОВЕ ЧИСЛЕННОГО РЕШЕНИЯ

© А. В. Тропин

Башкирский государственный университет Россия, Республика Башкортостан, г. Уфа, 450074, ул. Фрунзе, 32.

Тел./факс: + 7 (34 7) 273 6 7 78.

E-mail: TropinA V@rambler.ru

В статье, первой из трёх по данной теме, предложен алгоритм построения глобальных приближенных формул для решения сложной прямой кинетической задачи. Алгоритм использует численное решение. Глобальные формулы получаются в результате комбинирования локальных. Каждая локальная формула является точным решением задачи, упрощённой по определённому правилу на некотором интервале времени. Относительная погрешность глобальных формул слабо зависит от значений кинетических параметров.

Ключевые слова: прямая кинетическая задача, жесткость, численное решение, точное решение, приближенные формулы, правило упрощения.

Введение

Сложные химические реакции, как правило, протекают с участием реагентов различной активности: от стабильных реагентов до радикалов. Это обусловливает разное влияние каждой из одновременно протекающих элементарных стадий на изменение концентрации того или иного вещества, причём само влияние обычно значительно изменяется в ходе реакции.

С одной стороны, это приводит к жёсткости соответствующей системы обыкновенных дифференциальных уравнений (СОДУ) химической кинетики и трудностям при её численном решении. Жесткость, прежде всего, выражается в неустойчивости алгоритма численного интегрирования по явным схемам, его чувствительности к малым погрешностям каждого шага. Жёсткие системы требуют достаточно малого шага интегрирования, и отсутствует возможность его увеличивать, несмотря на малое изменение концентраций. Увеличение шага сопровождается скачками в численном решении. Основная причина - сильная устойчивость точного решения жёсткой системы, любое отклонение от которого сопровождается резким изменением некоторых производных, стремящихся "вернуть отклонившееся решение в русло точного решения" [1].

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

В работе [3] были попытки использовать асимптотические методы (а именно метод сращивания) для более тонкого исследования СОДУ на различных временных масштабах. Эти методы широко используются в физике при решении систем дифференциальных уравнений с малым параметром [4].

В данной работе предлагается алгоритм получения приближённых формул для решения прямой кинетической задачи, основанный на неравнозначности влияния элементарных стадий на изменение концентраций различных веществ в различное время. Этот алгоритм (и соответствующая программа на Мар1е) появился как результат многочисленных попыток (например, [5]) автоматизировать асимптотические методы применительно к прямой задаче химической кинетики, но без использования малого параметра.

Постановка задачи

Рассматривается абстрактный механизм реак-

ции

k.

n

X а. . ■ A.

j, i i

i = 1

k -

i =1

j, i

■ A j = 1,..., r,

(1)

где Л, (/=1,...,п) - наименование веществ в системе; г - число элементарных стадий; а,-,,-, в,,, - стехиометрические коэффициенты (количества молекул /'-го вещества, участвующего в ,-ой реакции; обычно равны 1, 2, реже 3, иногда допускается 1/2 и другие); / к/ - константы (коэффициенты) скорости/ой стадии в прямом и обратном направлении соответственно.

тг 0 / 0 0 0 \

Пусть задан вектор а = (а ь а 2, ... ,а п) начальных концентраций веществ Л,. Тогда (при некоторых естественных допущениях [6]) в качестве модели химического процесса (1) можно рассматривать задачу Коши (СОДУ с начальными

16

раздел МАТЕМАТИКА

условиями)

^ = г^ • w(a(t)), а(0) = а0, ^

где а(/) = (а^/), а2(/),..., ап(/)) - вектор-функция концентраций веществ Л, в момент времени /; Гт -транспонированная стехиометрическая матрица размера п*г с компонентами у,, = в,/ - а/ w(a(/)) -вектор-столбец скоростей элементарных стадий, компоненты которого, согласно закону действующих масс [6], вычисляются по формуле

w (a(t)) = k + ■ П a (t) j 1 - k . ■ П a (t)

b

j, 1

j

j

i = 1

j = 1,..., r.

(3)

i =1

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

Алгоритм и обсуждение его реализации в системе Maple

1. По заданной схеме реакций и начальным концентрациям строится задача Коши. Эта задача предварительно численно интегрируется каким-либо методом для решения жёстких систем. Шаг интегрирования желательно переменный, обратно пропорциональный наибольшей из скоростей изменения концентраций. Наиболее удобно брать равномерный шаг по логарифмическому времени. Встроенные в Maple процедуры численного решения жёстких систем оказались непригодными, поэтому использовался обычный явный метод Рунге-Кутта с переменным шагом, но с неявной схемой Эйлера для “быстрых” переменных. Опыт автора показывает, что такое сочетание наиболее эффективно для СОДУ химической кинетики. Точность счёта проверяется различным дроблением шага интегрирования. Этот шаг алгоритма полностью автоматический, поэтому соответствующая часть программы не требует вмешательства человека.

2. Задаётся порог малости, например е=0.01. На основе численного решения выясняется, на каких интервалах времени и какие члены в СОДУ могут быть отброшены так, чтобы при этом возмущение решения, вносимое этим упрощением системы, было порядка s. В работе [7] выяснено, как именно следует упрощать систему. Правило упрощения формулируется следующим образом.

Если для некоторого i, при te [a,b) выполняется неравенство

X g i • w. (a(t))

j = 1 ^j j

< e ■ max

j

g . ■ w . (a(t))

i, j j

то ,-е дифференциальное уравнение на [а,Ь) можно (но не обязательно) заменить на алгебраическое г

0 = Е Г- , ■ (а({)).

/ = 1 Ь/ /

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

g . ■ w . (a(t)) i, j j

g . ■ w . (a(t)) i, j j

<£-тах;

/

то слагаемое Уу^,(а(/)) в этом алгебраическом уравнении можно удалить.

Если для некоторых , и / при / из некоторого интервала времени выполняются неравенства

X g ■■w ,(a(t)) i I, j j

j = 1

> e ■ max

j

g

w j (a(t))

g . ■w ,(a(t)) < e^max<

\1’ j j I i

X g f- w (a(t))/max(a. (t)n- max(a. (t)), f = 1 .J J t ‘ I t ‘

то на этом интервале времени в i-м дифференциальном уравнении слагаемое yij-wj(a(t)) можно (но не обязательно) удалить. В теории дифференциальных уравнений с малым параметром это соответствует регулярному возмущению уравнения. Данный шаг алгоритма также полностью автоматический.

3. Исходя из данных, полученных в предыдущем пункте, интервал времени [0,да) разбивается на цепочку подынтервалов, на каждом из которых исходная система упрощается до подсистемы дифференциально-алгебраических уравнений, удовлетворяющей следующим двум условиям:

a) возмущение решения, вносимое упрощением системы, мало;

b) упрощённая подсистема допускает точное интегрирование.

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

4. Подсистемы цепочки интегрируются и корректируются. К сожалению, система Maple практически всегда представляет результаты символьных вычислений в неудачном виде, а иногда отказывается выполнять элементарные упрощения, поэтому результат приходится доводить вручную.

а

r

r

r

Maple очень плохо решает системы дифференциально-алгебраических уравнений, так как алгебраические подсистемы имеют множество решений, поэтому в программе алгебраические подсистемы решаются отдельной процедурой, учитывающей, что решения должны быть положительными.

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

6. Выясняется область значений коэффициен-

тов скоростей kj+, kj и других параметров, при которых полученные приближённые аналитические формулы хорошо аппроксимируют точное решение исходной задачи. Если при каком-то варьировании параметров резко возрастает относительная погрешность приближенных формул, то иногда их удаётся подправить, повторив пункты 3 - 5. Для этого следует выяснить: на каком интервале

времени ошибка велика и для какой неизвестной? Далее, на соответствующем интервале, уточнить упрощённую подсистему (например, оставить член в СОДУ, который по правилу упрощения можно удалить). Затем, откорректировав комбинированные приближённые формулы, опять проверить их зависимость от параметров.

Выводы

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

диапазонах, несмотря на то, что само решение при этом существенно изменялось. Если при каком-то наборе значений параметров формулы для решения оказывались непригодными, то нетрудно было подобрать другие или исправить старые. Для этого достаточно снова применить алгоритм с новым набором значений параметров.

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

Разумно предположить, что данная методика применима не только в химической кинетике, но и в других областях, где моделями являются сложные системы нелинейных обыкновенных дифференциальных уравнений, и, как показывает практика, предпочтительно жёсткие.

Методика позволяет вести качественный анализ процессов, потому что при упрощении теряется лишь часть информации. В отличие от численного решения приближённые формулы не являются просто количественной аппроксимацией информации о процессе. Например, приближённые формулы удобно использовать при решении обратной кинетической задачи. Совершенно очевидно, что при этом резко сокращается объем оптимизационных вычислений: вместо многократного численного решения прямой кинетической задачи -решается одна задача нелинейной регрессии.

ЛИТЕРАТУРА

1. Деккер К., Вервер Я. Устойчивость методов Рунге-Кутты для жёстких нелинейных дифференциальных уравнений. М.: Мир, 1988.- ЗЗ6 с.

2. Семёнов Н. Н. Цепные реакции. Л.: Госхимтехиздат, 19З4.

3. Калякин Л. А., Масленников С. И., Комиссаров В. Д. Анализ механизма жидкофазного ингибированного окисления методом асимптотических приближений. Уфа, 1987.- 41 с. Деп. в ВИНИТИ 1З.07.1987. №4998-в 87.

4. Найфэ А. X. Методы возмущений. М.: Мир, 1976.- 456 с.

5. Тропин А. В., Спивак С. И. // Кинетика и катализ, 1995. Т.З6. №5. С.658-664.

6. Вант Гофф Я. Г. Очерки по химической динамике. Л.: ОНТИ ХИМТЕОРЕТ, 19З6.- 178 с.

7. Тропин А. В. Редукция математических моделей механизмов цепных реакций: дис. ... канд. физ.-мат. наук. Уфа. 1998.- 104 с.

Поступила в редакцию 23.03.2007 г.