О двумерной модели реакции окисления оксида углерода на наночастицах палладия Текст научной статьи по специальности «Математика»

УДК 517.928.4:517.929.5

О ДВУМЕРНОЙ КИНЕТИЧЕСКОЙ МОДЕЛИ РЕАКЦИИ ОКИСЛЕНИЯ ОКСИДА УГЛЕРОДА НА НАНОЧАСТИЦАХ ПАЛЛАДИЯ*)

Е, А. Лашина, Н, А. Чумакова, Г, А. Чумаков

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

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

Работа выполнена при поддержке Сибирского отделения РАН (Междисциплинарный проект 80).

@ 2013 Лашина Е. А., Чумакова Н. А., Чумаков Г. А.

1. Кинетическая модель

В условиях автоколебаний скорости реакции окисления СО палладий на поверхности катализатора присутствует в двух фазах: металлической и оксидной. На поверхности металлического палладия присутствуют адсорбированные соединения СО[Рс1] и 0[Рс1], а также свободные активные центры [Рс1]; на поверхности оксидной фазы — СО[Р(Ю] и свободные активные центры [РсЮ]. Предполагается, что общее число активных центров палладия неизменно в ходе реакции, хотя в условиях автоколебаний соотношение числа центров типа [Рс1] и [РсЮ] меняется.

Для описания экспериментальных данных в работе [2] предложен механизм реакции окисления СО, состоящий из 6 стадий, причем стадии 1 и 4 обратимы:

1) СО(ё) + [Рс1]^СО[Рс1],

2) 02(ё) + 2[Рс1]^20[Рс1],

3) СО[Рс1] + 0[Р(1] ^ С02(ё) + 2[Р(1],

4) СО(ё) + [РсЮ] ^ СО [РсЮ],

5) СО [РсЮ] ^ С02(ё) + [Рс1],

6) 0[Рс1] ^ [РсЮ].

Здесь использованы следующие обозначения: [Рс1] и [РсЮ] — активный центр на поверхности металлического и оксидного палладия; СО[Рс1] и 0[Рс1] — оксид углерода и кислород, адсорбированные на металлической поверхности; СО[РсЮ] — оксид углерода, адсорбированный на оксиде палладия; СО^), 02^) и С02^) — оксид углерода, кислород и диоксид углерода в газовой фазе над поверхностью катализатора.

На основании данного механизма и закона действующих поверхностей построим кинетическую модель, описывающую динамику безразмерных концентраций х, у и V соединений СО[Рс1], 0[Рс1] и СО [РсЮ]. Отметим, что х, у и V определяются как доли активной поверхности, запятой соответствующим веществом, поэтому x,y,v £ [0,1]. Дополнительно предполагается, что доля £ металлической фазы [Рс1] па нано-частицах палладия изменяется в ходе реакции [2], так что £ = 0 соответствует случаю, когда вся поверхность палладия является оксидной, а £ = 1 — металлической. Поскольку в условиях автоколебаний на по-

верхности наночастиц палладия присутствуют обе фазы [Рс1] и [РсЮ], далее будем рассматривать 0 < £ < 1.

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

х ■

х = к1(1-х-у) - к-\(у)х - к3(у)ху£ -У = ЫуХ1 - х- у)2£ - к3(у)ху£ - аеу -

£ (1)

£ = ф(1 - £) - ау£),

V

V = /С4(1 — у) — к-4« — £У н---£,

1 - £

где к > О — так называемые константы скоростей стадий 1-6, г = ±1,2,3, ±4,5,6, £ = и а = къ/к§, точками обозначены производные по времени. Индекс г > 0 соответствует «прямой» реакции (стрелка

г<

(стадии 1 и 4).

В модели (1) предполагается, что активность катализатора определяется составом его поверхности, и для ] = — 1,2,3 константы скоростей стадий к^ = к(у) > 0 не являются постоянными, а зависят от у

к,-(у) = к™ + Щ ■ в(у), 3 = — ,2,3, (2)

где постоянные Дк^ и к^ положительны. Фу нкция в (у) непрерывно дифференцируема при у £ [0,1] и имеет вид

1, 0 < у < Ус — 5, '(у) = { «(у)> Ус — 5 <у <ус + 5, (3)

— 1, ус + 5 < у < 1,

при этом в (у) монотонна, так что ^^^ = 0 и < 0 при ус — 5 < у < ус + 5. Параметры 5 £ (0,0.5) и ус £ (5,1 — 5) определяют ширину

и расположение интервала значений у, на котором в (у) изменяется от —

ности металлической фазы палладия.

Правая часть системы (1) относительно вектора переменных (ж, y, v, f) рассматривается в области П = G х Iv х Ц, где

G= {(x,y) G R2:x >0, y >0, x + y < 1},

(4)

Iv = {v G R : 0 < v < 1}, I5 = {f G R : 0 < f < 1}.

Отметим, что согласно экспериментальным данным образование оксида палладия (стадия 6) и его восстановление (стадия 5) происходят более медленно по сравнению с остальными стадиями реакции. Поэтому кинетическая модель (1) содержит малый параметр 0 < е

Изучение системы (1) проводилось методами качественной теории обыкновенных дифференциальных уравнений и вычислительной математики [3—5]. Численное исследование выполнено в случае, когда

/ п \ kj - kj kr + kj

g(y) = Sm(-(yc-y)), А % = = (5)

где kj > kr > 0, j = — 1 , 2, 3, и принимают следующие значения [2]:

k^= 49.82, k^ = 0.19, kl2 = 105, kj = 2 • 104,

, к " (6)

kl3 = LW аО5, kr3 = í1.8, Ус = 0.37, ¿ = 0.1. При этом из (2), (3) и (5) следует, что kj(0) = kj и k¿(l) = kr.

2. Вырожденная система

е

(1). С точки зрения химии это условие означает, что на оксидной фазе [PdO] не происходит химического превращения (стадия 5), а возможна только адсорбция-десорбция СО (стадия 4). Активной для химических превращений веществ остается только металлическая фаза палладия (стадии 1-3). е

X = ki(l — x — y) — k_i (y)x — k3(y)xyf, y = k2(y)(l — x — yff — k3(y)xyf, (7)

f = 0

и уравнение

v k — v — k- v.

Общее решение линейного уравнения (8) имеет вид

v(t) = е-а4(щ — Ь) + Ь, а = к± + к_4, Ь = к±/а, 0 ^ щ ^ 1.

Поскольку из условия на параметры модели (1) следует, что а > 0 и О < Ь < 1, то для любых начальных данных ^^ = Щ £ 1у функция v(t) монотонна, v(t) ^ ^и £ ^ и v(t) = vs = Ь является асимптотически устойчивым положением равновесия уравнения (8). Отметим, что решение ^^^ ^^^адает в --окрестность значения vs за время

к4 + к_4 -

Следовательно, при больших к± + к_4 наблюдается быстрая стабилизация компоненты v(t) к стационарному значению, состояние поверхности оксидного палладия становится близким к равновесному.

Уравнение £ = 0 имеет постоянное решение £(£) = £о, £о £ (0,1), которое соответствует постоянству соотношения фаз [Рс1] и [РсЮ] на поверхности катализатора.

3. Однопараметрическое семейство двумерных динамических систем

Рассмотрим семейство динамических систем относительно переменных х и у в области О с параметром 0 < £о ^ 1:

х = к1{\ — х — у) — к— (у)х — к3(у)ху£0 = Р(х,у,£0), ^

у = к2(у)(1 — х — у)2£0 — к3(у)ху£0 = £0^(х,у).

О

т. е. если (хо,уо) £ О, то решение системы (9) с начальными данными х(0) = хо и у(0) = уо определено для всех 0 ^ £ < + то и (х(£), у(£)) £ О.

Заметим, что для доказательства этого свойства без ограничения общности можно положить £о = 1 и обозначить Р(х, у, 1) = Р(х,у). Векторное поле системы (9) па границе дО области О всюду за , , О

касается дО в точке (1,0): Р(1,0) = — к^ < 0 и ф(1,0) = 0. При всех

значениях параметров Р(0,1) = ф(0,1) = 0, т. е. (х8, у8) = (0,1) стационарное состояние системы (1). Матрица Якоби

Рх(0, 1) Ру(0,1) - (кг + к— + кг3) -кг

А \Зх(о,1) 3у(0,1)1 V -кг о

_-(к1 + кг_1 + ц)±л/щтжТ+Ц?

Л1 2 — -

имеет характеристические корни

2

и А1 < 0, > 0. Откуда следует, что стационарное состояние (0,1) системы (9) является седлом. Собственные векторы (01,61) и (02,62) матрицы А, соответствующие А1 и А2, таковы, что Ъ\/а\ = кГ| > 0 и 62/12 = —кГ/А < —1- Следовательно, в малой окрестности неподвижной точки (0,1) одна го выходящих сепаратрис лежит в С, а входящие сепаратрисы не принадлежат области С.

Лемма 1. Существует область Сд С С, 0 < 9 С 1, такая, что граница дСд области Сд является циклом без контакта для системы (9),

т. е. во всех точках множества дСд векторное поле системы (9) наирав-Сд

Доказательство. Рассмотрим точку (1,0). Так как Р(1,0) < 0, существуют 0 < 9\ С 1 и множество

идх= {(х,у) : 1 - 91 < х < 1,0 < у < 9ьх + у < 1}

такие, что Р(х,у) < 0 для всех (х,у) € С/д1. Кроме того, поскольку Ф(х,0) > 0 для 0 ^ х < 1 и ^(х,у) < 0 для х + у=1и0<х<1, векторное поле па границе дид щи х = 1 - 91 направлено внутрь области С\ид1.

Рассмотрим главные изоклины системы (9), т. е. множества точек области С, в которых Р(х, у) = 0 или ф(х, у) = 0. Пусть (х,у) € Си Р х, у

"*'<»>-1ГТШтщ-« (ю)

Уравнение ф(х, у) = 0 при каждом у € [0,1) является квадратным х

корня

XI,2Ы = 1 - У + А(у)/2 ± у/А(у)(1-у) + А2(у)/4,

где Л(у) = кз {у)у/к2 (у) ■ Очевидно, что х ,2(0) = 1, а при уф 0 имеем х(у)+у > 1 и (ху) £ О, в то время как х2{у)+у < 1 и (х(у), у) £ О для всех у £ [0,1]. Следовательно, вторая главная изоклина определяется функцией хц(у), где

X = хд(у) = 1 -у + А(у)/2 - у/А(у)( 1 - у) + А1(у)/4, у £[0,1]. (11)

Рассмотрим функции хР(у) и хц(у), задаваемые в (10), (11). Так

как

хР(1) = ха(1) = о, ^(1) = о, ^(1)<о, ау ау

существует 0 < 62 ^ 1 такое, что хР(у) > хц(у) для всех у £ [1 — 62,1)-Выберем 63 = 63(62) £ (хц(1—62), хР(1—62)) такое, что Р(#з , 1 — > 0, Р(62,1 — 62) < 0 и ф(х, 1 — 62) < 0 для 63 < х < #2- Отметим, что поскольку графики функций хР (у) и хц(у) для у £ [0,1] лежат в области О, то 63 < 62 и 63 ^ 0 при 62 ^ 0. Кроме того, так как Р(0, у) > 0 для 0 ^ у < 1, то 62 и 63 можно выбрать так, что Р(63,у) > 0 для 0 < у < 1 — 62.

Пусть Од = О\{ид1 и П^ и ид,2}, где 6 = тах 6^ и

¿=1,2

Пд2 = {(х,у): 0 < х < 63(62), 0 < у О — 62}, {(х, у) : 0 < х < 1 — 62 < у < 1, х + у < 1}, дОд Од

Од

Од

{х$,у$) £ Од, то решение (х(1),у(1)) задачи Кошн для системы (9) с начальными данными х(0) = хо, у(0) = уо существует, единственно н (х(Ь), у(Ь)) £ Од для всех 0 < I < + то.

6>

справедливо для всех 0 < 6 < 60 н дОд ^ дО при 6 ^0.

Од

На основании теоремы существования и единственности решения задачи Коши, теоремы о непрерывной зависимости решения от начальных данных и следствий 1 и 2 из леммы 1 справедлива следующая

Теорема 1. Пусть (x0,yo) G G. Тогда решение x = x(t), y = y(t) системы (9), удовлетворяющее условию x(0) = щ, y(0) = yo, существует, единственно н (x(t), y(t)) G G для всех 0 ^ t < + ж.

4. Анализ положений равновесия при 0 < £0 < 1

Рассмотрим главные изоклины системы (9) для 0 < £о < 1, а именно множества точек (x,y) G G, в которых P(x,y,£o) = 0 либо Q(x, y) = 0. Выражая x из уравнения P(x, y, £о) = 0 как функцию от y и получим формулу, аналогичную (10):

х = хР(у,£о) = -—, , klj~\ , У} , ч е • (12)

ki + k— (y) + h(y)yb

Отсюда следует, что для любого 0 < £о < 1 изоклина (12) лежит в области G, причем значения xP( xP(0, £о) = ki/(ki + k_i(0)) < 1 не зависят от

Вторая изоклина Q(x,y) = 0 не зависит от £0- Представив это уравнение как (1-х — у)2 — (VÄ ■ л/х)2 = 0, получим формулу, аналогичную (12):

xq{y) = /Ы(1 - y), (13)

где

/Ы=уЦ(1-*) + ЛЫ-уЗД. (14)

у/4(1-у) + А(у)+у/Щ Отметим, что /{y) G (0,1) при y G (0,1), ^^ 1 и /(1) = 0. Тогда 0 < xq xq(0) = 1 и xq(1) = 0.

Стационарные состояния (xs,ys) системы (9) являются точками пересечения главных изоклин, т. е. решениями системы уравнений P(xs,ys,£o) = 0, Q(xs,ys) = 0. Тогда выполнены равенства

xs = xQ{ys), xs = xp (ys,b), (15)

т. e. xq{ys) = xp(ys,£o)-

Индекс особой точки [3]. Пусть C — простая замкнутая крн-вая, v — заданное на ней векторное поле, ad — некоторая прямая на x, y

точек М^ (к = 1,2,... , п) кривой С, в которых вектор у(М) направлен параллельно прямой Предположим, что кривая С обходится точкой М в положительном направлении, и пусть р есть число точек М^, при прохождении через которые вектор у(М) проходит через направление прямой двигаясь против часовой стрелки. Пусть, далее, д — число точек М^, в которых вектор у(М) проходит направление прямой двигаясь то часовой стрелке. Точки М^, в которых ве ктор у(М) достигает направления двигаясь, скажем, по часовой стрелке, а потом начинает двигаться в обратном направлении (или наоборот), мы не будем принимать во внимание. Тогда индекс кривой С равен (р - д)/2.

Индексом (или индексом Пуанкаре) изолированной особой точки векторного поля V, соответствующего динамической системе, называет-

С

СС

точек поля V.

Сд

Лемма 3. Сумма индексов особых точек системы (9), располо-

Сд 9

Аналогично [6] сформулируем два утверждения, вытекающие из свойств главных изоклин и леммы 3.

С

точки с индексом, по модулю большим единицы.

Принцип нечетности. Если Сд 9

- Сд

мере ещё две особые точки с индексом +1, причем особые точки с ин-

-у

так что +1, -1, +1,.... Таким образом, число стационарных точек в Сд с ненулевым индексом всегда нечетное, и число особых точек типа узла (фокуса) всегда на единицу больше числа седел.

Заметим, что для всех значений параметров система (9) имеет стационарное состояние (х8,у8) = (0,1). Рассмотрим матрицу Якоби первых производных правых частей системы (9) в точке (0,1):

(Рх(0,1,Ы Ру(0,1,Ы й = 0 -к1 - к— - к3Ъ -к й и<>Зх(о,1) ^(о,1)^ ^ -к3г& о ).

Она имеет собственные значения

_ кг + кг_! + Т + + *з£о)2 + 4М5£о

" -2-'

Собственный вектор = (-А1,к3£о) соответствует собственному значению А1 < 0, и ^ = (-^,к3£о) — собственному значению А2 > 0. Стационарное состояние (хз, уз) = (0,1) системы (9) является седлом, причем поскольку -АДк3£о) > 0, обе входящие сепаратрисы лежат вне области О. Кроме того, -АДк3£о) < -1) и в окрестности точки (0,1) одна выходящая сепаратриса лежит вне области О, в то время как вторая располагается в О.

Рассмотрим у € [0, 1). Выразим £о из (15) как функцию от уз:

& = £Ы = 1

, 1 - Куз) , ,

«1 —ТГ~\--к-г уУ*

Цуз)

(16)

к3(у8)Уз

Выражение (16) в неявном виде задает зависимость координаты уз стационарного состояния (хз, уз) системы (9) от параметра£о> координата

хз

Если ф 0 для всех уа (Е [0,1), то по теореме о неявной функ-

ции можно определить зависимость уз = ^(£о) так, что £о-= 0. В этом случае для всех значений параметра £о система (9) в области О\дО имеет единственное положение равновесия (хз, уз), где уз =

Пусть = 0 при некотором у3 = у3. Тогда если й ф 0, то

при уз = уз функция £(уз) имеет точку экстремума, и при изменении £о вблизи £о = £(уз) изменяется число корней уз уравнения £о-£(уз) = 0, и в системе (9) происходит бифуркация стационарных состояний. Таким образом, точка М = (уз, £(уз)) на кривой

Ф = {(уз,£о) :£о = £Ы, уз €[0,1)}

является точкой поворота.

На рис. 1 кривая Ф, обозначенная через АКВСР, имеет две точки поворота ^В. Здесь коэффициенты уравнений (9) задаются соотно-

к

Ф имеется 5-образный гистерезис, т. е. при £ < в системе (9) суще-

уз

0.6 0.5 0.4

0.2 0.1 0

0 0.2 0.4 0.6 0.8

Рис. 1. Зависимость y-координаты стационарных состояний системы (9) от параметра £о: CD и KB\{K} — устойчивые стационарные состояния, BC\({B} U(C}) — седлa, AK — неустойчивые фокусы и узлы. Значения параметров соответствуют (6) и ki = 4 • 103.

к 0; при £ > £с стационарное состояние тоже определяется однозначно, а при £B < £ < £с существуют три стационарных состояния.

Одним из параметров, которым можно управлять в ходе эксперимента, является PCO — парциальное давление CO(g) в газовой фазе над поверхностью катализатора. От этой величины линейно зависит коэффициент к\ = k\PCO (k® = const). С целью исследования зависимости ys от двух параметров к\ и £о, рассмотрим множество

S = {(£о, h,ys) : £о = Ziysih), кг >Q, ys € [0,1)},

где для каждого к\ > 0 функция £{ys,ki) = £(ys). Тем самым кривая

k

довапие свойств поверхности Е будем проводить с помощью методов теории катастроф [4]. Из рассуждений, проведенных выше, следует.

k k

k

y

Рис. 2. Поверхность = {(&1,у) ■. у — уз е [О, 1),£ = £(уз)}, где £(уз) определено в (16).

стационарных состояний над плоскостью параметров к и £о является многолистной.

В дальнейшем будем предполагать, что значения параметров системы (9) соответствуют (6). Тогда поверхность Е, показанная на рис. 2, имеет сборку Уитни. Проекция П поверхности Е на плоскость параметров (к, £о) не взаимно однозначна. Существует множество значений к и £о, являющееся объединением двух кривых, которым соответствуют два прообраза отображения П (рис. 3, кривые 1), являющиеся образами складок на поверхности Е. Точка пересечения кривых 1 соответствует сборке на поверхности Е. Тем самым если при возмущении параметров точка (к,£о) пересекает кривую 1, то в системе происходит бифуркация, связанная с изменением числа стационарных состояний.

Таким образом, если значения (к,£о) рассматриваются из области Г^; ограниченной кривыми 1, то система (9) имеет три стационарных состояния. Для значений (£0,^1) Г2звв из^ система (9) имеет одно стационарное состояние.

5. Бифуркация стационарного состояния седло-узел

Предположим, что значения параметров (6). определяющих зависимости кДу), ] = — 1 ,2,3, фиксированы. Определим множество точек (к{п, £8п) на плоскости Для которых система (9) имеет стацио-

нарное состояние (х8п,у8п), являющееся седло-узлом, т. е. выполнены равенства

Р(х8п, у8п, £о ) = ^(х8п,увп) ~ О,

РХ (х8п , у8п , ^8 ) х8п , у8п) Ру( х8п, у 8п , ^8 ) Qx(Х 8п , у 8п) - О,

и. кроме того, седловая величина отлична от нуля:

Р-х{х8п, у8п, + ^^х8п, у8п, ^О ) Ф 0.

Тогда при малом возмущении параметров в системе (9) происхо-

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

4500 3900 3300 2700 2100 1500 900 300

0

Рис. 3. Ветви бифуркационной кривой, выделяющие на плоскости параметров область существования трех стационарных состояний системы (9).

Учитывая, что хэп = хц(узп) (13), получим из системы (17) зависимости к{п = кГ(Уап) И £Г(Узп):

£3„ _ к^1(узп)/2(у3„) - 92(Узп)11(Узп) ^

л (у /у вп) VУ5п)УЭп Му п

7,бп _ к-1(у8П) + къ(узп)узп^п 4,1

1-1Ып) ПУЗП1 { '

где /(у) определена в (14),

Л(у) = 2^Ы(1 - у)[1 - /Ы] + Ыу)у,

/2(y) = (1- y)

2к2(у)[1-/(у)]-^М(1-у)[1-/(у)]2

ыу) + ^у)/(у)

йЫ = My)y/(y) + (1 - y)[i - fy)]/(y)

, / Ч , dh(y)

ыу) = (у)/(у) + - у)[1 - ду)шу)-

Таким образом, для 0 < ys n < 1 выражения (18) и (19) определяют множество значений £о = Соn и ki = kfn, при которых система (9) имеет негрубое стационарное состояние седло-узел.

Пусть коэффициенты системы (9) определяются соотношениями (2), (3), (5) и (6). Тогда множество значений £о и к, при которых система имеет стационарное состояние седло-узел, показано на рис. 3. В этом случае если точка (£о, ki) лежит внутри (или вне) области, ограниченной кривыми 1, то соответствующая система (9) имеет три стационарных состояния (или одно).

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

состояния порождает однопараметрическое семейство периодических решений. Различные максимальные семейства периодических решений будет подробно изучаться при моделировании сложных автоколебаний, поскольку наличие таких максимальных семейств является одним из достаточных условий для существования сложной динамики (см. [7]) в полной модели (1).

ЛИТЕРАТУРА

1. Slavinskaya Е. М., Stonkus О. A., Gulvaev R. V., Ivanova A. S., Zaikovskii V. L, Kuznetsov P. A., Boronin A. I. Structural and chemical states of palladium in Pd/AhOe catalysts under self-sustained oscillations in reaction of CO oxidation // Appl. Catalysis A: General. 2011. V. 401. P. 83-97.

2. Lasbina E. A., Slavinskaya E. M., Gbumakova N. A., Stonkus O. A., Gulvaev R. V., Stadnicbenko A. L, Cbumakov G. A., Boronin A. L, Demidenko G. V. Self-oscillations in CO oxidation on Pd0/Al203 catalyst // Chem. Eng. Sci. 2012. V. 83. P. 149158.

3. Андронов А. А., Леонтович E. А., Гордон If. И., Майер А. Г. Качественная теория динамических систем на плоскости. М.: Наука, 1966.

4. Арнольд В. И. Теория катастроф. М.: Наука, 1990.

5. Баутнн Н. Н., Леонтович Е. А. Методы и приемы качественного исследования динамических систем на плоскости. М.: Наука, 1976.

6. Чумаков Г: А. Динамика нелинейной системы дифференциальных уравнений // Сиб. мат. журн. 2007. Т. 48, № 5. С. 1180-1195.

7. Чумаков Г: А., Слннько М. Г. Кинетическая турбулентность (хаос) скорости реакции окисления водорода на металлических катализаторах // Докл. АН СССР. 1982. Т. 266, № 5. С. 1194-1197.

г. Новосибирск

15 октября 2013 г.