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

130 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

MSC 80A30

АНАЛИЗ СТОХАСТИЧЕСКОЙ МОДЕЛИ ХИМИЧЕСКОЙ КИНЕТИКИ БИНАРНОЙ АВТОКАТАЛИТИЧЕСКОЙ РЕАКЦИИ

Фам Минь Туан, Ю.П. Вирченко

Белгородский государственный университет, ул. Победы, 85, Белгород, 308015, Россия, e-mail: virchObsu.edu.ru

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

Ключевые слова: уравнения химической кинетики, стехиометрические коэффициенты, стохастическая модель, уравнение Фоккера-Планка, стохастический дифференциал, плотность распределения, фазовый переход, бифуркация.

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

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

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

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 131

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

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

[2], касающиеся одномерной динамической системы с двумя свободными параметрами, которая названа ими в [1] генетической моделью. По замыслу авторов, она связана с описанием кинетики некоторого класса бинарных автокаталитических химических реакций. В настоящей работе мы сохраняя идеологию работ [1,2], расширим результаты их исследований, добавив в эту систему еще один параметр, который допускает естественную физическую интерпретацию. Следующие второй и третий разделы работы посвящены, соответственно, определению исследуемой модели в детерминированном случае и выводу уравнения для стационарной плотности распределения фазовой координаты системы при мультипликативном возмущении параметра модели белым шумом. В заключительных 4-м и 5-м разделах мы подробно исследуем равновесную плотность распределения и, в частности, ее качественную перестройку при индуцированном шумом фазовом переходе. 2

2. Детерминированная модель. В этом разделе мы рассмотрим модель кинетики химической реакции в случае идеально термодинамически равновесной среды (без учета случайных флуктуаций).

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

132 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

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

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

A + X + Y k2 2Y + A*

ki

B + X + Y 2X + B*

кз

где X, Y, A, B, A*, B* - символы химических реагентов и при этом вещества, обозначаемые символами A, B, A*, B*, выполняют роль химической среды, в которой возможно протекание прямой и обратной реакции со сравнимыми друг с другом скоростями ki, i = 1, 2, 3, 4.

На основании базовых уравнений химической кинетики [3], описывающих динамику этих пар одновременно протекающих реакций, имеем

Nt(X) = k2Nt2(Y)Nt(A*) - кгNt(X)Nt(Y)Nt(A) + k3Nt(X)Nt(Y)Nt(B) - k4N2(X)Nt(B*),

(1)

Nt(Y) = kiNt(X)Nt(Y)Nt(A) - k2Nt2(Y)Nt(A*) + k4N2(X)Nt(B*) - k3Nt(X)Nt(Y)Nt(B) ,

(2)

где Nt(A), Nt(A*), Nt(B), Nt(B*), Nt(X),Nt(Y) - зависящие от времени t числа частиц соответствующих реагентов. Из этой системы уравнений следует закон сохранения суммарного числа молекул обоих реагентов в каждом физически малом объеме термодинамической системы, так как сумма двух уравнений приводит к d(Nt(X) + Nt(Y))/dt = 0. Тогда Nt(X) + Nt(Y) = N = const.

Обозначим посредством x(t) = Nt(X)/N, y(t) = Nt(Y)/N концентрации частиц, соответственно, реагентов X и Y в момент времени t. Согласно указанному закону сохранения, должно иметь место соотношение

x(t)+ y(t) = 1 ■ (3)

Так как кинетические уравнения для Nt(X) и Nt(Y) связаны, то в дальнейшем мы будем рассматривать только первое из них.

Типичная постановка химической задачи такова, что реагенты A, A*, B, B*, составляющие среду реакции, имеются в большом избытке. В этой ситуации их изменениями со временем можно пренебречь, так как они очень малы по сравнению с самыми величинами Nt(A), Nt(A*), Nt(B), Nt(B*). Поэтому далее мы опускаем при их записи аргумент t. По порядку величины числа N (A), N (A*), N (B), N (B*) совпадают, но намного превосходят числа Nt(X) и Nt(Y). В этой ситуации частоты химических превращений согласно указанным схемам реакций очень велики и по порядку величины пропорциональны произведению N(A)N. Поэтому мы перейдем к другому масштабу времени в уравнении (1), а именно, заменим N[k2N(A*) + k4N(B*)]t на t, разделив (1) на N2[k2N(A*) + k4N(B*)]. При такой замене временной параметр становится физически безразмерным. Тогда из (1) получается следующее уравнение для концентрации

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 133

x(t):

2 ( hN(B) - kiN(A)

х = ay — xy (

,k2N (A*) + k4N (B*) где введен безразмерный коэффициент

k2N (A*)

a =

k2N (A*) + k4N (B*)

j — (1 — a)x2

a e [0,1].

(4)

(5)

Заменим, в полученном уравнении x2(t) на x(t)(1 — y(t)), а y2(t) на y(t)(1 — x(t)), учитывая (3), и, введем параметр

Л

k3N(B) + k4N(B*) - hN(A) - k2N(A*) k2N(A*) +k4N(B*)

(6)

являющийся также характеристикой реакции, который может принимать любые по величине значения, так как зависимость параметра Л от управляемых извне переменных N (A), N (B), N (A*), N (B*) может быть определена в широких пределах, в принципе от — то до то. После таких преобразований, уравнение (4) принимает вид

x = a — x + Лx(1 — x).

(7)

Оно является основой для всех дальнейших построений. Оно имеет две стационарные точки, которые являются решениями квадратного уравнения a — x + Лx(1 — x) = 0. Одной из этих точек является

x

Л - 1 + у/(Л - l)2 + 4Ло:

2Л

(8)

Она всегда находится внутри интервала [0,1]. Так как значение производной правой части (7) по x в этой точке равно (-V(A- 1)2 + 4аА), где подкоренное выражение всегда положительно в силу a e [0,1], то это значение отрицательно, и поэтому стационарная точка x устойчива. Другая стационарная точка

Л- 1 - у/(Л - I)2 + 4Ло:

2Л

обязательно должна быть неустойчивой, на основании геометрических соображений, так как, в противном случае, между двух устойчивых стационарных точек автономного дифференциального уравнения первого порядка обязательно должна находиться еще одна неустойчивая стационарная точка. Она лежит вне интервала [0,1]. Следовательно, любое решение с начальным положением xo внутри интервала [0,1] обязано стремиться к устойчивой стационарной точке x.

В заключение конструкции модели заметим, что значения параметра a = 0, либо 1 являются особыми, так как при этом N(A*) = 0, либо N(B*) = 0. Неясно, имеет ли физический смысл такая ситуация.

134 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

Проинтегрируем уравнение (7)

dx

а — x + Ax(1 — x)

где xo G [0,1]. Выпишем явно решение при A > 0. Тогда, вычисляя интеграл, находим

у/(Х — I)2 + 4аА

ln

A i_УДДДунЫ 2/2 2

П 1 У(Л - I)2 + teA

2/2 2

= —t + ln C ,

где 0 < C < 1 - произвольная постоянная. Ограничение постоянной C сверху единицей связано с тем, что нам нужно решение дифференциального уравнения (7) с начальным условием x0 G [0,1] при t = 0 при A > 0. Наконец, явный вид решения

x

1 1 у/(А - I)2 + 4аА

2 ~~ 2А + А

1

1

1 _ Се~1^“ Х)2 + 4с^ 2

(9)

Это решение стремится при t ^ оо, как уже отмечено выше на основе рассуждений общего характера, к стационарной точке (8).

Так как точки а = 0,1 являются особыми с физической точки зрения, то это отражается и на явном решении (9) модели. При а = 1 имеется притягивающая точка {1, еслиА > — 1; — A-1, еслиА < —1}. При а = 0 единственной притягивающей точкой является {1 — А-1, еслиА > 1;0, еслиА < 1}.

Физическая интерпретация проведенного анализа математической модели (7) указывает на то, что с термодинамической точки зрения асимптотическая глобальная устойчивость стационарных состояний системы сохраняется при любой интенсивности связей, наложенных средой на систему. Иными словами, в идеально термодинамически равновесной среде не происходит качественных изменений динамики системы так, что у нее имеется единственное устойчивое стационарное состояние. Для рассматриваемой схемы химических реакций это означает, что при любых, даже очень больших отклонениях реального значения отношения продукт/субстрат (A*B*)/(AB) от его равновесного значения, определяемого законом действующих масс,

N(A*)N(B*) N(A)N(B)

равн

hh

fc2&4

(10)

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

X

t

1

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 135

3. Стохастическая модель. Модифицируем модель, построенную в предыдущем разделе, с точки зрения учета влияния термодинамических флуктуаций среды на протекание химических реакций. В модели, предложенной выше, предполагалось наличие идеального термодинамического равновесия среды и флуктуациями на фоне больших величин чисел N(A) и N(B), а N(A*) и N(B*), имеющихся в большом избытке, мы пренебрегали, следуя [1,2].

Так как эти флуктуации носят случайный характер, то их учет в рамках предполагаемой модификации должен приводить к тому, что такая модель должна быть, с необходимостью, стохастической. В полной мере учет флуктуаций среды может быть произведен только в рамках идеологии неравновесной статистической механики. Поэтому предлагаемая ниже стохастическая модификация феноменологической модели носит ограниченный характер. Простейший прием такого ограниченного учета флуктуаций среды, не отказываясь от термодинамических представлений, принятый в неравновесной термодинамике, состоит в том, чтобы ввести случайные возмущения в уже имеющиеся в описании системы термодинамические параметры (в нашем случае, параметры Л, а). Естественно, что такое возмущение нужно вводить посредством аддитивной добавки в термодинамические параметры некоторого стационарного случайного процесса. Этим мы предполагаем, что несмотря на наличие флуктуаций, система в целом находится в состоянии равновесия с точки зрения статистической механики. Кроме того условие стационарности стохастического возмущения, вообще говоря, является достаточным условием для того, чтобы система, испытывающая его воздействие, по истечении достаточно большого интервала времени, вышла на стационарный режим эволюции. В этой работе мы будем изучать стохастическую модель, которая получается из введенной детерминированной феноменологической модели возмущением только параметра Л. Этим самым мы предполагаем, что флуктуации чисел N(A),N(B) оказывают влияние на протекание химической реакции, в то время как флуктуации чисел N(A*),N(B*) пренебрежимо малы. (Числа N(A),N(B) не входят в определение параметра а.) Таким образом, мы далее будем считать, что

At = Л + (t, (И)

где внешний шум (t имеет нулевое среднее значение и его интенсивность равна а2. Здесь и далее, мы, для четкого смыслового выделения тех математических объектов, которые являются случайными, вводим для них специальное обозначение в виде постановки знака «тильда» над соответствующим объектом.

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

136 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

образом, что Zt = Тогда случайные колебания термодинамического параметра Л выражаются следующим образом

At = Л + &Zt ■ (12)

Заменяя в (7) параметр Л на At и умножив после этого уравнение на dt получим

dxt = [а — xt + ЛXt(1 — xt)] dt + axt(1 — xt)dwt, (13)

где искомая неизвестная функция xt переобозначена на xt, так как она, по своему смыслу, стала случайной функцией, и положено формально dwt = ,dt. Здесь wt - траектории случайного в обычном смысле (не обобщенного) процесса, который называется стандартным винеровским процессом. Однако, траектории винеровского процесса не имеют поточечно понимаемых производных в любой точке t, так как A - обобщенный случайный процесс. По этой причине, выражение dwt в (13) не является дифференциалом в обычном принятом в математическом анализе смысле. Это так называемый стохастический дифференциал. Этот дифференциал не имеет однозначного толкования. Его точное математическое определение и, соответственно, определение дифференциального уравнения (13), которое называется стохастическим дифференциальным уравнением, зависит от того, для описания каких физических ситуаций конструируемая математическая модель в виде такого стохастического дифференциального уравнения применяется.

В нашем случае, естественно, рассматривать уравнение (13) как идеализацию дифференциального уравнения

dAt = [а — xt + ЛXt(1 — At)] dt + xt(1 — At)dZ(t, (14)

где траектории случайного процесса ,t дифференцируемы в обычном смысле так, что d,t - обычный дифференциал функции ,t, и уравнение (14) представляет собой обычное дифференциальное уравнение, в котором, однако, имеются случайные параметры. Это означает, что представленная запись дает не одно дифференциальное уравнение, а, напротив, целое их многообразие, параметризованное множеством траекторий ,t и на этом множестве задано распределение вероятностей, которое индуцирует распределение вероятностей на множестве решений уравнения (14). В этой ситуации имеет место теорема Вонга-Закаи (см., например, [4]): Если в распределении вероятностей случайного процесса ,t перейти к пределу в слабом смысле, то есть перейти к пределу в характеристическом функционале процесса, к характеристическому функционалу обобщенного случайного процесса стА, то распределение вероятностей решений дифференциального уравнения (14) переходит в распределение вероятностей решений стохастического дифференциального уравнения (13), в котором стохастический дифференциал dwt должен пониматься по Стратоновичу [5].

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

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 137

по Стратоновичу, эквивалентно некоторому видоизменению уравнения (13), в котором стохастический дифференциал понимается по Ито [6]. Это измененное уравнение строится из (13) по фиксированному правилу и, в нашем случае, имеет вид

dxt

а - Xt + AXt(1 - Xt) +

а

2

Xt(1 - Xt)(1 - 2Xt)

dt + axt(1 - Xt)dwt.

(15)

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

а2

f(x) = а — х + Аж(1 — х) + —ж(1 — ж)(1 — 2ж) , (16)

g(x) = x(1 - x) , (17)

в терминах которых уравнение (15) принимает вид

dxt = f (Xt)dt + ug(Xt)dwt.

(18)

Тогда указанное выше эволюционное уравнение имеет вид уравнения Фоккера-Планка [1]:

dp(x,t) _ d[f(x)p(x,t)} a2 d2[g2(x)p(x,t)}

dt дх 2 дх2 1 J

Так как каждая отдельно взятая случайная траектория процессов, определяемых (14) (соответственно, (15)) не имеет физического смысла, то не имеет смысла потраектор-ный поиск решений этих уравнений. Физический смысл имеют только статистические характеристики указанных процессов и, в частности, плотность распределения p(x,t). Поэтому базовым уравнением для анализа сконструированной нами модели будет, в дальнейшем, служить именно уравнение (19).

4. Стационарная плотность распределения. В этом разделе мы вычислим стационарную плотность распределения p(x) - независящее от времени решение уравнения (19). При этом мы оставим в стороне вопрос о том, является ли эта стационарная плотность равновесным (то есть устойчивым по отношению к малым возмущениям) решением, и не будем решать вопрос о том, при каких условиях решение p(x, t) уравнения (19) стремится при t ^ оо к этой стационарной плотности.

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

dp(x,t) dJ(x,t)

dt дх

где плотность потока вероятности J (x) выражается формулой

J (x, t)

f (x)p(x,t)

a2 d[g2(x)p(x, t)]

2

dx

138 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

Тогда дифференциальное уравнение, которому удовлетворяет стационарная плотность p(x) имеет вид

dJ (x)

dx

0.

Общим решением этого уравнения является J(x) = J = const. В частности, поток вероятности на границе пространства состояний нашей термодинамической системы -отрезка [0,1] равен J.

Введем вспомогательную функцию

q(x) = g2(x)p(x).

С учетом постоянства потока вероятности, для нее получается замкнутое дифференциальное уравнение первого порядка

dq(x)

dx

2 f(x) а2 д2(х)

2

q{x) =------J .

а2

Общее решение этого уравнения

2 f f(x)

q(x) = A exp

a2

g2(x)

dx

2J

a2

exp

1 fJM.

a2 J g2(x)

dx

dx

зависит от одной произвольной постоянной A. В этой записи все интегралы понимаются как неопределенные. Таким образом, общее выражение для искомой стационарной плотности вероятности имеет вид

p(x)

A

g2(x)

exp

2

_ f(x) a2 J д2(х)

dx

2J

exp

2

_ f(x) a2 J д2(х)

dx

dx .

(20)

a2g2(x) j

Мы будем далее рассматривать это решение только для случая J = 0. Основанием к такому ограничению является то, что только при этом условии границы фазового пространства - точки 0 и 1 являются естественными границами для диффузионного процесса в смысле Гихмана-Скорохода [7] во всей области допустимых значений параметров Л и а, то есть диффундирующая частица не может преодолевать эти границы. Это свойство связано с расходимостью следующих интегралов

г 1

exp

2

а2

/Ы

д2(и)

du

dx

2

exp

: /Ы

_а2 J д2(и)

du

dx

для функций f (x) и g(x), определяемых (16) и (17). Здесь а и b - произвольные точки из [0,1], сколь угодно близкие, соответственно, к 0 и 1. Итак, общее выражение для стационарной плотности распределения имеет вид

p( x)

A

1 [Ш

д2(х) ^ [a2 J д2(х)

■ exp

dx

(21)

Подставим в общее выражение (21) явные выражения для функций f (x) и g(x), соответствующие нашей модели (13)

p(x)

A

x2(1 — x)

2 exp

—Я(ж)

а2

(22)

о

b

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 139

H (x)

а — x + Ax(1 — x) + a2x(1 — x)(1 — 2x)/2 ----------------------------------------*

Вычисление неопределенного интеграла H (x) дает

H(x)

а — 1 1x

а

—Ь (2а + А — 1) In

x

x

1x

2

+ у 1п[ж(1

х)].

и, после этого выражения в (22), формула для стационарной плотности принимает вид

p(x)

А / ж \ ^ ( 2 /а — 1

ж(1 — ж) \ 1 — х) Р \ <т2 \ 1 — ж

а

x

(23)

где

в

2(2а + А - 1) и2

(24)

Плотность распределения (23) теряет смысл при а = 0 и А < 1, так как в этом случае в < 0 и плотность p(x) неинтегрируема в окрестности точки х = 0. По той же причине она теряет смысл при а =1 и А > — 1. При этом в А 0 и p(x) неинтегрируема в окрестности точки x = 1. В дальнейшем, мы исключим из анализа значения 0 и 1 параметра а.

Постоянная A в этом выражении определяется, естественным образом, из условия нормировки плотности p(x),

1

!p(x)dx=ь ■

0

то есть

A-1

1 я

f dx f x у f2 /а — 1

J ж(1 — ж) \ 1 — ж/ eXP{^ \1 -x

а

x

(25)

где интеграл, очевидным образом, сходится.

Для вычисления нормировочного интеграла (25) совершим в нем замену переменной интегрирования

x 1 — а

1 x а

u

e

5

A

-1

2

exp {----- — /3 In

a2

1а

а

exp

----Wa(l — a) eh и + /3u\ du

a2

Ж

то есть

A

1

ГХР

—;+ P In

a2

1—а

а

K-e

4 r— a2

а)

i -i

1

(26)

140 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

где К-в(•) - модифицированная функция Бесселя второго рода с показателем (—в), которая для любого показателя v £ C и положительного x определяется интегральным представлением

СЮ

Kv(x) = ^ J e~xchu~uudu.

— Ю

5. Фазовый переход, индуцированный шумом. Нашей дальнейшей задачей состоит в установлении качественного устройства плотности p(x) и, в частности, определения при каких условиях на параметры модели а £ [0,1], Л £ R, а2 > 0 возможен фазовый переход - качественная перестройка этой плотности.

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

p'(x)

a2g2(x)

f(x) — a2g,(x)g(x) p(x) = 0 .

Уравнение a2g(x)g'(x) = f (x) приводит к следующему

2

P{x) = a — x + Аж(1 — x) ——ж(1 — ж)(1 — 2ж) = 0 , x £ [0,1]. (27)

Так как это уравнение кубическое, то оно имеет хотя бы один корень. Этот корень расположен на [0,1], т.к. P(1)P(0) < 0. Таким образом плотность имеет хотя бы один экстремум. Заметим, что P(x) и f (x) совпадают при а2 = 0. Уравнение f (x) = 0 при а2 = 0 определяло неподвижные точки для детерминированной системы. При этом известно, что у детерминированной системы имелась ровно одна стационарная точка. По видимому такое положение должно сохраняться при очень малых значениях а2, то есть при очень малой амплитуде шума плотность p(x) должна иметь ровно один экстремум. Этот экстремум должен быть максимумом, так как на концах отрезка [0,1] выполняется p(0) = p(1) = 0 (см. (23)). Ввиду того, что уравнение (27) кубическое и, следовательно, оно не может иметь более трех вещественных решений, то плотность p(x) может иметь либо одну, либо три экстремальных точки. В последнем случае, эти экстремальные точки представляют собой два максимума и один минимум, который находится между ними.

С целью упрощения исследования уравнения (27) в общем случае, перейдем в нем к другой переменной x = y + 1/2, при этом (27) преобразуется в уравнение

а 2Р{у+ 1/2) = у* + \у2 +

а2

1

а2

1

4

y+

"1 —а 1 Л 0 y£ ' 1 1"

|_2 <72 4<т2 ’ [ 2’ 2J

. (28)

Это уравнение для исследования перехода не пригодно при точном равенстве а2 = 0.

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 141

При исследовании уравнения (28) будем различать два случая а) 12а2 — 3а4 — 4А2 = 0 и b) 12а2 — 3а4 — 4А2 = 0.

Так как 12а2 — 3а4 — 4А2 = 0 или, что то же самое 3(а2 — 2)2 + 4А2 = 12, то условие a) выполнено только на поверхности цилиндра в координатном пространстве трех параметров (А, а2, а), описывающих исследуемую модель. В основании этого цилиндра лежит эллипс с центром (0, 2) в плоскости а = 0 и полуосями (\/3, 2). При этом условие

a) может иметь место только при А2 < 3. Только при таком условии имеются точки на поверхности цилиндра.

При выполнении условия a) уравнение (28) переходит в следующее

А 3 1

b+vj + а 2

1

а2

А(1 + 2а2)

9Й

0.

(29)

Из этого уравнения следует, что имеется только один вещественный корень

у

*

А / А(1 + 2а2)

За2 \ 9а4

1

--а

1

а2

1/3

(30)

Этот корень может определять единственный экстремум в точке x* = y* + 1/2 плотности p(x), если | у * | < 1/2. Покажем, что это неравенство всегда имеет место. Для его выполнения необходимо и достаточно, чтобы

А

З^2

А(1 + 2а2) _2

S 9Л +СТ

3

а~\

<

А 1 За2 2

3

Производя явное возведение в кубическую степень и пользуясь условием a) избавимся в крайних частях неравенств от кубической степени А,

А2 1

6а4 + 8

а

-2

а~\

<

А2 1

6а4 + 8

Избавляясь с помощью условия a) от А2 приходим к очевидному неравенству 0 < а < 1. Заметим, что в случае

1

2~а

А(1 + 2 а2) 9а2

0

(31)

уравнение (29) имеет утроенный действительный корень у* = —А/3а2. Таким образом, в этом случае три экстремальные точки сливаются p(y*) = pw(y*) = pw(y*) = 0. Условие |у*| < 1/2 приводит к ограничению а2 > 1, то есть только в таком случае указанный утроенный корень имеет место. Если же а2 < 1, то это означает, что такого утроенного экстремума нет, но при этом какой-то экстремум внутри отрезка [0,1] по x все же существует.

Рассмотрим теперь случай b). Мы не будем выписывать выражение для вещественных корней кубического уравнения (28) в общем случае, которое определяет расположение экстремальных точек. Вместо этого, найдем условие, при котором появляются кратные корни этого уравнения. Это условие, записанное в виде уравнения некоторой

142 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

поверхности У в пространстве параметров модели (Л, а2, а), которая разделяет области, с одной экстремальной точкой плотности распределения и, c двумя максимальными точками.

Наличие кратного корня в уравнении (28) приводит к необходимости существования общего решения у двух уравнений: (28) и условия обращения в нуль производной по у левой части уравнения (28),

а-2Р '(у +1/2)

3у2 +

2Л

~^У +

а2

1 1

~й2 4

0,

у е

1 1 2’ 2

(32)

Наш анализ основан на следующем простом утверждении, связанным с алгоритмом Евклида для полиномов.

Лемма 1. Пусть полином Q(z), z е C с единичным коэффициентом при старшей степени z имеет кратный корень z0, то есть он представляется в виде Q(z) = (z — z0)2R(z), где R(z) - полином степени deg Q(z) —2. Тогда полиномы Qi (z), l = 0,1,..., n—1, получаемые в результате применения алгоритма Евклида к паре полиномов Q(z) = Qo(z) и Qi(z) = Q'(z)/ai, Q'(z) = 2(z — z0)R(z) + (z — zo)2R'(z),

Qi(z) = (z — zi)Qi+i(z) + ai+2Qi+2(z), l = 0,...,n — 2; (33)

Qn- l(z) (z zn-l)Qn(z) , (34)

где коэффициенты ai, l = 1,...,n, выбираются таким образом, чтобы коэффициенты при старших степенях в Qi(z) равнялись единице, n < deg Q(z) — 1 и degQn = 1 так, что deg Qi+2(z) < deg Qi+i(z), l = 0, ...,n — 1 и, поэтому, degQi(z) < degQ(z) — l, обладают свойством Qi(z0) = 0, , l = 1,..., n.

□ Утверждение получается пошаговым переходом от значения (l + 1) к значению (l + 2). На нулевом шаге Q0(z0) = Ql(z0) = 0 в силу кратности корня z0. На всех последующих шагах, если Qi(z0) = 0 и Qi+l(z0) = 0, то из (33) следует, что Qi+2(z0) = 0. Процесс должен оборваться, на (n — 1)-м шаге, когда deg Qn = 1, и при этом остаток от деления Qn+l(z) на этом последнем шаге должен быть тождественно равным нулю. I

Следствие. Если полином Q(z; al, ...am) зависит от параметров al,..., am, изменяющихся в некоторой допустимой для них области Q С Cm и допускает существование кратного корня z0 для некоторых значений набора (al,..., am) из области Q, то требование тождественного равенства нулю не зависящего от z последнего остатка от деления Qn+l(al, ...,am) = 0 в алгоритме Евклида приводит к уравнению гиперповерхности У коразмерности 1 в области Q, для которой отношение (al,.., am) е У является необходимым и достаточным условием существования кратного корня у полинома Q(z; al,..., am).

Для применения сформулированного общего утверждения к исследованию наличия кратного корня у полинома а-2Р(у + 1/2) = Q(y) нужно положить степень полинома равной 3, то есть требуемое применение алгоритма Евклида обрывается на первом шаге,

Q(W; а,Л,а2)

у3+4++

а2

1

а2

1

4

у+

1 а 1 Л -0, уе ' 1 1"

|_2 J а2 4<т2 ’ [ 2’ 2J

1

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 143

„ , , 2\ 2 2Л 1

Qi{y’,OL,\,a)=y +7Г^У+ т,

так как на нулевом шаге имеем

Q(y; а,Х,а2) = Qi(y; а,Х,а2)

1 1

<т2 4

У +

Л

3^2

0, у е

1 1 2’ 2

+ a2Q2(y;а А, а2),

и на первом -

Qi(y; а, А,а2)

У + У* +

2А

Зет2

Q2(y; а, А, а2) + a3Q3(a, А, а2) ,

где

a2Q2(y; а, А, а2)

12а2 — 3а4

18а4

<+<Зз(<+ А, <т2) — -

4А2

----У +

1 У а2 4

"1 1 1 А(1 + 2а2)

2 а2 9а4

2 2А

+ У* + ~Х^2,У*

и введена величина

у*

2 А( 1 + 2<j2) + (18а: — 9)<т2 12а2 - За4 - 4А2

(35)

которая существует в силу условия b).

Условие существования двойного экстремума Q3(a, А, а2) = 0 стационарной плотности вероятности дает :

а 2Р'(у* + 1/2)

3y2 +

2А

“+У* + а2

1

<т2

1

4

0,

(36)

Это условие должно выполняться при подстановке явного выражения (35) y*. Оно может иметь место в том и только в том случае, когда дискриминант квадратного уравнения относительно y* положителен, то есть имеет место

4А2 + 3а4 - 12а2 > 0 . (37)

Это означает, поверхность У, разделяющая области, с одним и двумя максимумами лежит вне цилиндрической поверхности {(А, а2, а) : 4А2 + 3а4 — 12а2 = 0}.

Согласно утверждению леммы, y* является двойным корнем уравнения а-2Р(y + 1/2) = 0. Второй корень y0 полинома а-2Р(y +1/2), при котором реализуется максимум плотности p(x), p(y0 + 1/2) = max, выражается формулой

Уо = -2у* - ^ • (38)

а2

Для того, чтобы считать, что двойной корень y* является как раз тем корнем, при значении которого плотность p(x) испытывает качественную перестройку - возникновение второго максимума, нужно доказать, что |y*| < 1/2, то есть соответствующее значение х* = 1/2 + y* находится в области определения плотности p(x) и при этом

144 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

должно также иметь место |у0| < 1/2. Доказательство этого утверждения мы оформим в виде леммы.

Лемма 2. Если у*, определяемое (35), удовлетворяет уравнению (36) при выполнении условия (37), то при a = 0,1; а2 = 0 выполняется неравенство |у*| < 1/2.

□ Покажем, что в условиях леммы невозможен знак равенства |у*| = 1/2. Докажем, например, что у* = 1/2. Неравенство у* = -1/2 доказывается аналогичным образом.

Допустим противное у* = 1/2. На основании уравнения (35), имеем равенство

1 _ 2А(1 + 2а2) + (18сс — 9)<т2

2 “ 12а2 - За4 - 4А2 '

С другой стороны, из (36) следует а-2Р(1) = 1/2 + А/а2 + 1/а2 = 0, то есть —А = 1+а2/2, А < 0. Подставим это выражение в написанное выше равенство. Откуда следует (1 — а)а2 = 0 и, следовательно, наше предположение не выполняется.

При a Е (0,1), А > 1 и достаточно малых значениях а2 из (35) следует, что |у*| < 1/2 при |А| > 1. Так как у* при фиксированном значении а, как функция от А и а2, изменяется непрерывно в области вне эллипса {(А, а2) : 4А2+3а4 —12а2 = 0}, то для того, чтобы попасть в область |у* | > 1/2, оно должно при каких-то значениях параметров А и а2 обратиться в ±1/2, что невозможно. I

Замечание 1. Равенство у* = 1/2 возможно при а = 1, если 1 < а2 < 2 для того, чтобы второй корень у0 попал в отрезок [—1/2,1/2].

Замечание 2. В случае b) невозможно возникновение утроенного корня уравнения а-2Р(1/2 + у) = 0. В самом деле, для существования тройного корня нужно, чтобы двойной корень у* совпал с у0 = — 2у* — А/а2. Это дает у* = — А/3а2 и такое выражение должно быть корнем уравнения (36). Но так как решения этого уравнения имеют вид

у±

А ^4А2 + 3а4 — 12а2^

За2 6а2

то такое положение противоречит условию (37).

Таким образом, ввиду существования удвоенного корня у* в разрешенной для него области изменения, можно утверждать, что уравнение (36) определяет, на основании следствия Леммы 1, поверхность У в пространстве параметров (А, а2, а), при пересечении которой происходит фазовый переход. Уравнение этой поверхности, получаемое из (36), имеет вид

27а4 (а — 1/2^ 2 + А(9а4 + 18а2 — 4А2) (а — 1/2^ = А4 + А2 (1 — 5а2 — а4/^ — 4а2 (1 Отсюда следует, что поверхность может быть составлена из двух частей

а~\ = 54“! |а(4А2 - 18а2 - 9а4) ± ^ л/(А\2 + За4 - 12а2)31 ,

3

*74) .

(39)

склеенных по эллипсу в плоскости а = 1/2.

НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31 145

Множество общих точек поверхности (39) и поверхности {(Л, а2, а) : 4А2 + 3а4 5 6 7 — 12а2 = 0} эллиптического цилиндра приводит к уравнению (31). Однако, можно показать, на чем мы не останавливаемся в этой статье, что часть поверхности, определяемой уравнением (39), которая соответствует значениям параметра а £ [0,1] лежит в полупространстве {(Л, а2, а) : а2 > 4}.

Заключение. Наглядно, изученный нами фазовый переход состоит в следующем. Пусть частицы сортов X и Y окрашены в разные цвета, например, желтый и синий. Тогда реакционная смесь веществ в области изменения параметров модели c одним максимумом плотности распределения, имеет какой-то цветовой оттенок, соответствующий значению концентраций х* и 1 — х*. Вне этой области, в которой произошел фазовый переход, наблюдаются колебательные изменения цвета реакционной смеси смещаясь то в сторону концентраций х+, 1 — х+, то в сторону концентраций х-, 1 — х-, где х± - точки двух максимумов плотности р(х), пребывая в каждом состоянии некоторое случайное время, но, в среднем, одинаково долго. При фазовом переходе, в реакторе наблюдается два хорошо различимых макроскопических состояния.

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

Литература

1. Хорстхемке В., Лефевр Р. Индуцированные шумом переходы: Теория и применение в физике, химии и биологии / Пер. с англ. /М.: Мир, 1987. - 400 с.

2. Arnold L., Horsthemke W., Lefever R. White and coloured external noise and transition phenomena in nonlinear systems // Zs. Phys. - 1978. - B29. - P.367-373.

3. Яблонский Г.С., Быков В.И., Горбань А.Н. Кинетические модели каталитических реакций / Новосибирск: Наука (Сиб. отделение), 1983. - 256 с.

4. Wong E., Zakai M. On the convergence of ordinary integrals to stochastic integrals // Ann. Math. Stat. - 1965. - 36. - P.1560.

5. Stratonovich R. L. A new representation for stochastic integrals and equations // SIAM J. Control. - 1966. - 4. - P.362.

6. Ito K. Stochastic differential equations on a differentiable manifold // Nagoya Math. J.. -1950. - 1. - P.35.

7. Гихман И., Скороход А.В. Стохастические дифференциальные уравнения / Киев: Нау-кова думка, 1968.

146 НАУЧНЫЕ ВЕДОМОСТИ

Серия: Математика. Физика. 2013. №12(155). Вып. 31

ANALYSIS OF CHEMICAL KINETICS STOCHASTIC MODEL OF BINARY AUTOCATALYTIC REACTION

Pham Minh Tuan, Yu.P. Virchenko Belgorod State University,

Pobedy Str., 85, Belgorod, 308015, Russia, e-mail: virchObsu.edu.ru,

Abstract. The stochastic model of binary chemical reaction is studied on the basis of chemical kinetics equations with stochastically perturbed parameters. The perturbation is described by generalized random process named "white noise". The stationary first order marginal probability distribution density of the random process describing the evolution of concentration of fixed reaction component is investigated. In contrast to investigations of this model fulfilled earlier, our study of the phase transition induced by the noise is done in the three-parametric phase space of conditions mixture.

Key words: chemical kinetics equations, stoichiometric coefficients, stochastic model, Fokker-Plank’s equation, stochastic differential, distribution density, phase transition, bifurcation.