Автоколебания в системе аскорбиновая кислота – Дегидроаскорбиновая кислота в присутствии оксигенированных комплексов железа(II)

УДК 541.128.7

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

У Г. Магомедбеков

Приведены результаты по исследованию химических колебательных процессов, реализующихся при окислении аскорбиновой кислоты в присутствии оксигенированных комплексов железа(П) в гомогенной среде. Определены пределы концентраций реагента и катализатора, рН и температуры, при которых реализуются концентрационные колебания. На основе Фурье-преобразования и реконструкции динамики по временной последовательности данных определены частоты колебаний, а также размерности фазового пространства и аттрактора. Предложен возможный механизм протекания процессов и составлена математическая модель, описывающая кинетику. Выполнен термодинамический анализ процессов, на основе чего выявлены причины проявления химических нестабильностей. Проведен качественный анализ системы дифференциальных уравнений, описывающей кинетику реакций, и получено их численное решение. Определены характер стационарного состояния и возможность бифуркации. Установлено, что математическая модель полуколичественно описывает протекающие в системе процессы.

Исследования последних четырех десятилетий [1-3] показали, что во многих физических, химических, биологических и экологических системах возможно возникновение самоподдерживающихся упорядоченных структур. Это явление, характерное для протекания нелинейных неравновесных процессов в открытых системах при условии, что микропроцессы протекают кооперативно, получило название «эффект самоорганизации», а сами системы стали называться «диссипативными структурами» [3]. Возможность проявления пространственных, пространственно-временных и временных структур была признана после появления фундаментальных работ по неравновесной термодинамике, а также после введения понятий «диссипативная структура» и «самоорганизация» (И. Пригожин [1, 3]) и «синергетика» (Г. Хакен с сотр. [2]).

Особо следует выделить химические нестабильности - неравновесные фазовые переходы в химических системах, так как они дают много примеров реализации пространственных, пространственно-временных и временных структур [4]. В разных научных центрах проводят исследования по поиску новых колебательных химических систем, выявлению возможности появления колебательных режимов в ранее известных реакциях, а также разработке физико-математических основ описания наблюдаемых явлений [5].

Данная работа посвящена исследованию колебательных химических реакций, протекающих при окислении аскорбиновой кислоты в присутствии оксигенированных разно-лигандных комплексов железа(11) с диметилглиоксимом и гистидином в гомогенной среде, и рассмотрению общих подходов к интерпретации результатов такого типа исследований.

Экспериментальная часть

При проведении работы были использованы следующие реагенты: диметилглиоксим (ДМГ), L-гистидингид-рохлорид (His), аскорбиновая кислота (AH2), трис-(гид-роксиметил)аминометан (трис-буфер), гептагидрат сульфата железа(11) (FeSO4-7H2O), этанол-ректификат. Все растворы готовили растворением реактивов квалификации не ниже «ч.д.а.» в бидистиллированной воде, а в случае с ДМГ - в этаноле. Аскорбиновую кислоту очищали перекристаллизацией из воды, рН в системе устанавливали при помощи трис-буфера и контролировали при помощи иономера ЭВ-74. Катализатор в виде комплекса железа(11) с ДМГ и His готовили в соответствии с [6] при молярном соотношении Fe(II):ДМГ:His = 1:2:1, причем для каждой серии опытов использовали свежеприготовленный раствор катализатора. Исследование процесса проводили в стеклянном реакторе, состоящем из двух ячеек цилиндрической формы диаметром 37 мм и высотой 55 мм, соединенных между собой электролитическим ключом (насыщенный раствор KCl), и помещенном в ультратермостат (точность термостатирования составила ±0,01°), объем реакционной смеси соответствовал 20 см3. Регистрацию концентрационных колебаний в системе проводили путем измерения потенциала точечного платинового электрода относительно хлорсеребряного в течение определенного промежутка времени (интервал между условными единицами Ах = 12 с).

Результаты эксперимента

Типичная кривая зависимости относительного потенциала от времени приведена на рис. 1 . Окислительно-восстановительные процессы в системе аскорбиновая кислота-дегидроаскорбиновая кислота в присутствии

в качестве катализатора комплекса Ре(П)-ДМГ-Н1з-02 протекают в колебательном режиме.

Системы с регулярными колебаниями обычно имеют индукционный период (на рис. 1 не приведен). Важными характеристиками таких систем являются амплитуда и период (частота) колебаний [6].

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

Полученные данные показывают, что зависимость амплитуды колебаний от исходной концентрации реагента проходит через максимум (наибольшие значения амплитуды наблюдаются при Ск = 7,5-10-3 моль/л). Величина индукционного периода имеет тенденцию к повышению с ростом концентрации аскорбиновой кислоты. При этом не наблюдается однозначных зависимостей частот (периодов) колебаний от исходной концентрации субстрата.

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

Важным обстоятельством является то, что концентрационные колебания возникают только в присутствии ок-сигенированных комплексов железа(11). Тот факт, что именно связанный в оксигенированный комплекс кислород выступает в качестве окислителя, подтвержден нами экспериментально. В частности, использование солей железа(11) (Бе804) при значениях рН 7,96 приводит к тому, что эти ионы выпадают в осадок в виде гидроксосоедине-ний железа, а в отсутствие ионов металла окисление субстрата протекает в «гладком» кинетическом режиме.

Согласно литературным данным [8], катализатор имеет следующее строение:

0=0" А- О—Н

НС\ 0 Н'/ т " /СН3 С-NN-с

1С01<:О1

С-¥ / N-С

у/ъ-О+ъ-О

РеА

молекулярного кислорода к разнолигандному комплексу металла. Наряду с этим отмечено [8], что связанный в комплекс молекулярный кислород более реакционноспо-собен и процессы окисления с участием оксигенирован-ных комплексов протекают в «мягких» условиях.

Основные характеристики исследуемой колебательной системы в зависимости от температуры определяли при постоянных значениях концентрации аскорбиновой кислоты, катализатора и рН. Из полученных результатов следует, что с увеличением температуры величины амплитуды колебаний проходят через максимум, индукционный период в целом не меняется, а по основным частотам определенной зависимости проследить не удается. Важным обстоятельством является тот факт, что протекание реакции в колебательном режиме начинается при температуре не ниже 37°.

Таким образом, в результате эксперимента определены условия (пределы концентраций аскорбиновой кислоты Ск = 10-10-2 моль/л; катализатора Скат = 10-5- 10-4 моль/л; температуры в диапазоне Т = 3765°, рН 7,45-8,20), при которых имеет место протекание изучаемых процессов в колебательном режиме.

Обсуждение результатов

Анализ Фурье-преобразования временного ряда

В результате экспериментального изучения колебательных процессов обычно получают временной ряд в виде зависимости определенного аналитического сигнала от времени. Одной из основных задач в исследовании такого типа процессов является определение, какого рода временной эволюцией обусловлен этот сигнал. Существует несколько способов, позволяющих идентифицировать динамический режим и устанавливать его характеристики. Одним из приемов, использованных при выполнении данной работы, является метод, основанный на дискретном преобразовании Фурье [9].

Преобразование Фурье дискретного временного ряда х^ определяют [9-10] как операцию, порождающую соответствующий ему дискретный ряд х , такой, что

к

1 £ ( -2 —к

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

где к = 1, ..., п; - = -1

При интерпретации экспериментальных результатов

важно выяснить, какая информация о сигнале х(1) содер-

I I2

жится в спектре мощности, т.е. в зависимости х от

частоты V (V = кА1). На практике такая последовательность чисел конечна, содержит п значений и ограничена общей продолжительностью эксперимента 1тах = п-А1

Функция X, которая зависит от q независимых переменных 12, ...Дг, относится к периодической с периодом 2п по каждому из своих аргументов в случае, когда

к

значение ее не изменяется при увеличении любой из своих переменных на 2п. Это условие можно записать в виде:

Х(11Д2,...Д.р...Дг) = Х(у2,...;1 + 2я,... Дч), j =

Такая функция называется квазипериодической по времени, если все ее q переменные пропорциональны времени 1:

1 = од 1,

j =1,

Такого типа квазипериодическая функция имеет q ос-

новных частот V и такое же количество периодов Т.:

V = 2п,

j =1,

Т. = — = 2я/ю.

Спектр Фурье квазипериодической по времени функции в общем случае имеет довольно сложный вид. Поэтому при объяснении тех или иных экспериментальных фактов прибегают к определенным упрощениям. В частности, если квазипериодический сигнал х^!,...,«)^) является суммой периодических функций

х(ю11,.,юг1) = £ х. (соД

то, используя линейность соотношения (3), можно представить спектр мощности в виде суммы q спектров каждой из функций х^юД). Следовательно, суммарный спектр Фурье квазипериодической функции х(1), нелинейно зависящей от периодических функций переменных юД, будет состоять из множества пиков на частотах V,, V,. ^ и их гармониках, т.е. он содержит компоненты со всеми частотами вида I + т^2 + .+т^ I где т. - произвольные целые числа.

Исходя из этих соображений, спектры идентифицируют, пытаясь найти две основные частоты v1 и v2, которые позволят представить частоты линий с большой амплитудой в виде простых комбинаций 1т^ + т^2 I. Например, когда q = 2, каждая ненулевая компонента спектра сигнала х(ю11, ю21) является пиком с абсциссой |т^ + т^21 с небольшими т1 и т2 (0, ±1, ±2, ...).

Анализ такого типа был проведен для всех временных зависимостей, полученных в работе. Спектры мощности, соответствующие случаю, приведенному на рис. 1 , представлены на рис. 2. Выделенные в результате данного анализа основные частоты приведены в таблице. Получены три разных типа спектров Фурье (спектры с одной частотой и двумя частотами, где отношение V1/v2 иррационально и рационально). Так как идентификация спектров с одной частотой не представляет трудностей, рассмотрим анализ в последних двух случаях.

Как показано на рис. 2, а, в этом случае можно выделить три пика 1, 2 и 3 с частотами v1 = 2, v2 = 7 и v3 = 9 (в условных единицах). Переход к частотам в реальном времени проводится при помощи соотношения V. = 3600 ^/1024 т. Представив пик 3 как комбинацию частот V1 и v2 в виде |т^ + т^21, т.е. 2т1 + 7т2 = 9, получаем, что т1 = 1, т2 = 1. Поскольку отношение V1/v2 = 2/7 иррационально, то спектр сигнала с двумя периодами всюду плотен. Отсюда следует, что в спектре будет наблюдаться ограниченное число частот (в нашем случае равное двум), соответствующее линиям с большей амплитудой. Аналогичная картина наблюдается и для случая, представленного на рис. 2, б ^ = 4, v2 = 14 и V3 = 16). Если же отношение частот рационально (таблица), то спектр не везде плотен и квазипериодический сигнал в действительности является периодическим с периодом Т = п1Т1 = п2Т2, происходит синхронизация частот v1 и

ю

ДЕ, мВ 807060 50

200

400

600 т, усл. ед.

ДЕ, мВ 100 -

80

60 -

40 20

10 мин

200

I 1 г

400 600

т, усл. ед.

Рис. 1. Зависимость относительного потенциала процесса окисления аскорбиновой кислоты от времени: а - Ск = 7,5-10 моль/л, С = 2,5-10-5 моль/л, г = 40°, рН 7,96; б - С = 110--3 моль/л, С = 110-4 моль/л, г = 40°, рН 7,96

б

Основные характеристики колебаний при окислении аскорбиновой кислоты

Скат = 1-10-4 моль/л; t = 40°; рН 7,96

С^ 103,моль/л 0,5 1,0 2,5 5,0 7,5 10,0 15,0

Амплитуда, mV 18±2 25±3 30±3 38±4 45±5 25±3 15±2

Индукционный период, мин 6 8 12 15 25 22 10

Частота УЬ ч-1 V2, ч-1 У1/У2, усл.ед. 1,76 2,92 6/10 1,17 4,10 4/14 1,76 2,93 6/10 2,34 0,88 2,05 3/7 1,76 4,69 6/16 2,93 4,69 10/16

СR = 110-3 моль/л; t = 40°; рН 7,96

Си -105,моль/л 2,5 7,5 10,0

Амплитуда, mV 16±2 20±2 25±3

Индукционный период, мин 6 6 8

Частота V1, ч-1 V2, ч-1 Vl/V2, усл.ед. 1,76 2,93 6/10 1,17 4,10 4/14 1,17 2,34 4/8

СR = 1-10--3 моль/л; Скат = 1-10-4 моль/л, t = 40°

pН 7,45 7,96 8,20

Амплитуда, mV 8±2 25±3 26±3

Индукционный период, мин 8 8 6

Частота V1, ч-1 V2, ч-1 Vl/V2, усл.ед. 1,76 2,93 6/10 1,17 4,10 4/14 1,76 4,10 6/14

СR = 1-10-3 моль/л; Скат = 110-4 моль/л; рН 7,96

Ъ °С 37 40 50 60 65

Амплитуда, mV 15±5 25±3 25±3 15±2 12±2

Индукционный период, мин 9 8 10 10 8

Частота V1, ч-1 V2, ч-1 Vl/V2, усл.ед. 1,76 4,69 6/16 1,17 4,10 4/14 1,17 3,52 4,69 12/16 3,52 5,86 12/20

V2 или затягивание частоты. Все линии спектра Фурье при этом являются гармониками низшей частоты:

V = 1/г = ^/п1 = ^

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

Анализ полученных результатов позволяет сделать следующие выводы: 1) наблюдаемые осцилляции являются следствием протекания процесса окисления аскорбиновой кислоты в колебательном режиме, т.е. указывают на детерминированный характер периодических явлений;

следует отметить, что если бы эти колебания носили случайный характер, то спектр Фурье был бы сплошным; 2) в рассматриваемых системах число основных частот соответствует двум, следовательно, колебаниям подвергаются как минимум два компонента (исходные вещества, интермедиа-ты или продукты) реакционной смеси.

Вместе с тем Фурье-преобразование не позволяет произвести различие между динамическим хаосом и случайным сигналом. Это ограничение рассматриваемого подхода вынуждает обратиться к другим существующим методам, позволяющим охарактеризовать особенности динамики протекающих процессов [9].

20—

10 —

х

к

120

80

40

Л-

х

к

Ал лМ- „л-Л^лл

-АЛ-

50

100

V, усл. ед.

50

100

V, усл. ед.

Рис. 2. Спектры мощностей временных рядов: а - Ск = 7,510 моль/л, Скат = 2,5-10 моль/л, г = 40°, рН 7,96; б - Ск = 1-10 моль/л,

Ск

1-10 моль/л, г = 40°, рН 7,96

Восстановление аттрактора по временной последовательности данных

Важными характеристиками для описания динамики процессов, в которых наблюдаются автоколебания и вообще динамический хаос, являются величины размерностей фазового пространства и аттрактора [7]. Размерность фазового пространства п соответствует размерности химической системы (числу компонентов, реагирующих в системе) и поэтому по его размерности можно судить о минимальном числе переменных, необходимом для описания динамики. Аттрактор - это траектория в фазовом пространстве после окончания переходных процессов. Размерность аттрактора d указывает на тип колебаний: при d = 1 реализуются незатухающие периодические колебания, при d = 2 налицо квазипериодические колебания с двумя несоизмеримыми частотами, а в случае d > 2 и если при этом d является нецелым числом можно ожидать появления в системе хаотических колебаний. Последнее условие реализуется в большинстве случаев при d > 3. Следует также отметить, что значение размерности фазового пространства всегда больше величины размерности аттрактора (п > d).

Для описания динамики любой «сложной» системы необходимо знание ее фазового пространства. Фазовые портреты, являющиеся проекциями на определенные маломерные подпространства полного фазового пространства, колебательных химических реакций позволяют отразить их характерные особенности более четко, чем зависимости концентрации от времени. Обычно при построении

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

В данной работе в качестве такого параметра было использовано изменение относительного потенциала исследуемых систем от времени. При интерпретации полученных результатов нами был использован предложенный в [11, 12] и теоретически обоснованный в [13] эффективный метод оценки характеристик аттрактора в сочетании с возможностью восстановления траекторий в фазовом пространстве по задержкам времени. Использование данных подходов с целью описания динамики процессов позволило построить фазовые портреты и на их основе аттракторы с использованием диаграмм типа хп - х р где хп- значение функции х в определенный выбранный момент времени, а хп+1 - значение той же функции через

2

2

а

б

1

2

1

3

2

0

0

АЕ

82

78

74

70

70

74

78

82 АЕ

АЕ

92

84

76

68

72 76

84 88 92 АЕ

б

а

Рис. 3. Зависимости АЕ(г+г) от АЕ(г): а - Ск = 7,5-10 3 моль/л, Скат = 2,5-10 5 моль/л, г = 40°, рН 7,96; б - Ск = 1-10 3 моль/л,

Си = 1-10-4 моль/л, г = 40°, рН 7,96

определенный промежуток времени Аг. Исходя из одной, зависящей от времени, переменной можно восстановить траекторию в п-мерном фазовом пространстве, выбирая в качестве координат величины х^), х^+т), х^+2т),..., x(t+(n-1)x), где т - выбранная временная задержка. При выборе соответствующего значения т можно полагать [12], что эти переменные будут линейно независимы. В результате можно получить серию п-мерных векторов, представляющих фазовую диаграмму динамической системы. Существенным в этом подходе является то, что имеющаяся у исследователя информация достаточна для того, чтобы развернуть динамику системы в многомерном фазовом пространстве и отобразить фазовый портрет системы.

На рис. 3 приведены графики зависимости АЕ(г+т) от АЕ(т), где АЕ - изменение потенциала системы, г - время, а т - интервал времени между последовательными выборками для типичных случаев, представленных на рис. 1. Вид этих рисунков свидетельствует о сложном характере процессов, протекающих в исследуемых системах, а проявление фазового портрета в виде замкнутой кривой (предельного цикла) подтверждает заключение о реализации колебательных процессов.

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

Данных по такой обработке экспериментальных результатов слишком мало, чтобы на основе вида этих рисунков можно было сделать какие-либо заключения. С целью более четкой интерпретации экспериментальных результатов было проведено исследование динамики при помощи корреляционной функции аттрактора [11, 12]. Возможность построения корреляционной функции обо-

сновывается следующим образом: пусть X. (точка фазового пространства с координатами (х^)..., х0(^ + (п - 1)т)}) является началом отсчета для вычисления расстояния до остальных N - 1 точек (|X- X. |). Это позволяет сосчи-

VI , л /

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

С(г) =Nlim^0(1/N2).[число пар точек

для которых | Xi-Xj |< г],

где ^ ] - индексы, упорядочивающие точки вдоль траектории, на которой расположено N точек. Более строго функция С(г) определяется в виде:

1 N , ,

С

N у=1

где 0 - функция Хевисайда, равная по определению единице (при положительных х) и нулю (при остальных значениях х). Отклонение С(г) от нуля служит мерой влияния точки X. на положение других точек. Поэтому функцию С(г) рассматривают как интегральную корреляционную функцию аттрактора.

Если зафиксируем некоторое малое значение е и воспользуемся им в качестве своеобразного «метра» для зондирования структуры аттрактора, представляющего собой линию, то число пробных точек, расстояние которых до заданной точки не превышает г, пропорционально г/е [1]. В более общем случае, когда аттрактор представляет собой ^мерное многообразие, число точек будет пропорционально (г/е)а Поэтому при сравнительно малых г функция С(г) должна изменяться как

С

Иными словами, размерность аттрактора d определяется наклоном зависимости 1п С(г) от 1пг в определенном диапазоне г:

1п С

Эта функция называется корреляционной функцией аттрактора.

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

Для значений п = 2, 3, 4, 5, 6 и 7 были вычислены С(г) и определены значения размерности аттрактора d по тангенсам углов наклона касательных к кривым, построенным на основе этой корреляционной функции.

Использование полученных значений размерностей аттрактора позволило построить зависимости этих величин от размерности фазового пространства (рис. 5). Величины размерностей аттрактора d перестают зависеть от размерности фазового пространства п, как только п становится больше 5, а размерности аттрактора принимают значения, равные 3-4.

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

Математическое моделирование закономерностей протекания процессов в колебательном режиме

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

Кинетическая модель является основой математического моделирования механизма химических реакций. Кине-

тика сложной химической реакции, отвечающая определенной, заданной схеме превращений

vlгXl кг ; ^Х., 1, . = 1,., К; г = 1,., Я,

описывается системой дифференциальных уравнений вида:

ДГХ ] К V

^=£

где X. - вещества, Vlг - стехиометрические коэффициенты; Vlг = 0, если Х. не входит в число веществ, вступающих в реакцию на г-й стадии, vjг = 0, если веществ нет среди продуктов, образующихся на г-й стадии.

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

При моделировании кинетических закономерностей первостепенное значение приобретают как физико-химическое обоснование модели, так и исследование ее физической корректности [15]. В частности, если реакция осуществляется в закрытой системе, к математической модели предъявляются такие общие требования, как неизменность общей массы, невозможность существования отрицательных концентраций, выражение скорости реакции через концентрации по определенному кинетическому закону, существование положительной точки детального равновесия, стремление к равновесию с любого начального состояния, существование функции состояния системы, обладающей свойством диссипации. Основными же задачами анализа динамики неравновесных процессов являются исследование возможности существования стационарных состояний и определение их числа, определение устойчивости этих состояний с помощью какого-либо подхода (например, по Ляпунову), исследование проведения системы в различных пространствах (фазовом, динамическом, параметрическом) [14].

В разных условиях проведения процесса реакция окисления аскорбиновой кислоты в присутствии оксигениро-ванных комплексов железа(11) может обнаруживать разные типы поведения.

Для понимания химической основы этих сложных временных явлений как на качественном, так и на количественном уровнях на основе литературных данных [16, 25-27] был выполнен анализ возможных механизмов исследуемых процессов и из них выбран наиболее адекватно отражающий физическую картину исследуемого явления.

0.2 -0.4 -1 -1.6 -2.2

-2.8 -0

1п С

X. п=2 X п=3 X п=4 X п=5 X п=6 X п=7

0.4

1.6

2.2

2.8

1п ,г 3.4

0.5 0

-0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4

1п С

п=2 п=3 п=4 п=5 п=6 п=7

1п г

-0.2

0.4

1.6

2.2

2.8

3.4

1

Рис. 4. Зависимости 1пС(г) от 1пг: а - Ск = 7,5-10 3 моль/л, Скат = 2,5-10 5 моль/л, г = 40°, рН 7,96; б - Ск = 1-10-3 моль/л, Скат = 1-10-4 моль/л, г = 40°, рН 7,96

— а

0 4

а

4

4

2

2

п

0

0

0

4

8

8

Рис. 5. Зависимости d от п:

Ск = 7,5-10-3 моль/л, Скат = 2,5-10-5 моль/л , г = 40°, рН 7,96; б - Ск = 1-10-3 моль/л,

Ск

1-10 моль/л , г = 40°, рН 7,96

В качестве основных реакций, протекающих при окислении аскорбиновой кислоты, были использованы схемы из [16] с учетом того факта, что окислителем выступает молекулярный кислород, химически связанный со сме-шаннолигандным комплексом железа(11) с диметилглиок-симом и азотистым основанием. Выбор рассматриваемого ниже механизма, на наш взгляд, можно обосновать следующим образом: а) процесс проводили при рН 7,58,2; в этих условиях аскорбиновая кислота находится преимущественно в форме АН- [16]; б) экспериментально подтверждено, что именно связанный в оксигенирован-ный комплекс кислород выступает в качестве окислителя; в частности, как уже было отмечено, при использовании солей железа(11) (например, Ре804) при изучаемых значениях рН соединения железа в виде гидроксосоединений выпадают в осадок и в отсутствие комплексов металла окисление субстрата протекает в «гладком» кинетическом режиме; в) опытным путем подтверждено обратимое взаимодействие кислорода с рассматриваемого типа

координационными соединениями [6] (по этому вопросу имеются литературные данные [8]); г) образование частиц типа А , Н02- и Н02- описано в литературе [16]; д) математическому анализу подвергаются данные по концентрационным колебаниям, полученные в условиях реализации установившегося колебательного режима в исследуемой системе.

На основании изложенного схему процесса окисления аскорбиновой кислоты можно представить в виде:

1. АН- + Еек02 ^ Рек + А' + Н02-,

• - к2 2. А' + Н02 ^ А + Н02-

3. А + АН- ^ 2А + Н+

4. А + Н02- ^ А + Н02

а

к

5. Н02 + Н02^ Н202 + 02 ,

6. Ре, + + 02 ^ Рек0,2

Необходимыми условиями реализации колебательного режима являются наличие автокаталитических процессов (стадии 2 и 3), а также регуляция по типу обратной связи (стадии 2, 4). Тогда дифференциальные уравнения, описывающие кинетические закономерности изменения концентраций веществ в реакционной смеси, имеют вид:

а[АН-]/Л = -к1[АН-[Рек022+] - к3[А][АН-],

^0/+]/^ = -к1[АН-][Еек022+] + к6[Рек2+][02],

d[Fek2+]/dt = к^АЩ^е^] - к6[Рек2+][02],

d[AЛ]/dt =к1[АН-][Рек022+] - к2[Ал][Н02] + + 2к3[А][АН-] + к4[А][Н02-],

d[H02-]/dt = к1[АН-][Рек022+] - к2[Ал][Н02] + + к4[А][ Н02-] - к5[Н02-]2,

d[A]/dt = к2[Ал][Н02] - к3[А][АН-] - к4[А][ Н02-], d[H02-]/dt = к2[Ал][Н02] - к4[А][Н02-], d[H202]/dt = к5[Н02-]2 , d[02]/dt = к5[Н02-]2 - к6[Рек2+][02].

Из соотношений (12) видно, что в качестве математической модели получается система из девяти автономных обыкновенных дифференциальных уравнений нелинейного типа.

Прежде чем провести качественный анализ данной системы нелинейных дифференциальных уравнений и численное ее решение, целесообразно провести термодинамический анализ процессов, протекающих по схемам (1-6), так как он позволяет указать на возможность и причины возникновения химических неустойчивостей.

Термодинамический анализ

О термодинамической устойчивости системы вдали от равновесия судят по знаку производной второй вариации энтропии (дЭ^р528 (р - средняя плотность), которая для стационарного состояния является функцией Ляпунова [17]. Если эта величина для стационарного состояния имеет положительный знак, то состояние системы устойчиво, а в случае, когда она принимает отрицательное значение, - неустойчиво, и система может перейти в качественно новое состояние. С целью определения причин и возможности потери устойчивости состояния системы и перехода в новое качественное состояние на основе функции Ляпунова был проведен термодинамический анализ протекающих процессов.

Если рассмотреть реактор, в котором протекают необратимые химические реакции, то аналитическое выражение для производной второй вариации энтропии имеет вид [18 ]:

í

где цк - химический потенциал, хк - концентрация к-го компонента, V11 - объемный расход реагентов, q1 - тепловой поток через поверхность Ат - химическое сродство, Т - температура, V - объем реактора.

Первый член, характеризующий изменение производства энтропии, обусловлен обменом энергии с окружающей средой, он равен

/

где р - средняя плотность, СТ - теплоемкость, КТ - коэффициент теплопередачи.

Изменение избытка энтропии, характеризующее обмен массы с окружающей средой, соответствует:

{5

где х0 - концентрация реагентов на входе в реактор.

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

Х1 + Х2 + ... ^ У; + У2 + ...

то выражения для скорости и химического сродства имеют вид:

= kiПxi ; к = к0е

-Е /КГ

А" = КТ1пКр + КТ1п

где к- и К - константы скорости и равновесия реакции, а х1 и yJ - концентрации веществ Х1 и YJ соответственно. Учитывая, что

51пК 1 = -Ю/КГ2 )8Т,

где - тепловой эффект реакции, вариация химического сродства, отнесенного к температуре, соответствует

5

вариация же скорости реакции запишется в виде:

5\¥1 = "1

С учетом этих замечаний и рассматриваемых в работе кинетических схем было построено выражение для термодинамической функции Ляпунова.

Введем следующие обозначения: [АН ]исх = а; [Рек022+] = Ь; [Ал] = х; [Н02] = у; [А] = 7; [АН-] = а-х-ъ = И. Тогда вариации химического сродства, отнесенного к температуре, для каждой стадии схем (1-6) будут равны:

8(А^/Т) = -Д/Т^Т + К(8И/И) - К(8у/у) + К(8х/х);

8(А"2/Т) = /Т2)8Т + К(8х/х) + К(8у/у) - К(8ъ/ъ);

8(А"3/Т) = /Т2)8Т + К(8И/И) + К(8ъ/ъ) - 2К(8х/х);

S(AW4/T) = —(Q4 /T2)ST + R(Sz/z) - R(Sx/x) - R(Sy/y); S(Aw5/T) = -(Q5 /T2)ST + 2R(Sy/y);

S

а вариации скоростей реакции соответствуют:

Sw1 = w1( Sh/h) + (E1/RT2)w1ST; Sw2 = w2(Sx/x) + + w2( Sy/y) + (E2/RT2)w2ST;

Sw3 = w3(Sh/h) + w3(Sz/z) + (E3/RT2)w3ST; Sw4 = = W4(Sz/z) + (E4/RTV4ST;

Sw5 = 2w5

Тогда производная термодинамической функции Ляпунова имеет вид:

ÎSwtS(Awi/T )dV = V{[(Sx/x)2 Rw2 + (Sy/y)2 R(w2 +

+ 4w5) + (S z/z)2 R(w3 + w4) + (S h/h)2 R(w1 + w3) +

+ (Sx/x)(Sy/y)2Rw2 + (Sx/x)(Sz/z)R(-w2 - 2w3 - w4) +

+ (Sx/x)(Sh/h)R(-w1 - 2w3) + (Sy/y)(Sz/z)R(-w2 - w4) +

+(Sy/y)(Sh/h)R(-w1) + (Sz/z)(Sh/h)2Rw3] +(S T/T2).

.(Sx/x)[-WiEi - W2(Q2 - E2) -2W3E3 - W4E4)] +(ST/T2).

.(Sy/y)[-W1E1 - w2(Q2 - E2) - W4E4 - 2w5 (Q5 + E5)] +

+ (S T/T2)(S z/z)[-W2E2 - W3(Q3 - E3) - W4 (Q4 - E4)] +

+ (ST/T2)(Sh/h)[ -w1(Q1 - E1 - w3(Q3 - E3)] + [(ST2/T).

.

Анализ полученного выражения позволяет указать на причины, вызывающие потерю состояния устойчивости и побуждающие к переходу ее в иное качественное состояние. В соотношении (24) слагаемые, отвечающие за протекание прямых реакций, являются положительно определенными квадратичными формами, и это обстоятельство способствуют стабилизации системы. Однако наличие автокаталитических процессов (стадии 2, 3) приводит к тому, что в выражении производной функции Ляпунова появляются отрицательные слагаемые

-[(Sx/x)(Sz/z)R(w2 + 2w3) + (Sx/x)(Sh/h)Rw2 + + 2(Sz/z)(Sh/h)Rw3].

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

Необходимо также рассмотреть вклад в производную второй вариации энтропии слагаемых, связанных с наличием обратных связей. В выражении (24) за это отвечают слагаемые -[(Sx/x)(Sz/z)R(w2 + w4) + (Sy/y)(Sz/z)R(w2 + w4)], имеющие отрицательный знак. При преобладании вклада данных членов также вероятно возникновение химических неустойчивостей.

Рассмотрим слагаемые в производной термодинамической функции Ляпунова, отвечающие за возникновение

термокинетических колебаний, связанных с автокатализом: -(8Т/Т2)(8х/х)^2^2 - Е2) + 2-№3Б3)]. Если повышение температуры будет способствовать увеличению скорости стадий (2) и (3), то в системе будет происходить накопление анион-радикала А-'- (5Т > 0 ^ 5х > 0). Следовательно, рассматриваемые члены будут принимать более отрицательные значения, вследствие чего будет усиливаться дестабилизация системы. Наличие обратных связей также может привести систему к термокинетической нестабильности. В реализации обратных связей принимают участие члены

-(5T/T2)(5x/x)[w2(Q2 - Е2) + 2w3E3)] - (5T/T2)(5z/z)[w2E2 +

+ W4 «4 - Е4)].

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

Анализ математической модели

Большая часть данных по изучению колебательных химических реакций, экспериментальному определению условий их реализации и определению параметров концентрационных колебаний получена в процессе численного эксперимента [19]. Как и во многих направлениях современных теоретических исследований, в области описания диссипативных структур сложилось фактически два взаимодополняющих подхода: формальный математический анализ на основе абстрактных понятий современной математики, осуществляемый на том же уровне строгости, что и собственно математические исследования, но ориентированный на специфические особенности изучаемых объектов, и численное исследование сравнительно простых моделей, демонстрирующих ключевые особенности протекающих процессов.

Одним из подходов к объяснению динамических особенностей является исследование поведения системы в фазовом пространстве. Построение фазовых портретов систем позволяет представить решение в целом и определить как характер стремления системы к стационарному состоянию и тип особой точки [14], так и области притяжения устойчивых периодических решений, отвечающих предельным циклам. При этом получение решения системы дифференциальных уравнений в фазовом пространстве в виде предельных циклов указывает на возможность возникновения автоколебаний. Для ответа на эти вопросы необходимо использование методов современной качественной теории дифференциальных уравнений и решение их численными методами.

Для описания данной схемы реакций необходимо учитывать изменение во времени концентрации девяти компонентов (математическая модель представляется в виде

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

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

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

Так как нам известны исходные концентрации реагента и катализатора, то через уравнения материального баланса можно выразить концентрации веществ, в которые они превращаются согласно приведенной кинетической схеме. Ограничением при составлении математической модели было и то обстоятельство, что концентрация молекулярного кислорода в ходе реакции не изменяется (она определяется растворимостью кислорода в воде и составляет примерно 10-3 моль/л [20]); постоянство концентрации частиц НО2- определяется обратимостью реакций их образования (реакции 2 и 4 в кинетической схеме), т.е. для рассматриваемого изменения концентрации частиц НО2-можно использовать принцип стационарности Боденштей-на, в соответствии с ним ^Н02~]/Л = 0 и, следовательно, [НО2] постоянна.

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

¿Сх/Л = к1Св(СА - сх - с2) - к2СхС¥ + 2к3С2(СА - сх -

- Сг) + к4^

аСу/Л = ЦСв(СА - Сх- С2) - к2СхС¥ + к4С2 - 2к5Су2,

= к2СхСУ - к3Сг

где Сх = [А-], Су = [Н02 ], С = [А], СА = СК(исх),

Св = Си(исх).

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

Стационарные состояния определяются из системы алгебраических уравнений [21]

дх, у, г) = 0,

а тип устойчивости определяется корнями уравнения

X3 + А^2 + А2Х + А3 = 0,

где А1, А2, А3 вычисляются через элементы матрицы Якоби системы в стационарном состоянии. Нас будут интересовать случаи, при которых стационарное состояние неустойчиво. С этой целью в соответствии с существующей процедурой [15]:

1. Проводят линеаризацию системы. Для этого: а) вводят новые переменные £ = с - с, характеризующие отклонения текущей концентрации с от концентрации с в особой точке; б) в уравнении (26) от х переходят к £; в) проводят разложение правой части в ряды по степеням £ и пренебрегают нелинейными членами. Полученную линеаризованную систему в векторном виде можно пред-ствить следующим образом:

Л

= А(£),

где А - матрица частных производных; А =

ЭДе,К) до .

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

2. Получают характеристическое уравнение, играющее важную роль при анализе устойчивости в малом,

¡А - ёА| = 0,

где Е - единичная матрица.

3. Исследуют корни характеристического уравнения. Изучение устойчивости особой точки в малом сводится к определению знака корней алгебраического уравнения п-й степени. В соответствии с теоремой А.М. Ляпунова тривиальное решение систем обыкновенных дифференциальных уравнений неустойчиво, если хотя бы один корень характеристического уравнения имеет положительную вещественную часть [21].

При анализе сложных реакций, когда уравнений три и более, исчерпывающий качественный анализ затруднителен [14]. Схемы сложных реакций, приводящие к трем уравнениям, обсуждались в различных работах [5], однако системы рассматриваемого типа ранее не исследовались. Предпринята попытка нахождения области параметров обсуждаемой модели, где существуют автоколебания. Для решения поставленной задачи система вначале была приведена к безразмерному виду с помощью соотношений:

x = юСх, y = nCY, z = yCZ, t = St', Z = k3k4/kjk2, M = k^/k^2, ю = k^^2, n = к^/Ц2^2, Y = k3/kj, S = p = k42CB/k2k3CA, о = k2k3/2k5k4.

В безразмерном виде система уравнений выглядит так: Zdx/dt =

dz/dt = pxy - (a - Zx - z)z - p-1z,

а стационарная точка является единственной с координатами

X = a/Z, Y = 0, Z = 0.

Характеристический многочлен имеет вид:

Zb - х Zb

0

(p/Z)a

(p/Z)a - X

-(p/Z)a

b - p-1 b - p-1

- p-1- X

= 0.

Коэффициенты характеристического уравнения (35) равны:

Л! = -[р-1 + ар/ж + жЬ]; А2 = [аЬр/ж + ажр-1 + 2аЬр];

А3 = -2аЬ2р, или Д1 = А1 < 0; Д2 = Л1Л2 - Л3 < 0; Д3 =

= АзД2 > 0.

Для определения устойчивости особой точки воспользуемся критерием Раусса-Гурвица [21]. Для этого составим матрицу из коэффициентов характеристического уравнения

i А 1 0 >

A3 А2 А1

0 V 0 А3

Условием асимптотической устойчивости является положительность всех главных миноров матрицы (33). В обсуждаемом случае Д1 = А1 < 0; Д2 = Л1Л2 - А3 < 0; Д3 = АзД2 > 0, т.е. особая точка является неустойчивой. Для более детального анализа характера неподвижной точки (стационарного состояния) в случае системы из трех дифференциальных уравнений необходимо знать выражение [18]:

Q = -A12A22 + 4A13A3 + 4A3 - 18A1A2 A3 + 27 A32.

В нашем случае О > 0, А2 > 0, А3 < 0, Л1Л2 - Л3< 0, следовательно, действительные части корней характеристического уравнения положительны и поэтому реализуется неподвижная точка типа неустойчивый фокус.

Таким образом, проведенный качественный анализ системы обыкновенных дифференциальных уравнений пока-

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

Вторым этапом исследования было проведение численного эксперимента.

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

Исследуемая нами модель сформирована в виде трех взаимосвязанных нелинейных дифференциальных уравнений. Компьютерное представление частного решения этих уравнений получается после приближенного интегрирования системы с использованием определенного численного алгоритма. При выполнении настоящей работы для реализации решений систем обыкновенных дифференциальных уравнений была использована программа МаЛеай с подпрограммой Rkadapt, внутренне адаптированной к шагу интегрирования [22]. Достоинства этой программы по сравнению с существующими заключаются в том, что имеет место хорошее взаимодействие с аналитическими методами, развитыми в качественной теории дифференциальных уравнений, разумное соответствие физическим явлениям в системах с сосредоточенными параметрами, а также ясность и наглядность получаемых результатов в терминах траекторий динамической системы и ее фазовых портретов [23, 24].

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

0,04

X

у.

ъ

0,01

110

2000 4000 6000

-4

2000 4000 6000

Рис. 6. Зависимость изменения концентрации частиц [А\ ], [Н02\] и [А] от времени

X.

у.

ъ.

0,01

0,005

1 000

1 000 2000 ti .

2000

1 .1 0-4

1 000 2000 ti

Рис. 8. Зависимость изменения концентрации частиц [А\ ], [Н02.] и [А] от времени

0

0

-4

4-10

5-10

-4

2-10

2-10

0

t

констант скоростей брутто-реакций (например, окисления аскорбиновой кислоты в дегидроаскорбиновую).

При выполнении настоящей работы значения констант скоростей были заимствованы из [16, 25-27]. При использовании этих значений констант в оценке параметров ц и р величина каждого из них составила 0,1. На основе этих данных была предпринята попытка определения пределов концентрации реагента и катализатора, при которых возникают концентрационные колебания.

С этой целью при численной реализации математической модели были определены условия: а = 10-3-10-2; Ь = 10-4-10-3; число точек принято равным 5000, а величина шага соответствует 1 . Вместе с тем при получении решений учитывается и тот факт, что в большинстве случаев подвергаются колебаниям концентрации как минимум двух веществ (как было отмечено ранее, на это указывают результаты Фурье-анализа экспериментальных данных). Кривые в координатах «концентрация реагентов (х = [Ал]; у = [Н02]; г = [Л]) - время и фазовые портреты» представлены на рис. 6, 7. В указанных условиях проведения численного эксперимента реализуются незатухающие колебания,

а фазовые портреты имеют вид предельного цикла.

При изменении параметров а и Ь вид кривых соответствует затухающим или возрастающим колебаниям. Например, когда а = 10-2 и Ь = 5-10-4, имеет место проявление затухающих колебаний (рис. 8, 9).

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

Вместе с тем полученные численные результаты находятся в удовлетворительном согласии с экспериментальными данными и области реализации концентрационных колебаний соответствует СА = 10 -10 , Св = 10-4-10-2 (моль/л).

В заключение отметим, что исследована ранее не описанная каталитическая система, в которой реализуются концентрационные колебания.

СПИСОК ЛИТЕРАТУРЫ

1. Николис Г., Пригожин Н. Познание сложного. М., 1990.

2. Хакен Г. Синергетика. Иерархии неустойчивостей в самоорга-

низующихся системах и устройствах. М., 1985.

3. Николис Г., Пригожин Н. Самоорганизация в неравновесных

системах. М., 1979.

4. Жаботинский А.М. Концентрационные колебания. М., 1974.

5. Колебания и бегущие волны в химических системах /Под ред.

Р. Филда, М. Бургер. М., 1988.

6. Магомедбеков У.Г. // ЖНХ. 1977. 42. С. 277.

7. Яцимирский К.Б. // Теорет. и экспер. химия. 1988. 24. С. 488.

8. Братушко Ю.И. Координационные соединения 3d переходных

металлов с молекулярным кислородом. Киев, 1 987.

9. Берже П., Помо И., Видаль К. Порядок в хаосе. О детермини-

стическом подходе к турбулентности. М., 1991.

10. Отнес Р., Энонсон Л. Прикладной анализ временных рядов. Основные методы. М., 1 982.

11. Grassberger P., Procaccia I. // Physica. 1983. 9D. P. 189.

12. Ben-Mizrachi A., Procaccia I., Grassberger P. // Phys. Rev. 29A. P. 975.

13. Takens F. Lecture Notes in Math. Heidelberg - N.Y., 1981. P. 898.

14. Быков В.И. Моделирование критических явлений в химической кинетике. М., 1988.

15. Химическая и биологическая кинетика / Под ред. Н.М. Эмануэля, И.В. Березина, С.Д. Варфоломеева. М., 1983.

16. Mushran S.P., Agrawal M.C. // S. Scient. Ind. Res. 1977. 36. P. 274.

17. Кафаров В.В., Дорохов И.Н., Кольцова Э.М. Системный анализ в химической технологии. Энтропийный и вариационный методы неравновесной термодинамики. М., 1988.

18. Кольцова Э.М., Гордеев Л.С. Методы синергетики в химии и химической технологии. М., 1999.

19. Ермолаев Ю.Л., Санин А.Л. Электронная синергетика. Л., 1989.

20. Досон Р., Эллиот Д., Эллиот У., Джонс К. Справочник биохи-

мика. 1991. С. 476.

21. Баутин Н.Н. Поведение химических систем вблизи границ области устойчивости. М., 1 984.

22. Pojman J.A. Studying nonlinear chemical dynamics with numerical experiments. Department of chemistry & Biochemistry. University of Southern Mississippi, 1997. P. 339.

23. Херхагер М., Парголль Х. Mathcad 2000. Полное руководство. Киев, 2000.

24. Strizhak P., Menzinger M. // J. Chem. Ed. 1996. 73. P. 868.

25. Сычев А.Я., Исак В.Г. // Успехи химии. 1995. 64. С. 1183.

26. Grinstead R.R. // J. Am. Chem. Soc. 1960. 82. P. 3464.

27. Khan МЖТ, Martell A.E. // J. Am. Chem. Soc. 1967. 89. P. 7104.

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