‘p-‘u'

., 3.‘ V
W' ”'25' ?l 1.73.; :‘5
. ' kit...

fi-fi'fl'
fim- ,

. rm
.nmu.
5"“ ‘V 2am

$32.55
'1 J -

w

c
‘v
‘11 ""

£34k,“ ,1 ’ LL, :9 “14' C __ . .. .
‘ v“! I atx'gfiiii ‘ ‘th. . r L
. 'v ‘I- 1w

. ‘3‘: A3533 ““331 x4; m‘c-“Z—é'ig‘fiféia: ”‘

“w“‘i‘éeat‘ .zfi'a’léimx‘ I7. {-1.25 .. y

L ' "1"" . :I‘ 7

3.1%“. ~54 .

_: “jg \; .‘r.“'"
u I”,

‘ .. '1'"li "‘4'";

'0' '.".":%:]f , z;

.. "I W: W;

IIUIIIqul' WI. ' ,

I . I I- nms .
l

" W-
I

" ': r ,I I'JIOVJ...‘|

H ‘r'; - -- u h
‘ .‘ ‘1." “My“. we! {_
‘I’;“ L‘”‘"’ 7 ‘ Iv' "I‘D Hit-F vim?

 

THE—251:5

This is to certify that the

thesis entitled
Hlt’pt' 15') be and that noel/LEW» bite-U cc- YL
Witcflo stud“ /\t¢ktifi3 at (ML 54* “ (53L
V ' .
MC (M k presented by

J ~ Kctt'a'lxkh}

has been accepted towards fulfillment
of the requirements for

 

rk . 1) degree in Ph");>t(dt CNnL‘LDt—CLX

l

{mi— 5- fir

Major professor

Date 3/7/3‘f

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

 

MSU

LIBRARIES
.-:—.

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

 

HYSTERESIS AND TRANSITIONS BETWEEN MULTIPLE STEADY STATES

OF THE SCHLOGL MODEL

BY

J Kottalam

A DISSERTATION

Submitted to
Michigan State University

in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Chemistry

1984

ABSTRACT

HYSTERESIS AND TRANSITIONS BETWEEN MULTIPLE STEADY STATES
OF THE SCHLOGL MODEL

BY

J Kottalam

Several interesting phenomena are associated with
multiple steady states in reacting systems far from
equilibrium. To study these phenomena in simplest form,
Schlogl introduced a cubic reaction model for autocatalytic
isomerization. In the parameter space for this model there
is a region in which two stable steady states and an unstable
steady state are found. In this dissertation spontaneous
transitions between the stable nonequilibrium steady states
of the Schlogl model and hysteresis in induced transitions
are analyzed.

The stochastic dynamics are described by a birth-and-
death master equation. The time-dependent solutions of the
master equation can be constructed using the eigenvalues
and the eigenvectors of the transition matrix in the master
equation. If the system parameters lie in the interior of
the multiple steady-state region, then the system undergoes
Kramers relaxation. The Kramers rule is valid for the Schlogl
model only if both the ratio of the two eigenvalues closest
to zero is substantially different from unity and the shapes
of the two peaks in the slowest decaying mode match the corres-

ponding peaks of the stationary mode. The region where the

Kramers rule breaks down is determined. The eigenvectors and
the time-dependent solutions are displayed for points close
to the boundary separating the single and multiple steady-
state regions and for points well in the interior.

If the system is prepared at a single steady state and
an externally controllable concentration parameter is increased
at a fixed rate and then decreased at the same rate, the
response of the system generates a hysteresis loop. This
effect has been simulated both deterministically and
stochastically. The deterministic results indicate a static
contribution to hysteresis in the multiple steady-state
region, while this contribution disappears when fluctuations
are taken into account. In the single steady-state region
both deterministic and stochastic results indicate that
hysteresis is purely dynamic.

Noyes has proposed a criterion for relative stability
of stable steady states in multistable systems. This criterion
is based on coupling two reacting systems by material exchange
and determining the final state of the coupled system. It is
shown in this dissertation that such a coupled system exhibits
more steady states than do the individual systems, and that
the mixing experiments are not suitable for a study of
relative stability. It is also shown that variation in the
rate of material exchange shifts the steady states of the
coupled system; these shifts are analyzed for cubic and

quartic models.

TO MY MOTHER
who showed me the pleasure of learning
AND FATHER

who showed me the value of education

ii

ACKNOWLEDGMENTS

I am grateful to all who assisted me during this
endeavour. The contributions of only a few are mentioned here.

Needless to say that this work and much of what I have
learnt would be impossible without the excellent guidance of
Dr. Katharine Hunt. She is a perfect example for the definition
of a guide: "Let not the wise disturb the minds of the
ignorant; let him, working with devotion, show them the joy
of good work." (Gita 3:26).

Sunil taught me computer programming. He was always there
to pick me up whenever the computer and other facets of life
trampled me down. The few friends I did have during my stay
at MSU served the purpose of having millions of friends.

Of all contributions from my wife, the only thing I find
words to express is that she typed the manuscript.

My interest in mathematics was created by my mother in my
childhood. My education is a manifestation of my parents'
intelligence and wisdom. Therefore, my achievements, if any,

are really theirs.

iii

LIST OF
CHAPTER

CHAPTER

2.3

CHAPTER

3.6

CHAPTER

TABLE OF CONTENTS

FIGURES

I INTRODUCTION

II MODERN METHODS FOR MACROSCOPIC CHEMICAL
KINETICS

The Schlogl model and the deterministic
analysis

The stochastic formulation

Methods of solution

III A REVIEW OF STUDIES ON THE SCHLOGL
MODEL

Further properties of the steady states

Coexistence of two phases

Approximate solution of master equation

Relative stability

Mean field theory and multivariate master
equation approaches to include local
fluctuations

Critical phenomena

IV A STUDY OF ASYMPTOTIC RELAXATION IN THE
SCHLOGL MODEL USING EIGENVECTORS OF
THE TRANSITION MATRIX

The Schlégl model and its stochastic

formulation

iv

vii

11

20

32
32
36
44

47

59

73

83

84

4.2 Transition between the two stable steady
states 87
4.3 The validity of Kramers relaxation 98
CHAPTER V HYSTERESIS IN TRANSITIONS BETWEEN MULTIPLE

NONEQUILIBRIUM STEADY STATES OF THE

SCHLOGL MODEL 127
5.1 Schlbgl model and its steady states 127
5.2 Deterministic simulation of hysteresis 129
5.3 Stochastic simulation 134

CHAPTER VI DETERMINISTIC STUDY OF MULTIPLE STEADY
STATES IN COUPLED FLOW TANK REACTORS 141

6.1 Multiple steady states in a single

continuously stirred flow-tank reactor 142

6.2 Coupled flow-tank reactors 146

6.3 Comparision with the Noyes analysis 168
CHAPTER VII FUTURE WORK AND DEVELOPMENT 173

APPENDIX A Proof that the transitions generating
hysteresis cannot occur before a marginal
stability point is reached according to the
deterministic rate equation 175
APPENDIX B Method and program to find the eigen-
values of the transition matrix 178

APPENDIX C Programs to construct the eigenvectors

of the transition matrix 183
APPENDIX D Deterministic simulation of hysteresis 204
APPENDIX E Stochastic simulation of hysteresis 208

V

APPENDIX F Special features of the cubic mechanism
in a flow tank reactor 213

REFERENCES 2 1 7

vi

LIST OF FIGURES

Figure 2.1 The steady states of the Schlbgl model
Figure 4.1 Eigenvectors of the transition matrix A
for the Schldgl model, when the parameters
correspond to a point in the interior of the

multiple steady-state region

Figure 4.2 Time evolution of probability distri-
bution when the parameters correspond to a point
in the interior of the multiple steady-state
region. Initial distribution is 6(x-500).

Figure 4.3 Time evolution of probability
distribution. Initial distribution is 6(x-50).

Figure 4.4 Time evolution of probability
distribution. Initial distribution is 5(x-253),
peaked at the unstable steady-state X-value

Figure 4.5 Time evolution of probability
distribution. Initial distribution is 6(x-230)

Figure 4.6 The first three nonzero eigenvalues of
the transition matrix A. c4 = 3.33333.

Figure 4.7 The first three nonzero eigenvalues of
the transition matrix A. c4 = 1.7

Figure 4.8 The first three nonzero eigenvalues of
the transition matrix A. c = 1.3 (no multiple

4
steady states)

vii

88

94

95

96

97

104

105

106

Figure 4.9 Regions of the parameter space where the
conditions on the eigenvalues (equation (4.31))
is satisfied.

Figure 4.10 Eigenvectors of the transition matrix Q
near a marginal stability point

Figure 4.11 Time evolution of the probability
distribution near a marginal stability point

Figure 4.12 Eigenvectors of the transition matrix
near the critical point

Figure 4.13 Time evolution of probability
distribution near the critical point. Initial
distributions are a) 6(x—120), b) 6(x-300), and
c) 6(x-182)

Figure 4.14 Time evolution of probability
distribution near the critical point and near
a marginal stability point. Initial distribu-
tions are a) 6(x-300), b) 6(x-100), and
c) 6(x-190).

Figure 5.1 Results from deterministic (curves with

arrows) and stochastic (noisy curves) simulations

of hysteresis in the multistable region (c4 = 1.7L

The S-shaped curve is the set of steady states
Figure 5.2 Results from deterministic (curves with

arrows) and stochastic (noisy curves) simulation

of hysteresis in the Single steady-state region.

(c4 = 1.3).

viii

108

110

115

116

120

123

131

132

Figure 5.3 Area inside the hysteresis loop versus
the rate of change of B. (Deterministic results) 133
Figure 5.4 Comparision of deterministic (upper line)
A and stochastic simulation (lower line) results.
(24 = 1.7. 139
Figure 5.5 Comparision of deterministic (upper line)
and stochastic simulation (lower line) results.
8 = 50. 140
Figure 6.1 Steady states for the cubic model in a

single flow tank reactor. 147

Figure 6.2 Steady states of coupled reactors

system for Noyes cubic model. k0 = 1.5><10_ss-1 150
Figure 6.3 Steady states of coupled reactor

system for the cubic model. k0 = 2.7><1o"5s‘1 152
Figure 6.4 Steady states of coupled reactor

system for the cubic model. k0 = 2.516963x10'5s‘1 153
Figure 6.5 The potential function m for the cubic

model with k0 = 2.517><lo"5s,‘1 and kx = o. 156
Figure 6.6 The potential function ¢ for the cubic

model with k0 = 2.517x10"5s'1 and kx = 4.0x10'5

s"1 157

Figure 6.7 The potential function m for the cubic

model . k0 = 2.517x1o‘55'1 and kx = 9.061x1o’5

3-1. 158
Figure 6.8 The potential m for the cubic model.

k0 = 1.5xlo'5s'1 and kx = 3.9x1o'ss'1 159

ix

Figure 6.9 The potential m for the quartic model.

k0 = 3.374x1o'5 s'1 and kx = 1.28x1o’5 3'1. 160
Figure 6.10 The potential m for the cubic model.

k0 = 3.5x1o"5 s"1 and kx = 1.1><1o'4 s'l. 161
Figure 6.11 Comparision of steady state patterns of

the single and coupled flow reactors. 172

CHAPTER I

INTRODUCTION

All efforts to understand natural phenomena involve
observation, deduction, induction, and investigation, which
are characteristic of human intellect. One of the methods is
mathematical modeling, in which a physical system is first
observed, mathematical relations are deduced, these relations
are applied inductively (or generalized) to other similar
systems, and these systems are investigated to verify the
relations. Since models are deduced from limited observations,
they need not be complete, and we should incorporate
modifications as the scope of the models widens to include
new fields of investigation.

A Chemist's ambition is to completely understand all
aspects of chemical reactions. One such aspect is the
variation in the concentrations of different species during
a reaction. These variations can be investigated at several
levels of precision. One extreme is to view matter as
continuous bulk material with the amounts of each substance
varying smoothly in time, governed by the law of mass action.
In this case the equations of motion for chemical concentra-
tions are ordinary differential equations. On the other hand,
the exact (classical) treatment would involve the study of
encounters between various molecules using the equations of

motion for the position and velocity vectors of each atom in

2

each molecule. Although such "molecular dynamics" calculations
are being carried out, a typical laboratory situation involves
0(1023) molecules each having 0(101) atoms or more. Moreover,
the number of observations actually made on a system in the
laboratory is very small compared to the number of independent
data elements that would be available from an exact theoreti-
cal treatment. Hence in statistical mechanics the encounters
between molecules are viewed as frequently occurring random
events and the attempt is made to derive the expected values
of the few observable quantities and the expected deviations
from those values. Thus we are often satisfied with the
macroscopic description. When the deviations are large the
system may exhibit phenomena which cannot be treated by
deterministic methods. Then we need to incorporate stochastic
character into the theory.

Let us consider a spatially homogeneous mixture of n

reacting chemical species whose particle numbers X(1), X(2),

, x‘n)

are components of X. In macroscopic chemical
kinetics the time evolution of X is given by a first-order

ordinary differential equation:

—— = Egg). (1.1)

This equation is nonlinear except in special cases and is

usually autonomous (i.e., {(X) is usually independent of t).

The steady states of the system are those points XS in
state-space at which the numbers of all chemical species are

stationary with respect to time.

.__ = 5(E.) = o, (1.2)

If the system is closed, we call the steady state an
equilibrium state. For systems operating at or near equilib-
rium conditions, the steady state is unique, and it is globally
and asymptotically stable. Asymptotic and global stability
mean that the system will eventually arrive at the steady state
starting from an arbitrary initial condition. This is a
consequence of the minimum entropy production principlel.

However, an open system forced away from the linear
regime near equilibrium by a matter flux across its boundaries
can exhibit exotic phenomena such as the presence of several
steady states and periodic oscillation. If one of these
structures (steady-states or oscillatory solutions of equation
(1.1)) is unstable, then the solutions of the deterministic
differential equations do not describe the system adequately.
If a one-variable system has more than one steady states, not
all can be stable. Suppose there are two stable steady states
X1 and X3 for the Single-variable system; then in the neighbor-
hood of X1 the system moves towards X1 and therefore away

from x3 (and similarly for X3). This implies the

4

existence of another point between X1 and X3 from whose
neighborhood the system can only move away; this is an unstable
steady state. By the same argument,2 two stable steady
states:h1a two-variable system must be accompanied by an
unstable steady state, a saddle point,or an unstable limit
cycle.+ In these situations we need to study the system in

more detail.

In the stochastic description, sometimes called the
mesoscopic description, the system can occupy any possible
state with a finite probability at a given time. Thus we use
a probability distribution over the species-number space in
place of the species number to specify the state of the system
at a given time. We again use the law of mass action, but we
use it to calculate the probability of transition from one
state to another and not to calculate the rate of the reaction
directly. Using the transition probability we arrive at a
partial differential equation - called the master equation —
governing the evolution of the probability distribution in
time.

For example, when the parameters take certain values,

the Schlégl model3 has three steady states.one of which is

 

I Near a saddle point the motion is towards the point in one
direction and away from it in a perpendicular direction. A
limit cycle is a periodic solution of a set of differential

equations in two dependent variables.

unstable.4 According to the deterministic rate equations only
one of the stable steady states can be reached from a parti-
cular initial point, but according to the master equations
both stable steady states and their neighborhoods can be
reached with significant probability independent of the
initial distribution. Moreover, when an external parameter is
varied beyond a threshold value, the system undergoes a transi-
tion from one stable steady state to another. When we reverse
this variation, the threshold for the reverse transition is
different. The extent of this hysteresis effect is markedly
different in the deterministic and stochastic descriptions.

In the next chapter the stability analysis of determinis-
tic steady states and the stochastic formulation are illustra-
ted using the Schlfigl model as an example. Methods for solving
the stochastic equations are also explained. Chapter III
provides a review of the literature on topics related to
multiple steady-state systems; coexistence of two phases,
relative stability of steady states, critical phenomena, and
alternate formulations of stochastic kinetics are discussed.
In chapter IV results for the time evolution of the probabi-
lity distribution as described by the master equation are
presented. In chapter V deterministic and stochastic simula-
tions of hysteresis in the SchlBgl model are presented. These
studies are based on a master equation which takes into account
uniform fluctuations in the particle number, but excludes

spatial fluctuations. For a more realistic description,

particularly in the vicinity of the critical point, both
kinds of fluctuations should be considered. In chapter VI
the effects of coupling two Open reacting systems are

studied. Chapter VII contains suggestions for future work

on these problems.

CHAPTER II

MODERN METHODS FOR MACROSCOPIC CHEMICAL KINETICS

2.1 THE SCHLOGL MODEL AND THE DETERMINISTIC ANALYSIS

A system is said to be within the linear regime around
an equilibrium state if the thermodynamic fluxes (e.g., rate
of reaction) are linear combinations of the forces (e.g.,
chemical affinity). In this region, it is well known that
there can be only one stable steady state.1 Several models
have been proposed to study chemical kinetics far away from
equilibrium. One model exhibiting multiple steady states is

Schlagl's termolecular model.3 It consists of the reactions

2 X + A 3 X

B X. (2.1)

Here the numbers of molecules of A and B are kept at
predetermined constant values by contact with external
reservoirs or by appropriate feeding into or removal from

the reactor. We assume the reactor to be efficiently

stirred so that the system is always homogeneous and diffusion
effects need not be considered. The constants k1, k2, k3, and

k4 are specific rate constants independent of the size of the

system. While the rate constants are characteristic of a
reaction mechanism, the concentrations of A and B can be
controlled externally. For convenience we will also set the
numbers of A and B molecules equal and denote the number by

B. Then the deterministic rate equation is
QX = k bxz - k2x3 + k b - k4x. (2-2)

where X = X/V and b = B/V are the concentrations of X and B
(and A) respectively, and V is the volume of the system.
The steady states of the system are the real solutions

of the algebraic equation
BX- X+kB-kX=0. (2.3)
$7 {:7

The steady states are displayed in Figure 2.1 for various
values of the parameters. It can be shown5 that for

k1k4

k2R3

< 9, there is only one steady state for all values

 

n:

of B, while for n > 9 there is a range of B values in which
there are three steady states (see Figure 2.1).

1, X2, and X3 be the roots of this cubic equation.
The stability of each steady state can be deduced by a linear

Let X

stability analysis. The rate equation can be written as

ISOO

ISSD-

1200L

I050b

750*
600~
450-
I 5 K) l5 0 7F30
300-

50

 

 

'0.0 02 0.4 0.6 03 1.0 1.2 L4 LG La 2.0
mefl

Figure 2.1 The steady states of the Schldgl model.
The number of X molecules at steady state as a function of
the number of B molecules for various values of the

rate constants.

10

k
dX _ _._2 _ _ _
5E - v2(X X1)(X X2)(X X3). (2.4)

Let us consider an arbitrarily small deviation 0Xi from Xi

(i=l,2,3). Substituting éxi = X - Xi in equation (2.4) gives

ddxi k2
_;:_ = - 57 (5xi+xi-xl)(6xi+xi-x2)(6xi+xi-x3)
= KiOXi (2.5)

to first order in 6X1, where

k2
K1 ’ ' 52(X1'X2’(X1'X3)

k2
K2 = ' ;7(x2-x1)(x2-x3)
k2
and K3 = - ;§(X3-Xl)(X3-X2). (2.6)

These linear differential equations have the solutions

axi = axi e 1 . (2.7)

The steady state Xi is stable, marginally stable,or

11

unstable according to whether an arbitrarily small deviation
dxi decays, remains stationary,or grows, and this in turn
depends upon the sign of Re(Ki). From the definition of Ki
one concludes the following: If all roots are real, with

X < X2 < X3, then X and X are stable and X is unstable.

1 3 2

1 or X3 (or both) it is

marginally stable. If there is a complex pair of roots, then

1

Whenever X2 coincides with either X
the real root is always stable.

2.2 THE STOCHASTIC FORMULATION

In order to study transitions between two stable steady
states where the average value of X reaches and passes through
the unstable steady state value, the deterministic analysis
is not adequate. For a better description, we consider the
number of molecules of X in a fixed volume at a given time
as a discrete random variable. Accordingly we focus on the.
probability P(x,t) that the random variable X has value x at
time t. The evolution of this distribution is determined by
the occurrence of chemical reactions at random as explained
below.

A stochastic process is completely determined by
specifying an initial distribution and a hierarchy of
conditional distributions; i.e.,

P(xn+1,tn+1lxo,t0;xl,tl;...;xn,tn) denotes the conditional

probability that the system has Xn+1 molecules of X at time

tn+1 given that it had x at time t

O 0' x at time t and so on,

1 1’

12

where x0, x1, ..., xn are d-dimensional vectors for a process

with d variables. If it happens that

P(xn+1'tn+1IXO’tO;X1’t1;'"7Xn'tn) = P(xn+1'tn+llxn'tn)
(2.8)

for all n, then the process described by this conditional
probability is called a Markov process. In other words, of

all conditions specified in the past, only the most recent is
relevant in determining the future evolution of a Markov
process. This does not imply that the future is independent of
the past, but this dependence is entirely contained in the
present. If the Markov process satisfies the additional

condition

P( It) = P(X

I
Xn+1'tn+1lxn n n+l'tn+1+Tlxn’tn+T)' (2'9)

then the process is said to be a stationary (or time-homo-
geneous) Markov process. The Markovian character of stochastic
processes corresponds to the first order nature of differential
equations whereas the stationary property corresponds to the
autonomous behavior of deterministic differential equations.
Just as an n-th order deterministic differential equation can
be transformed to a set of n coupled first order equations?

a non-Markov stochastic process requiring the specification

of n conditions for its complete description can be

transformed to a Markov process is nd variables?

.13

A Markov process is completely determined by specifying
an initial distribution and a transition probability
P(x,t(y,s). Using elementary theorems of probability theory

it can be shown that

P(x,t+s) = X P(x,t+s|y,t) P(y,t). (2.10)
y

This equation is called Kolmogorov's equation. Using first

principles of differentiation, we obtain the master equation8

3 _
§EP(x,t) — E Axy P(y,t). (2.11)

where

lim P(x,t-l-s Iy,t) - 0

s—+0 3

xy (2.12)

is the probability of transition per unit time from state y
to state x. AXy is called the transition probability rate or
infinitesimal transition probability when x # y. Using the

normalization condition on P(y,t+s|x,t) namely,

P(x,t+s|x,t) = 1 - 1 P(y,t+s|x,t), (2.13)
va‘x

we get

14

A = - 2 A . (2.14)

The postulates of stochastic chemical kinetics are
1) Chemical reactions are Markov random processes.
2) Given that the system has x molecules at time t, the
probability that a reaction step will occur in the interval

(t,t+6t) is
.1—
am(x,t)6t + o(6t), (2.15)

where am(x,t) is proportional to the number of reactant
combinations leading to the m-th step at time t. For example,
if the second step in a list of all elementary reactions is
3X ——> 2X + A, then a2(x,t) = c2x(x—1)(x-2)/3!

The second postulate gives the transition probability as

Axy(t) = ; am(y,t), (2.16)
where m ranges over all reaction steps which result in a net
change of y->x in the species number. The Markov property

is implied in the deterministic formulation also. The second
postulate is the analog of the law of mass action. Thus the

random variation of X is the only addition to the theory.

 

f The meaning of 0(6t): If 6a = 0(6t), then by definition

'1im 5a _ . _ 6‘
0t-901—E - 0, or equ1valently, o(5t) — 0(5t ), e > 1,

15

For example, the (infinitesimal) transition probability

 

for the Schldgl model is given bys’9
( a+(y) for y = x-l
a_(y) for y = x+1
Ax=1
y - a+(y) - a_(y) for y = x
L 0 ' for Ix-y|> 1, (2.17)
where
a+(x) = cle(x-l)/2! + c3B, (2.18)
a_(x) = c2x(x-l)(x—2)/3! + c4x, (2.19)

x and y range over nonnegative integers,and the ci are the
proportionality constants which are related to the rate
constants ki. The relation is made by taking the deterministic

limit in the master equation.10 We find

k. = v1 c./n.1 (2.20)

where n1 is the molecularity of the i-th step in X and mi is
the total molecularity of the i-th step (e.g., n1 = 2, m1 = 3).

The master equation for the Schldgl model is

§%P(x,t) = a+(X-1)P(X‘1rt) + a_(x+1)P(X+l’t)

- [a+(x) + a_(x)]P(x,t). (2.21)

16

The stationary distribution Ps(x) is the solution of
iLP(x t) = 0
at I I

i.e., a+(x-1)PS(x-l) + a_(x+1)Ps(x+l)

= [a+(x) + a_(x)]PS(x). (2.22)
From this the detailed balance relation
a+(x-1)Ps(x-1) = a_(x)PS(x) (2.23)

follows by induction. The detailed balance relation gives the

expression for the normalized Ps(x) as

X
W a+(y-1)/a_(y)

 

._ y=1
’Ps(x) - z . (2.24)
E a+(y-1)/a_(y)
zy=1
This can be written as5
Ps(x) = ps(0)e‘¢(x’ (2.25)

in standard form, where the "stochastic potential" ¢(x) is

given by

x a (y)
_ - . (2.26)
¢(X) _ ygllna+(Y'1)

17

Although the master equation should be regarded as
the fundamental equation in this treatment, it is often
convenient to work with other approximate equations such as
the Fokker-Planck equation and the Langevin equation which
may be solved by standard mathematical techniques. The master

equation can also be written as

§EP(X, ,t) = g AxyP(y,t) = E A(x—z,z)P(x-z,t),
(2.27)
where A(x-z,z) is the infinitesimal transition probability

for a change z ¢ 0 in the number of molecules in a system

containing x-z molecules, and A(x,0) is fixed by
2 A(x-z,z) = 0. Expanding each term on the right hand side
z

as a Taylor series about 2 = 0 in the first argument of

A(x-z,z) and in P(x-z) yields

5%P(x,t) = XXL-L1,“): z)P(x, t)](- 2:)m (2.28)
z m me
This is the Kramers-Moyal expansion.11 The Fokker-Planck

equation12 is obtained by neglecting terms corresponding to

m > 2. Thus it reads

 

Z

+
NIH
Ntv:

2[A(X,2)P(X.t)]. (2.29)

18

For the Schlogl model, from equation (2.17)

 

( =
a+(x) for z 1
- a+(x) - a_(x) for z = 0
A(x,z) = (
a_(x) for z = -1
(0 for |z| > 1. (2.30)

Hence the Fokker-Planck equation becomes

2
in t) - -—3—[ ( )p( t)) +l-a—ta(x)1>(x tn
at x' ‘ ax r X X' 2 ax2 ' ,
(2.31)
where
r(x) = a+(x) - a_(x), (2.32)
and a(x) = a+(x) + a_(x). ' (2.33)

In general neglecting terms with m > 2 cannot be justified.
Van Kampen13 and Kubo £5 3114 have developed systematic
expansion procedures in power series of a small parameter
(usually the inverse system Size). Gillespie15 has shown

that if A(x,z) = 0 for z > 1, then the Fokker-Planck

+
equation when discretized with unit step Size"

 

1-Discretizing with unit step size means approximating BP/ax

by AP/Ax with Ax = 1; then AP = P(x,t) - P(x-1,t).

1’9
reproduces the master equation exactly, whereas the conditions
A(x,12) # 0 and A(x,z) = 0 for |z| > 2 require four terms to
be kept in the Kramers-Moyal expansion. Matsuo et a116 have
discussed methods of constructing Fokker-Planck equations
from the deterministic equations describing physical processes.
Obtaining the coefficients r(x) and a(x) in the Fokker-Planck
equation from the Kramers-Moyal expansion is one of the methods.
Another method is to adjust the coefficients in such a way that
the resulting stationary distribution coincides with the
stationary distribution given by the master equation (2.24).
The two term Fokker-Planck equation is useful because it
describes the process specified by the stochastic differential

equationl.7

(I) dXt = r(Xt)dt + /a(Xt)th (2.34)

or, equivalently,

(S) dXt = [r(Xt) - % ax(Xt)]dt + Va(Xt)tha (2.35)

where (I) and (8) indicate ItO and Stratonovich calculi

respectively, W is the normal Wiener process and aX is the

t
derivative of a. Such processes are relatively well understood

mathematically. When a(x) is a constant independent of x,

the resulting stochastic equation

dX = r(Xt)dt + /Edw

t t, (2.36)

20

is the Langevin equation,which has played an important role

in the theory of Brownian motion.18

2.3 METHODS OF SOLUTION

Since the master equation is a linear differential
equation in P(x,th its solution can be written formally in
terms of the eigenvalues and eigenvectors of the transiton
matrix A, when A is independent of time.

Relation (2.14), namely

2 A = o for all x, (2.37)
y' yx

means that the rows of g are linearly dependent and this

implies det g = 0. Thus the equation

EAX

Ps(y) = 0 (2.38)
Y

Y

has a nontrivial solution under the constraint
2 P (x) = 1.
x s

This guarantees the existence of a normalized stationary
distribution. This distribution is unique if and only if all
states of the process communicate with each other; i.e., if
and only if any arbitrary state can be reached from any other

state in finite time. If this condition is not satisfied the

21

process degenerates into several disjoint processes.19 For the
SchlBgl model the requirement for uniqueness is satisfied,
since a+(x) > O for all x and a_(x) > 0 for all x > 0.

Since the stationary distribution is an eigenvector
corresponding to zero eigenvalue (see equation (2.38)), zero
is always an eigenvalue of A. Though detailed balance implies
the existence of a stationary distribution, the converse is
not true. However, detailed balance implies additional
properties of the eigenvalues.

In addition to (2.37) and AXy > 0 for all y # x, if
detailed balance also holds

( i.e., if Ax Ps(y) = AnyS(x) ), (2.39)

Y

then the nonzero eigenvalues of A are real and negative.8

Proof:

Consider the matrices
(2.40)
g = g I_\._____W. (2.41)

S is obviously real and it has the same eigenvalues
as A. Further
%

Ps(x) Axy Ps(y)

%

(I)
ll

XY

%

- -%
PS(X) AYXPS(X)PS(y)

22

5 8
A XPS(x)

= Ps(y) y

YX ; (2.42)
i.e., g is also symmetric. Hence the eigenvalues are

real. Let 9-j be an eigenvector of g corresponding to

an eigenvalue lj # 0.

Hm
M)
u

2 Q .S Q .
x,y X3 XY Y]

Q 2 + E Q

. Q
x3 x,y#x

2 S .s .
x XX X] xy y:

- - i A Q .2 + 2 Q

S .Q -
x,y#x YX X3 x,y#x XY X3 Y]

2 2
_ - 1 A 0 . — Z A Q -
le>x YX X3 X,Y<X YX X]

+ E Snyijyj + Z Snyij

x.y>x x.y<x Y]

 

 

 

 

 

Z Ps(y) L5 2
= - A Q - Q -
XIY>X xy yj PS(X) x3
5 0 since AXy > 0 for y ¢ x. (2.43)
By choice of Q-j'
Q S .
A = ——03 = Q.)
J 2 7
Q .
—.3
so (2.43) implies
A. g 0, (2.44)

3

which was to be proved.

Since A is not symmetric the eigenvectors of g and AT
are not the same, but a relation between them can be deduced
by use of g. Let the columns of Q be the normalized
eigenvectors of g in a particular order and A be the diagonal

matrix containing the eigenvalues. Then

(2.45)

k)
Hm
M)
u
H>
m
B
a
K)
)o
u
"H

P

m
M)
"2
u»
"2
K)

II
H>

or g A (2.46)

"w
II
>
s

24

 

where
g = E‘lg (2.47)
and .3 = .W Q (2.48)
Thus L j and R . are the j-th eigenvectors of §.and QT
respectively. Moreover,
f3 = .0de Q = 1
Eliminating Q between (2.47) and (2.48), we have
-2
g, = g: g . (2.49)
The solution of 8—851: = A}: can now be written (when A__ is
time-independent) as
ljt
P(x,t) = 20)]. ije , (2.50)
3
where
R..P(i,0)
= 13 . 2. 1
Otj l Ps(l) ( 5 )

and P(i,0) is the initial distribution. The eigenvalue

is

equation for the k-th eigenvector R-k

A.. . = ,, , ' =
z lJRJk AkRik l 0,1,...

25

Summing over i, it becomes

ZZA..R. =AER..
l J 13 jk k i lk
Since 2 Aij = 0, it follows+ that
i.
g Rik = 0 If 1k ¢ 0. (2.52)

This means that )P(i,t) = 1 provided {Ps(i) = 1.
i i

Thus we need to diagonalize the matrix A in order to
solve this problem completely. g is infinite dimensional, but
for a nonexplosive chemical reaction we can always find a large
region of the state space outside which there is no significant
probability to find the system. Since reaction steps of high
molecularity are not common, A is always a narrow band matrix.
For the Schlbgl model it is tridiagonal.

The matrix diagonalization method is valid only if the
transition probabilities are time-independent,i.e., if the
system is autonomous. Gillespie has devised a numerical
algorithm for stochastically Simulating any chemically

reacting system10 and he has extended it to include

 

.1-

assuming that ZAin. converges uniformly in i so that we

jk
can interchange the order of the two infinite summations.

26

nonautonomous systemszo. This simulation should exactly
reproduce the solution of the master equation, since it is
also based on the same postulates from which the latter has
been derived.

To derive the master equation we used the postulates to
arrive at the transition probability Axy(t)0t + 0(0t) that an
event occurs in (t,t+6t) resulting in a transition from y to
x. Instead we now focus on the probability P(T,m;x,t)dr that
the next event occurs in (t+1,t+1+d1) and that it is the m-th
reaction step, given that the system has x molecules of X

at time t;
P(I,m;x,t)d1 = PO(T;X,t) am(X,t+T)dT, (2.53)

where P0(T;X,t) is the probability that no reaction occurred

in (t,t+T) given that the system is at x at time t.
Decomposing the interval (t,t+s+ds) into (t,t+s) and

(t+s,t+s+és) we find that the probability that the m-th step

does not occur in (t,t+s+6s) is given by

P0m(s+és;x,t) = [l - am(x,t+s)és - o(6s)] P0m(s;x,t),

(2.54)

Therefore,

8 _ _ .
Epom(s,x,t) - am(X,t+S) POm(SIXIt)° (2°55)

27

Solving this first order differential equation with the

initial condition P0m(0;x,t) = 1, we obtain

t+T
P0m(T;X,t) = exp - am(x,s)ds . (2.56)
t

Thus the probability PO(T;X,t) that none of the steps occurs

in (t,t+T) is

t+T
P0(T;X,t) = I; P0m(I;x,t) = exp - [ a(x,s)ds ,
t
(2.57)
where
a(x,t) = Z am(x,t), (2.58)
m
and finally from (2.53)
t+T 1
P(T,m;x,t) = am(x,t+T) exp - I a(x,s)ds ). (2.59)
. J

For a given x and t, this is the density function of a joint
distribution in T and m. In order to simulate the process
determined by this density, we first obtain the probability

P1(T;x,t)d1 that the next reaction occurs in (t+T,t+T+dT)

28

irrespective of which step it is.

P1(I;x,t) = Z P(I,m;x,t)
m
t+T
= a(x,t+T) exp - J a(x,s)ds . (2.60)
t

Then the probability P2(m|1;x,t) that the next reaction is
the m-th step given that the next reaction occurs in

(t+T,t+T+dT) is

P(I,m;X.t) = am‘x'tm
PlTr;x,t) a(x,t+I)

 

 

P2(m|1;x,t) = (2.61)

The simulation algorithm consists of selecting a random
number I distributed according to Pl(1;x,t) and another
random number m distributed according to P2(m|r;x,t) when
the system is at x at t. The distribution functions

corresponding to P1(T;X,t) and P2(mII;x,t) are respectively
T

Fl(1;x,t) = J P1(s;x,t)ds (2.62)
0

P (v)1;x,t) (2.63)

and F2(m|1;x,t) = 1 2

">43

V

The range of the function F is [0,1] and the range of F

1 2

is a subset of the above interval. The probability that

T lies between T and T is given by

1 2

29

I2

I P1(s;x,t)ds = F1(T2;X,t) - F1(Il;x,t).
T1

Thus the distribution of I in [0,w] according to P1(T;X,t)
corresponds to the uniform distribution of F1 in [0,1]. Hence
to implement the selection of I, we choose a uniform random
number r1 in [0,1] and set

1 = F-1(r1;x,t) ; (2.64)

similarly we select another random number r2 and choose m

such that

F2(m-1|I;x,t) < r2 6 F2(m|I;x,t) . (2.65)
Having selected I and m, we advance time by I and carry out
the m-th step in the reaction scheme. The whole procedure is
then repeated for the new state. This iteration produces
a random path in time of the system. To get the mean path we
take a large number of these trajectories and average.

The state of a system close to equilibrium is given by
a Gaussian distribution for which the mean is the most
probable value;21 in systems far from equilibrium the mean
trajectory need not be the most probable one. The above
algorithm is suitable for finding mean quantities. An approach
that focuses attention on the most probable quantities and

deviations from them is the path integral approach.22 Let us

30

consider the stochastic differential equation

(I) dXt = - r(Xt)dt + dWt (2.66)
or, equivalently,

(S) dXt = - r(Xt)dt + th- (2.67)
The Green's function P(xf,tflx0,t0) (i.e., the parobability
that the system reaches xf at tf given that it is x0 at to)

is given by the integral of a functional over the class of

continuous functions x(t) subject to the end point conditions

x(tf) = xf and x(to) = x0,namely
tf
P(xf,tflxo,t0) = exp - J L(x,x.I)dI U[x(t)],
t0
(2.68
where
L(x,x,t) = 31mm + r(x(t)))2 -% -dd—x(r(x)). (2.69)

In other words, each path connecting (x0,t0) and (xf,tf) is

( tf
assigned a probability density exp ( - J L(x,x,I)dI and the
to J

total probability is obtained by summing the contributions

of individual paths. The most heavily weighted path is obtained

31

by minimizing the functional

fL(x,x,I)dI subject to the end

conditions. Thus the path withomaximum probability density

is given by the solution of

dzx dzr
.7 = r(x)—+7
dt dx

Similar path integral expressions have been developed23

(2.70)

in

special cases for the most probable path of a process satisfying

(I) dXt = r(Xt)dt + ¢a(Xt)th.

(2.71)

CHAPTER III

A REVIEW OF STUDIES ON THE SCHLOGL MODEL

Chemical reactions as a rule are nonlinear processes
from the thermodynamic point of view; i.e., the rate of
reaction is nonlinear in chemical affinity. Hence chemical
reaction models are convinient for studying phenomena that
are impossible in the linear regime. With this aim, Schldgl
introduced in 1971 two reaction models and showed the

3'4 He also showed the

presence of multiple steady states.
possibility of phase transitions and critical phenomena.

Since then these phenomena have been extensively studied by
various authors using both models, especially the termolecular

model. A summary of the important studies on the termolecular

model is given in this chapter.

3.1 FURTHER PROPERTIES OF THE STEADY STATES

The rate equation for the Schldgl model is

ix = - 3 7- -
dt kZX + k1[A]X k4X + k3[B]. . (3.1)
In terms of the scaled variables
3k2
n = mx, (3.2)
k12[A]2
= _______ 3.
and t8 9k t, ( 3)

32

33

it becomes

 

dn _ _ 3 2 _
EE- — n + 3n Bln + 82

s

= F'(n),
where

B = 9k2k4

1 k12[A]2

2
27k k [B]

B = 2 3

2 k13[A]3

The critical point is given by

(3.4)

(3.5)

(3.6)

1, (3.7)

Consider an arbitrary 82 > O. For this 82 there is a set of

concentration values nS
such that a system described by n =

is at steady state; i.e., F'(nS,Bl)

we get a relation between Bl and ns.

equation can be written as

F'(nS’Bl(nS)) = 0-

and corresponding 81 values 81(ns)

ns. 81 = 81(ns) and 8

"2
= 0. By solving this

Then the steady state

(3.8)

34

Differentiating this, we obtain

'dF' BF' aE' d5

dns ans 881 dnS

 

From the results following equation (2.7) we have

O?

F

W70
S

at any marginal stability point, while

fl... = —n ¢ 0
881

in general. Therefore,

1811. = o (3.9)
dn
5

at any marginal stability point. Since the two marginal
stability points merge at the critical point, it is an
inflection point; i.e.,

dzBl

_—7 = O (3.10)

dnS

at the critical point. From (3.8) we also obtain

35

d2F' 82F' 32F' dB aE' dZB d8 8 dF'
+ 1 + ——— 1 + 1 -
2 2 a as d as d d as d 0
dns ans ns 1 ns 1 ns ns 1 ns

 

 

 

 

2 _

(3.11)
The second , third, and fourth terms contain factors of

del dZB
-—— and

dnS dnS

 

% which vanish at the critical point. Thus equation

(3.11) leads to

j = 0 (3.12)

at the critical point. These relations can be generalized

to multicomponent systems.24
The above properties are possessed by the Schlogl

model when the concentrations of A and B are kept constant.

Escher and Ross25 have recently considered the Schldgl model

with constant flux of A, B, and X, and have studied the number

and stability of steady states. Let the input be a mixture of

these chemicals (of fixed composition) entering the reactor

at rate J and let the output have the same rate J. The

reactor is assumed to be well stirred so that the output

stream has the same composition as the material in the reactor.

If the input stream contains no A, then the system does not

exhibit instability. If the input stream contains only A,

then interesting steady state structures are obtained. For

example, there exists a particular concentration of A

36

in the input stream such that the following situation is
observed: For small values of J the system has two stable
steady states and one unstable steady state. As J increases
two of them disappear leaving a single stable steady state.

At a still higher value of J two more steady states appear,
but both are unstable. With further increase in J one of the
unstable states becomes stable, and eventually a single stable
steady state remains. Escher and Ross also studied a general
case where all species are in the input mixture. In this case
there are two disjoint regions of J where there are unstable
steady states. The instability in the higher J value range is
accompanied by two stable steady states, whereas the lower

range has a limit cycle oscillation in the concentration space.

3.2 COEXISTENCE OF TWO PHASES

Let the Schlégl reactions occur in a reactor which is
not being stirred and let the concentrations of A and B be
kept constant. If A and B diffuse very rapidly but X diffuses
slowly, then the rate equation becomes

Sat-x(E't) = Dvr2x(3:_.t) + F[x(_r_,t)1, (3.13)

where

F(x) = - k2x3 + k1[A]X2 - k4x + k3[B] (3.14)

37

and D is the diffusion coefficient for X. First we seek steady
states4 which are uniform along two spatial dimensions and

vary only along one coordinate r These steady states are

 

 

3.
solutions of
azn 3 2 dU
;;—5 = n - 3n + Bln - 82 = - dn , (3.15)
3
where
k1[A)
xi = ri i = 1,2,3. (3.16)
3/Dk2

For the analysis below, it is useful to note that this is the
equation of motion of a classical unit mass moving in the
potential

4 2

0(n)=-%—+n-B e

NIS
+

1 n (3.17)

if x3 is identified with the time coordinate and n with
position. To determine the steady-state coexistence criterion

for this case we further impose the boundary conditions

n (x3=—00)

ll
:3

and n(x =+w) (3.18)

ll
:3

3

so that the concentration of X takes on two different

38

homogeneous steady state values (n1 and n3) at the extremes

of the r3 coordinate. Equation (3.15) and boundary conditions
(3.18) together describe the motion of a classical body from

one maximum of the potential U(nl) where it was at rest at

t = -w to the second maximum U(n3) where it comes to rest

at t = +m. This motion is possible only if the potential

maxima are equal; i.e.,

U(n1) = U(n3). (3.19)

This condition is satisfied when the three steady states are

related by

n = (n

2 + n3) / 2 . (3.20)

1

Equation (3.20) is Schlogl's coexistence condition. The

homogeneous steady states become

= .. 1..
n1 1 7- 81
h2 = 1
n3 = 1 + ”3'81 (3.21)
and 81 = 82 + 2 at the coexistence point.

Next we seek steady states with Spherical symmetry

varying along the radial coordinate r = (E) and satisfying

39

the boundary conditions

n(r=0)

I)
:3

and n(r=w) (3.22)

ll
:3

Such states can be considered as droplets of the condensed
phase with composition n3 in a medium of composition n1.

The steady state equation becomes

 

'.a_._n—...2l)'.a_r_1_—diji
D 8r2 — 17 Sr dn (3°23)
9kZD
where D' = 2 2 . This corresponds to the motion of a
k A
1

mass in the potential U with dissipation of energy due to
friction. A solution exists if the energy lost is equal to the
potential difference; i.e., if

an) 2

- . J; __
U(n2) - U(n1) - 2D I r' [ r) dr . (3.24)
o

If we assume that the radius r of the droplet is substantially

0
larger than the thickness of the boundary layer over which
the concentration change is significant, then equation (3.24)
leads to the condition for the existence of droplets:

) B1

_ _ ii )IiI - __
8 — 61 2 + r0 20 1 1 3 }. (3.25)

40

Similarly, the boundary conditions

ll
:3

n(r=0)

and n(r=m)

ll
:3

(3.26)

correspond to a steady bubble of composition n in a medium

1

The condition for a bubble of radius r is

of composition n3. 0

(3.27)

Using a linear analysis, Schlbgl 33 3126 have recently
studied the effects of small fluctuations in systems with
coexisting steady states varying only along the r3 coordinate.
When Schlagl's coexistence condition is satisfied, we have

by the transformation y = n - l

2 2 2
£91 = VX y - y(y -y0 ). (3.28)
S

where y0 = 73-81 and t5 and xi are defined by equations (3.3)
and (3.16) respectively. The coexistence solution of this

equation is

i-

y = yo tanh (y0x3//2) . (3.29)

After the change of variable

41

g = tanh (y0x3//2) ,

the linearized equation in terms of the deviation

*
W = (y-y )/y0
becomes
y 2
d _ 0 _ 2 2
d?- - -§— (1 C )DZW + AW)

 

 

where
32 32
A .. +
a 2 a 2
x1 x2
and
v“ = -3—(1-'2) 1+ 9.()L+l) -—“—2 .
i a; 9 a: 1-C2

Introducing the separation of variables

w(xllX21CIts) = Q(xllxz) “(C) r(ts),
we get
A _
DZU - 0,

2
(A+k )Q = 0,

(3.30)

(3.31)

(3.32)

(3.33)

(3.34)

(3.35)

(3.36)

42

and P(ts) , (3.37)

II
(D

with

v = k + (4-02) -——-- (3.38)

The mode with smallest v is given by the associated Legendre

function,

0))
I‘d

 

3C
yo 3x (3.39)

P§(c) = 3(1-C2) =
3

with v = k = 0, and therefore the solution including the small

deviation is

*

"<
I

37200 d _,,( )
dx Y X3

*
y (x ) +
3 Y0 3

 

" *
O
{ l + 6x3 3x3 } y (x3)

*
Y (X3+0X3L (3.40)

where 6x = 3/200 / y Thus this lowest mode called a

3 0'

"Goldstone mode" 27 results in a slight shift along x The

3.
shifted solution is also stationary.

All time-dependent modes that vanish at infinity include

43

the factor

i:
I)

3’
P;(c) = - 3C(1-C2)2

3 sinh (y0x3//2)
= - 2 . (3.41)
[cosh (y0x3//2)]

 

This determines the long-time decay of perturbations along
the x3 coordinate. There is also a continuous spectrum of
solutions P2(;) that oscillate at infinity.

When the coexistence condition is not satisfied, there
is no stationary solution to equation (3.13). However, there
are solutions which move along the x coordinate with constant

3
velocity, c. The solutions are given by

y(t) = y0 tanh [(X3-.CtS)YO/2]’ (3.42)
where
c = /2 [n2 - (n1+n3)/2]. (3.43)

Here again there is a nondecaying Goldstone mode moving with

the same velocity, c. For

1 3 1
2 2 < 3 Y0 , (3.44)

 

44

the long-time regression is controlled by another discrete
mode. When the inequality (3.44) is not satisfied, the
long-time regression is controlled by the lowest of the

continuous Spectrum of modes.

3.3. APPROXIMATE SOLUTION OF MASTER EQUATION

Matsuo28 has developed WKB-type solution to the master
equation in the limit of large system size. As already seen
(equation (2.41)), the master equation can be transformed to

a form with symmetric matrix:

a _
330(x.t) - E sxy Q(y.t). (3.45)

Its solution is

A.

00 t
0(x.t) = 00(x) + 2 pj e 3 0j(x) (3.46)
i=1

where g Qj = ngj. Let us denote the diagonal elements of g
by So(x) and the off-diagonal elements by -S+(x). Then the

eigenvalue equation is

[80(X) - )j] Qj(X) = S+(x-l)Qj(x-l) + S+(X)Qj(X+l).
(3.47)

We look for solutions of the WKB type, namely

45
x
Q(x) = Re C(x) exp J w(y)dy , (3.48)

where C and w are assumed to be slowly varying functions of
x. Substituting this in the equation and its derivative and
neglecting small quantities such as C(x+l) - C(x), we get
expressions for C(x) and w(x). The functional form of the

solution depends on which of the following three situations

obtains:

a) - 28+(x) < So(x) - A < 28+(x)
b) So(x) - l > 25+(x)

c) So(x) - A < - 23+(X). (3.49)

We obtain quantum conditions on the eigenvalues by matching
the solutions in adjacent regions. For the Schlogl model the
state space can be divided into five regions separated by
x1, x2, x3 and x4, in which conditions a, b, a, b, a are
satisfied in order. If the middle region is macroscopic in
size, then the two a-regions can be considered independent
and the connection conditions are simplified. In this

approximation the eigenvalue 1 must satisfy one of the

following conditions:

46

x4

%Jq(x)dx = 9.2 + 35 ’ (3.50)

X3

where 21 and 22 are nonnegative integers and

_ So(x) - A
q(x) = cos . (3.51)
23+(x)

g...)

 

From this condition we can obtain the density of eigenstates

d1 d2

_ ___1_ __7;.

It is found that in the limit 81 = V-1 —9 0, the density of
eigenstates near A = 0 diverges at the marginal stability

points and the critical point. This divergence is given by

lim -1/4 . . . .
€1,A-+0 0(A) % A (marginal stability pOints)
-1/3 . . .
m A (critical pOint). (3.53)
Since the eigenvalues accumulate near A = O, the decay to the

stationary distribution becomes infinitely slow. This is called
"critical slowing down".

When there are three distinct steady states, it can be
shown that the first decay mode has one node and the second
has the largest of its three peaks in the vicinity of the

unstable steady state. Thus we can find the relaxation times

47

from the unstable steady states to be

1

A = [ k (XZ-Xl) (X3'X2) ] (3.54)

2

where k2 is the rate constant of the second step; similarly,

an expression can be obtained for the relaxation time from

the metastable state to the stable one.

3.4 RELATIVE STABILITY

We have already seen Schlogl's coexistence condition.
When this condition is not satisfied, it is clear from equation
(3.42) that the system will eventually become homogeneous.

If n is closer to n than to n then the final state will

2 1 3'
1; otherwise, it will be n3. This can be regarded as a
criterion for relative stability of the two stable steady

be n

states. Other proposed criteria are discuSsed below.
Noyes29 has considered a general reaction scheme with a
single reactant R producing a single product P,

R P ,

-—>
‘—
occurring in a continuous flow reactor. The Schldgl scheme with
constant flow boundary conditions is included in this class

of models. The rate equation is

dR_ __ ..
HE — kO(R0 R) v(R), (3.55)

48

where v(R) is the net rate of reaction, k0 is the rate of

flow of matter through the reactor and R0 is the concentration
of R in the input stream. It is assumed that the reactor is
efficiently stirred at all times; then the composition of

the output stream is the same as the composition in the
reactor itself. Functions for v(R) can be selected so that

there are two stable steady states Rd and R for a range of

8

k0 values. Let us prepare two adjacent flow reactors with

identical external conditions including the flux rate kORO,
but at different stable steady states. In order to determine
the relative stability of the two states, Noyes has proposed
making small holes in the wall between these reactors thereby
allowing hydrodynamic mixing with a rate constant kx. For
small kx the concentration of R in the two tahks*WillFbe
perturbed to the new values Ra' and R '.

8
Noyes has given a criterion, namely

0)

V

k0 + [—E (3.56)

g..__J
w
)l
O
\

at which a slight increase of kx is claimed to cause the RG'
state to disappear and the combined system to settle with
both reactors in the RB state indicating that R8 is more
stable. If

k0 + [33] = o (3.57)
R '

49

happens first as a consequence of mixing with steadily increa-
sing kx, then the higher stability is assigned to R0‘ As will
be seen in chapter VI, this criterion and the proposed mixing
experiment are unsuitable for determining relative stability.
Another criterion for relative stability is obtained by
taking the thermodynamic limit in the stationary distribution

30-32 To do this, we should first

of the master equation.
extract the system size dependence from all parameters and
take the limit as V —> m, X —9 w while X/V remains fixed.

The rate constants have the volume dependence
c. = k. n.! V , (2.20)

where ni is the molecularity of the i-th step in X and mi is
the total molecularity of the i-th step. The master equation

for the concentration x = X/V becomes
a — - -
fiP(X,t) - 3+(X €1)P(X Ellt) + a_(X+€l)P(X+€1rt)
- [a+(x) + a_(x)) P(x,t). (3.58)

Its stationary solution is given recursively by

a (x—E )
_ + .1- _
Ps(x) - Ps(x 81). (3.59)

a_(x)

where

50

81 = l/V’
b = B/v,
k1
a+(x) = Tbe (x-El) + k3bV 1
k2
and a_(x) = 7?VX(X‘€1)(X-2€l) + k4V. (3.60)

Ps(x) has local maxima at x1 and x3 and a local minimum at x3,

where x1 < x2 < x3 are the steady state concentrations. First

we evaluate

lim PS(X)
V—>°°_——— f0r0<x<xl;

Ps(x1)

in terms of the function

a (z+el)
0(2) = —2—————— . (3.61)
a+(2)
Ps(x)
ln —_T_—— = ln G + 1n G + ... + In G _
PS X1) X x+€1 X1 81
X1
= V Jln G(z)dz (3.62)
X

for large V. The integral representation is verified by writing

51

the integral as a Riemann sum with A2 = 81. Since Ps(z) is
increasing in 0 < z < x1, G(z) < l in this range and the
integral is negative. The required ratio thus goes to zero
in the thermodynamic limit. Following the same procedure

it can also be shown that

. P (x)
lim 5 =
v——>w P (x1) 0 for X1 < X < X2
5
. P (x)
lim 5 = m
and V m PS(X3) 6XX3 for x2 < x < . (3.63)

Thus the stationary distribution becomes sharply peaked at
either or both of the stable steady states as the volume of
the system increases without bound. Finally the relative

stability is determined by

X3
lim 3513;: = lim ex V ln G(z)dz (3 64)
V_—’m Ps(X3) V-—>w p '
X1

which is zero, one or infinity depending on whether the
integral is negative, zero or positive. Thus the equistability

condition is given by

X3
I 1n G(z)dz = O. (3.65)
X1

33

Procaccia and Ross have studied relative stability

52

based on the nonequilibrium thermodynamics developed by

Levine.34’35

Let us consider an open system whose microscopic
states X are numerous compared to the macroscopic information
available in terms of X. Then the distribution over the

microscopic states must satisfy

ixiP(X,t) = xi(t) for i 0,1,...,M (3.66)

X

where x0 = l imposes normalization. These equations are
insufficient to determine P(X,t). Of all the distributions
satisfying (3.66) the one with minimum information content
is postulated to describe the system, in the absence of any
other information. This distribution also maximizes the

entropy function

SIP] = - X P(X,t) 1n [P(x,t)/q(x)] , (3.67)
X

where g(X) is the degeneracy of state X. When X specifies the

number of molecules of each type, the degeneracy is

LG
)4
ll
——"‘“\
">43

11)

M
x.)1/|—|'x.1 (3.68)
\i j=1 J

The distribution thus obtained is

53

thz

P(x,t) = q(x) exp - ur(t)xr . (3.69)

3 0

where ur(t) is the Lagrangian multiplier corresponding to
the r-th constraint (3.66). The function

u
e 0 = 1 9(2) exp - _Z ui(t)xi (3.70)
X

is called the partition function. From the conservation of

total probability,

800 M Sui _

___ = _ Z Xi(t)-__" (3.71)

at i=1 at
In the following analysis the equilibrium state is used as a
reference state. (This state is attained in the SchlBgl model
only when the flux constraints are removed.) Denoting
equilibrium quantities with a superscript e, the entropy
deficiency is defined as

P(x,t)

K[P(t),Pe] = 2 P(X,t) in 'E'—'— ° (3.72)
X P (X)

Substituting (3.69) into (3.72) we find

M
x(t) = - 400(t) - .2

1-1Aui(t) Xi(t) (3.73)

where Aui(t) = ui(t) - pie. In analogy with classical mechanical

54

equations, We look for a function L({Xi},{Xi}) such that

 

 

 

 

 

d BL
86% [axi X k . (3°74)
j#i' k
This equation is satisfied if L is defined as
M dXi
L = EAUi —d—‘-t_ . (3.75)
To prove this, we first note that
SE - - dAuO - g x dAui - g A“ dXi
at dt i=1 1 dt i=1 1 dt
M dX.
= - {Mi—J. = -L (3.76)
i=1 dt
using (3.71). Now
31- _ - 3. __3K 1. .1
(3Xi]x i - dt[8xi]x i th“i (3°77)
j#i' k j#i' k

from (3.73).
The assumption of maximum entrOpy at all times has led
to the equation of evolution (3.77). This evolution equation

is now used to study the relative stability in the Schldgl

55

model. First since g(X) is a multinomial coefficient, we may

write the partition function as

 

‘11. N
e l] , (3.78)

2}
If
(D
H
(D
2
II
II M3
9*:
H.
U)
£1)
0)
U)
(:1
B
(D
0..
r1-
0
H
(D
B
D.)
H.
:3
O
O
:3
U)
('1'
D)
:5
‘1-
H.
:3
w
H.
3
(D

BY

-——— = - X . (3.79)

Differentiating (3.78) and comparing with (3.71),

du. dX.
—-1- = -3(--—- (3.80)
dt i dt

The evolution equation for the Schldgl model becomes

3.3-33
th 8x

SAuA dA BADX dX BAuB dB

 

 

 

= — + — + —- ’ (3°81)
8X dt 3X dt BX dt
dA dB .
where 5? and HE are changes due to the reactions alone.
Suppose we start the system at the stable state X at t = -w.

1
If the other steady state is more stable, it is expected that

56

the system will evolve to X given enough time. This transition

3
will occur only if the evolution equation is satisfied.

Multiplying by 3% and integrating from t = -w to t = w gives
00 2 00 X3
1 dX _ 3L dX _ 3X
.. (§[_d—t] dt -— {-5-}? aft- dt - I 3x dX (3.82)
-00 —00 X1

under the hypothesis that a transition occurs, i.e., X(t=m)=X3.
Since the left hand side is negative semidefinite, X3 is more

stable than X1 only if the integral on the right hand side
is negative. If it is positive, then the reverse transition

can occur. Thus the equistability condition is

8L _
I 3? dx — 0. (3.83)
x

Relative stability can also be studied using first passage
times based on the master equation.36 We consider the mean
time Iu of first arrival at X3 if the system starts from X
at time zero and we restrict attention to single variable

1

systems such as the Schlégl model. Let vi(x) be the average

number of x-exil transitions executed and tu(x) be the average

 

time spent in state x during the passage from X1 to X3. Then
from equation (2.60), we have
t (x)
u = -—l——-- (3.84)
VAX) + \)_(X) a(x)

for the mean residence time at state x and from (2.61)

57

 

 

vt(x) '_ ai(X)
_ (3.85)
0+(x) + v_(x) a(x)
Solving these equations for tu(x), we get
vim)
tu(X) = m. (3.86)

1

Since every trajectory under consideration has exactly one

X3-1 -9 X transition,

3

0+(X3-1) = 1. (3.87)

Every trajectory starting from XI must also satisfy certain
conditions in order to reach X3. At a state x < X1 it must
execute an equal number of x-éx+1 and x+1-9x transitions and
at any state between X1 and X3 it must execute one more
x-+x+l transition than the reverse. Since all trajectories
satisfy these conditions, so does the average.

v+(x-l) = v_(x) + h(x-Xl) for x < X3 (3.88)

where h(x) is the Heaviside step function. Thus we have a

recursion relation for tu(x):

a_(x) h(x-Xl)
tu(X-1) = m tu(X) + m , X < Xz‘l
(3.89)
with
1

t (X -1) -—-——.
u 3 a+(X3-1)

58

The mean first passage time is given by summing the mean times
spent in all accessible intermediate states.

X3-1

I = Z t (x), (3.90)
x=0 u

Similarly we can derive for the downward transition

1 = X t (x) (3.91)
d x=X +1 d '
1
where
t (X +1) 1 D
d 1 a_(X1+1)

and in general

a+(x) h(X3-x)
td(x+l) m td(X) + m° (3.92)

Then we can define relative stability based on the difference
Iu - Id. If this is negative, the upper steady state X3 is
more stable and conversely. Thus the equistability criterion

is

I = Td' (3.93)

A more rigorous and direct derivation of the same expression

59

for the first passage time from the master equation is
available in reference 37.

We have several coexistence conditions based on
different criteria for relative stability. Schlogl's criterion
includes the effects of diffusion but neglects fluctuations,
whereas the two thermodynamic criteria (equations (3.65) and
(3.83)) and the kinetic criterion (3.93) are based on
descriptions neglecting spatial fluctuations. These two kinds
of criteria need not agree. In an intermediate situation it
is possible38 that a local fluctuation large enough for a
phase transition arises due to the chemical reactions in the
system but it diffuses away rather than causing a transition.
On the other hand, when a slowly moving boundary is predicted
by Schlbgl's criterion, one may expect several dynamic patches
of the two phases due to local fluctuations. Thus the actual
coexistence point should lie between the points predicted by
Schldgl's criterion and the thermodynamic criteria. Local
fluctuations in nonequilibrium system can be investigated
using the multivariate master equation described in the next

section.

3.5 MEAN FIELD THEORY AND MULTIVARIATE MASTER EQUATION
APPROACHES TO INCLUDE LOCAL FLUCTUATIONS.

The master equation based on uniform concentrations of
the species in the reacting mixture is called the birth and

39

death equation. Malek-Mansour and Nicolis have advanced

60

several criticisms of this method. They have stated that
the prescription for calculating transition probabilities
treats the reactive encounters of molecules as equally
likely irrespective of the position of these molecules in
the reaction vessel. This statement is simply not true. What
appears in the master equation is the probability that a
reactive collision will occur somewhere in the reactor, which
is different from the probability that a given set of molecules
will undergo a reactive collision. This point has been amply

40

clarified by Gillespie. Another criticism is based on a

study of the apparently pathological reaction scheme

k
C+Y—i>2Y
k6
ZY -__—'> Do

The deterministic analysis gives an unstable steady state

Y = 0 and a stable steady state,

Y = k5C/k6' (3.94)

but the stationary distribution of the birth and death

master equation is

P (Y=y) = 8 . (3.95)

According to equation (3.95) the system settles into the

61

state predicted to be unstable by the deterministic analysis.
From this result Malek-Mansour and Nicolis concluded that the
transition probability calculated using the birth-and-death
model is wrong; this conclusion is based on the expectation
that the deterministic analysis always gives the correct
stability. A careful consideration of the reaction scheme
reveals that the result obtained from the master equation
represents a genuine phenomenon consistent with the assumed
mechanism.41 In the vicinity of the null state, the first
reaction is much more likely to occur than the second due to
the presence of C. Nonetheless it is possible that the few
remaining Y molecules will react, eliminating Y from the
system. Once the system has reached the null state, Y
molecules can no longer be produced. This means that the null
state will be eventually reached due to fluctuations and it
is the only asymptotically stable state. Thus the model system
is expected to behave as predicted by the birth-and—death
master equation. It is the deterministic equation that gives
wrong results due to the neglect of fluctuations.

However the master equation does have limitations, as
correctly pointed out by the authors mentioned above.39 They
have developed a mean field theory to remove the limitations
arising from neglect of spatially inhomogeneous fluctuations.
Although the birth-and-death master equation includes
fluctuations in particle number, it assumes the concentration

of the chemical species to be uniform throughout

62

the reactor. This need not be the case when stirring and

diffusion are not efficient. Hence a realistic system can

have properties such as a finite correlation length which

is a measure of the effective wavelength of local fluctuations.
The mean field approach starts by considering a subvolume

AV which is comparable in size to the cube of the correlation

length. The birth-and-death description is assumed to be valid

within this subvolume. Then the probability distribution for

the entire system P(xi,xo,t) evolves according to

_3_

8tP(Xi'XO't) = Rch(AV) + Rch(V-AV) + T(AV,V-AV),

(3.96)

where xi is the number of molecules in AV and x0 is the number
of molecules outside AV. Rch(AV) is the reaction probability
term arising from birth and death processes in AV and T arises
from the transport of matter across the boundary of AV.
Assuming that the transport of heat is very fast (so that the
system remains isothermal), T may be calculated using the
position and velocity distributions from the kinetic theory

of gases. Next the mean field assumption is introduced and the
evolution equation is summed over the external variables, xQ

and V-AV. The system is required to be homogeneous on the

average, i.e.,

< xO/(V-AV)> = < xi/AV > . (3.97)

63

Moreover, the particle distributions in AV and in V-AV are

assumed to be uncorrelated, so that

P(xi.xo,t) = PAV(xi,t)PV_AV(xO,t), (3.98)

Using these relations, a closed equation is derived for the

evolution of the probability distribution in V:

a _ -
3E PAV(Xi't) - Rch(AV) + D<X>[PAV(X1+1't) PAV(Xi't)]
+ D[(Xi+1)PAV(Xi+1't)—xiPAV(Xi't)1'

(3.99)

In equation (3.99) D is the effective frequency of diffusion
passage of particles across the boundary of AV and is related

to the diffusion coefficient D and mean free path by

_ D
D - mean free path x coherence length of fluctuations

(3.100)
If the characteristic time scale of macroscopic changes

due to chemical reaction is much longer than that of diffusion,

then there is an expansion parameter.42 For example, in the
Schlogl model if
c k
2 2
s = - = ——————- (3.101)
2 6D (AV)20

is small, then the "nonlinear" master equation becomes

64

a ' _ '1 _ y _ I
5ETp(X,t ) - Rch + 82 {<X>[P(X 1,t ) P(th )]

+ (x+l)P(X+1,t') - XP(x,t')} , (3.102)

c k

 

where t' = —-2-t = 2 2t and
6 (AV)
RCh = ii { a+(x-1)P(x-l,t') + a_(x+1)P(x+1,t')
- [a+(x) + a_(x)] P(x,t') }. (3.103)

We can seek a solution of the form

P(x,t') = P0(x,t') + E (x,t') + ... . (3.104)

2P1

Since P(x,t') must be normalized independent of e , we have

2

X Po(x,t') = 1 and 2 Pn(x,t') = 0 for n>1 . (3.105)
X X

The solution is obtained by use of a probability generating

function defined by

F(s,t') = sXP(x,t'). (3.106)

"MB

x O

F(s,t') can also be expanded as a series in 82

F(s,t') = Fo(s,t') + s Fl(s,t') + ... ,

2

65

subject to the conditions

F0(s,t') = l and Fn(s,t') = O for n>l .

s=1

 

(3.107)

Then the evolution equation becomes

8
as

 

———F(s,t') = R + €2-1(S-1){<X>F(S,t') - F(s,t')} ,

(3.108)

sXRC . By expanding <x> in powers of 62

<x> = <x> + e <x> + ... (3.109)

and then equating the coefficients of powers of 82 in equation
(3.108), a set of evolution equations for Fn(s,t') are
obtained; these can be solved successively.

As an example, let us determine the stationary distribution
for the Schlagl model, which is expected to be different from
that obtained from birth and death master equation. For the

Schlbgl model R is

3 2

2 d F 2 d F dF
R = (l-s) 5 ——§ - bls ‘-—Z + b3 —— - b2F (3.110)
ds ds ds

where

66

 

k1 3c1
bl = 'k— bAV - — bAV,
2 2
k 6c
b2 = k2 (AV)2 = C ;
2 2
k 6c
b = —3- b(AV)3 = ——3 bAV.
3 k2 c2

At steady state, then

2 d3F 2 d2F dF dF
E s ——— - b 5 ——5 + b —— - b F + —— - <x>F = 0.
2 ds3 1 ds 3 ds 2 ds
(3.111)
Equating the coefficients of 820 gives
dFO
———-- <x>0F0 = 0. (3.112)
ds

Solving equation (3.112) subject to the condition F0(s=1) = l,

we get

F0 = e . (3.113)

Equating the coefficients of €21 gives

67

 

2 d3FO 2 szo dF dF
s - b 5 1——1r + b ——— - b F + ——— - (X) F
ds3 1 ds 3 ds 2 0 ds 0 1
.. < =
x>1F0 0 . (3.114)

Setting s=1 and using (3.106), we obtain a moment equation

<x(x-1)(x-2)>o - b1<x(x-1)>O + b3<x>0 - b2 = 0.
(3.115)
Since F0 generates a Poisson distribution, <xn> = <x>n and
therefore <x>O is the solution of the macroscopic rate equation.
<x> 3 - b <x> 2 + b <x> — b s 0 (3 116)
0 l 0 3 0 2 ' °

Thus for the Schlogl model three steady states may exist

when e -90. To find F1(s), we substitute

2

F1(s) = F0(s)wl(s) (3.117)

in equation (3.114) to obtain

(3.118)

with 01(1) = 0. Solving (3.118) yields

68

_ _ _ l 3 3_ l 2 3_
01(5) - <x>1(s 1) 3<x>0 (S l) + 3b1<x>0 (s 1)

- b3<x>0(s-l) + b2(s-l) . (3.119)

To find <x> we equate the coefficients of 622:

1)

 

 

 

2 d3F1 2 d2El dPl sz
S 3"b1s 2+b3—"b2F1+_'<X>0F2
ds ds ds ds
-<x>1F1 - <x>2F0 = 0. (3.120)
dF2
Setting s=1 and noticing that-——— = <x>2, we find
ds s=1
d3P1 d2)?1 dFl
- b -_7T + b ——— - b F = 0. (3.121)
ds3 1ds 3ds 2 1 s=1

Substituting the expressions (3.117) and (3.119) for F1 and (1

we finally obtain

2
<x> = - 6b3<x>0 + <x>0[6b2 + 2b1b3 2b3] 2b2(b1 1).

1 2
3<x>O - 2b1<x>O + b3

(3.122)

This procedure clearly breaks down at the marginal stability

points where the denominator in (3.122) vanishes.

69

A more satisfactory, but more difficult treatment is the
multivariate master equation approach.43 In this approach too,
the entire system volume is divided into a number of cells
of equal volume AV. A birth and death process is assumed to
occur in each cell, in addition to migration between adjacent
cells. In contrast with the mean field approach, the equation
of motion for the detailed distribution P(§,t) is considered,
where the j-th component of § is the number of g molecules
in the j-th cell. P(§,t) is a n-dimensional joint distribution
when there are n cells and only one chemical component is
allowed to vary. The multivariate master equation for the

Schlogl model is

5%P(§.t) E a+(xr-1)P(§-§r,t) + a-(xr+1)P(§+§r't)

- [a+(xr) + a_(xr)] P(§.t)

+ D g [(xr+1)p(§+§_Q-_§r,t) - er(§,t)]
(3.123)
for all g in Zn where Z = { 0,1,... }, fir is a row vector

with r-th component unity and all other components zero,

1 runs over all nearest neighbours of the r-th cell,and.D is
the diffusion coefficient of x. The generating function

for the joint distribution is

X
F(Ert) = X W Sr r P(gg.t) (3.124)
X r

70

Since the variables in the Schlogl model can be scaled leaving

only two parameters we may define the parameters w and w' by

k 6c

w = -i-- 3 = -———3—3 - 3, (3.125)
k2 c2(AV)
k3b 6c3b

w' = -——— - 1 = - 1, (3.126)
k2 c2(AV)2

and we may select [A] such that

kllA]

 

= 3 (3.127)

k2

without loss of generality. Then w=w'=0 is the critical point.

The steady-state equation for the generating function becomes

3 2

 

 

 

2(1-sr)sr2 1: 8 E; - ii-jij;
r (AV) 38 AV as
r r
3F
+ 2 (1-3 )[(3+w) - (1-w')AV F]
r r 3s
r
8F
+ D X 2(sg-sr)fi = o. (3.128)
r 2 as

r

71

Nicolis 23 31 have assumed that the steady state distribution

is of the form

F = 2 exp AV E (Sr-1) , (3.129)
r

and that the volume of the cells is sufficiently large that
the solution can be obtained as a perturbation expansion in

a = (AV)’l << 1 . (3.130)

Considering the approach to the critical point along the line
w=w', we scale the quantities w and D near the critical point.
Since the generating function and its derivatives are needed

only at gel, we also put

sr = 6 gr + ... 0 < e < 1 (3.131)
along with

w = ijl + ... (3.132)
and D = efDl + .. . (3.133)

Substituting equations (3.129)-(3.133) into equation (3.128)

shows that the steady state equation has the leading terms

72

 

3
3 Z
{gr 2(18) 3+Z4€22612+2_€€laarz
ear r r 5r
+ZiarfleD[—3€——§8€—r]+...=o. (3.134)
r

The applicability of different kinds of theories near the

critical point is determined by the relative importance of

these four terms. Thus, for

f > max {j,2(e—1),2e-1}

and j = 2(l-e) = 2e-l, (3.135)

the birth and death description neglecting spatial fluctuations

is valid. On the other hand, for

j < min {j,2(1’€),2€-1} , (3.136)

diffusion dominates and chemical reaction effects appear as

perturbations. All terms have equal importance if

e = 3/4 and j = f = 1/2 . (3.137)

The resulting third-order equation in this case has the

solution

73

oo oo

Z({€r}) = [... J {der} [eXp( Eirer)] R({er}),

—oo -(X>
where

2 4 D

8 9
r 1 2
+ 4 + 7? §(62 er) '

_£_
1 2

R({6r}) exp -%2 w

r

(3.138)

The functional in the exponential is the Ginzburg-Landau

potential. If
e < 3/4 and 2e—1 = j = f, (3.139)

the third order term can be dropped and consequently the 6r4
term is missing in the solution. This corresponds to the

Gaussian approximation.

3.6 CRITICAL PHENOMENA
The steady state structure of the Schlogl model is similar
to that of a van der Waals gas or a ferromagnetic material.
A close analogy exists between the parameters of these
systems. The transition between multiple steady states of
the Schlogl model is similar to the gas-liquid transition.

For the Schlogl model the pump parameter B plays the role

74

Of pressure and a change in the ratio of rate constants affects
the shape of the X versus B curve as a change in temperature
affects pV curves. The order parameters are the concentration
of X and the density respectively. A less complete analogy
can be drawn with the magnetization of a ferromagnetic material
where the pump parameter and the rate-constant ratio correspond
to external field strength and temperature respectively.
Because of this similarity, the observed critical phenomena
in the case of gas-liquid and magnetic systems44 are expected
to be manifest in the Schlogl model also.

Inhomogeneous spatial fluctuations are very important
in critical phenomena. So, the birth-and-death description is
not sufficient for this study. Including the diffusion process

the deterministic rate equation can be written as4s'46

 

 

3gT-xqflt) = - x3 + 3Ax2 - (3+6)A2X + (1+6')A3 + szx,
(3.140)
where A, 6 and 6' are defined by
_ le _ clB
A - 3"— - -C__ ’
3 2
9k k BZV2 6c 0 B2
3 + 2 4 _ 2 4
5 _ ... )
k 2 C 2
1 1
27k22k3V2 6c22c3
1 + 6' = = -————-) (3.141)
3 2 3
k1 B c1

75

and t' = -— = —2t. X(£,t') is related to its Fourier

V 6

components Xq(t') by

X(£,t') L_d/2 2 x (t') eifl°£
g ..
L/2 L/2 -i .r
v _ 'd/2 o
and X (t ) L dr1 drd X(£,t )e ,
-L/2 —L/2

where L is the side of the d-dimensional cubic reaction vessel.

In terms of the Fourier components, the rate equation is

3 . _ 2 2 -d/2
SETXg(t ) - [(3+6)A + q DIXq + L BA g XS‘EXE

_ -d d/2 , 3
L 2 Z. xkxk,xq_k_k, + L (1+5 )A éq'o
13 ————— --
(3.142)
using
L/Z iklr1
e dr1 = L. (3.143)
-L/2

It is difficult to study the corresponding stochastic equations

in completely general form. The first improvement on the

76

deterministic equation is the Langevin equation obtained by
simply adding uncorrelated uniformly fluctuating ("white")
noise terms na to the deterministic equations; the averaged
mg and correlations satisfy

‘

< t > = O
nq( )

and < (t n (t > = 2F 6 6 t -t ) 3.144
n31 1) 32 2) 31 32 ( 1 2 . ( )

Making the change of variable to o

o X /A - 6 (3.145)

s a 3:9,

and adding the noise terms transforms equation (3.142) into

a _ _ 2 . d/2 ,_
at"°q [6+q D ]oa + L (6 5)Cg,9
+ L“d 2 2 ok ok, 0 _k_k. + n (t"), (3.146)
EE- - 51-— 3
u 2 2
where t = A t' and D' = D/A .

The corresponding Fokker-Planck equation giving the time
evolution of the probability distribution for {o } is a

second-order partial differential equation. Its stationary

solution is

PS({OE}) = N exp[-F({Gk})], . (3.147)

77

where N is a normalization constant and

F({ok}) = r’1 2 [6/2 + D'k2/2]|Ok|2
_ k _
+ % L-d X Z Z c o g o
k k' k" .12 £0 .12" _l<-_£l_-]SII
+ L d/2(6-6')o . (3.148)

This has the form of the Ginzburg-Landau Hamiltonian in the
Ising model of a ferromagnet.

Many physical properties diverge as the critical point
is approached. The strength of this divergence is characterized
by critical exponents. For example, when the critical point

6 = 6' = 0 is approached along the line 6 = 6'

<0 > m (-6) for 6 = 6' < 0 (3.149)

where Y1 is a critical exponent. Similarly other exponents

are defined by4s’46

1/Y2
<oo> N (6') for 6 = O, (3.150)
2 'Y3 _ ,
<00 > N 6 for 6 - 6 , (3.151)
‘Y '2Y
and <0 0 > m 6 3 D(q26 4). (3.152)

q “q

These exponents can be calculated using renormalization

78

group techniques.

Renormalization is a coarse graining transformation
followed by a change of scale.47 Suppose we have a joint
probability distribution over the microstates of the system.
From this we construct another distribution whose spatial
(and time) resolution is lower. To do this we have combined
several microstates together losing details within a single
collection. Let us parametrize the extent of this coarse
graining and scale change by s. Then the set {Rs:s>1} of
all such transformations is called the renormalization group
(even though it does not satisfy all the axioms of group
algebra). The basic hypothesis of renormalization group theory
is that at the critical point all length scales are equivalent;
that is, the system behaves identically at all levels
intermediate between microscopic and macroscopic scales
except close to the microscopic scale. This principle is an
observed fact in magnetization and condensation processes.
Accordingly the probability distribution corresponding to the
critical point becomes invariant asymptotically under a
renormalization transformation; i.e.,

lim

S_,m RSP = P* (3.153)

where P* is a fixed point of Rs for all s. Asymptotically the
distribution corresponding to a point in the neighborhood of

the critical point (but not on the critical surface) moves

79

away from P*. The functional dependence of this movement

on the parameters of the system can be related to the rate of
divergence of physical properties. Thus the critical exponents
can be calculated.

When we apply a renormalization to equation (3.148), we
find that the result has the same form with the coefficient
sq"-d appearing in front of the quartic term. Since we are
interested in large 3 behaviour, we can neglect the quartic
term in the potential when 4-d < 0 and obtain a Gaussian
distribution. The fixed point of the class of Gaussian
distributions and the class of Ginzburg-Landau distributions
are identical only at d = 4. The Gaussian distribution is
valid for d > 4. When d decreases from 4, this approximation
becomes increasingly inaccurate. Thus dC = 4 is a critical
dimension. This suggests a perturbation expansion in a parameter
with exponent 4-d. If the quartic term in the potential is
treated as a perturbation, then an expansion is obtained in
powers of 683 where 83 = 4-d. Thus even when d < 4, there is
a region (designated "the critical region") around the critical
point outside which the Gaussian approximation is valid. As d
increases the critical region shrinks and finally disappears
when d reaches 4. The Ginzburg criterion qualitatively governs
the size of the critical region.

From the perturbation expansion we obtain the following

critical exponents:

80

Y1 = l - 23 +62%; f 0(533), (3.154)
Y2 = 3 + .3 + @2532 + 0(533), (3.155)
Y3 = 1 + €53 + 3%632 + 0(833), (3.156)
and Y4 = % + f%€3 + 142 832 + 0(833). (3.157)

In equilibrium critical phenomena the heat capacity
diverges at the critical point. A property analogous to the
heat capacity exists in nonequilibrium systems. The entropy
(S) and entropy deficiency (K) functions defined by (see

equation (3.72))
S({Pi}) = - E pi 1n 8. (3.158)
e _ e
K({Pi},{Pi }) - E Pi ln (pi/pi ) (3.159)

(in appropriate units) can be viewed as the averages of the

random variables

b. = - ln Pi (3.160)
= — = e

r b. bi ln (Pi/Pi )/ (3.161)

called the bit number and relative bit number respectively.

81

Roughly speaking, bi is prOportional to the number of binary
digits neccessary to represent Pi' In information theory,
-S is called the Shannon information measure and K is called
the Kullback measure of information gain. In equilibrium
statistical thermodynamics, the probability distributions are

of the form (cf. equation (3.69))

Pi = exp — Z Aurxir ]. (3.162)
r

The heat capacity of a system is related to the variance of

energy in the canonical distribution. Since the bit number

corresponds to the energy of the microstates, (e.g.,
b. = - ln Pi = -——- (3.163)

for a canonical ensemble), a generalized heat capacity can be

defined in terms of the bit number variance, as

_ 2 2
Cr - (<b > - <b> ) / Aur . (3.164)
The divergence of bit number variance of the Schlbgl model
has been studied by Schlégl.48 When diffusion is effective and
6 = 6', the stationary distribution for the homogeneous mode
has the form

2

P(x) m exp [ - V (x +0)2 ] (3.165)

82

where o is a derived constant. When the bit number variance
is calculated from this distribution and plotted, it is found
that it dramatically changes at the critical point and that

this change approaches a divergence as V -+ w.

CHAPTER IV

A STUDY OF ASYMPTOTIC RELAXATION IN THE SCHLOGL MODEL USING

EIGENVECTORS OF THE TRANSITION MATRIX.

It is well known that deterministic methods of chemical
kinetics are not adequate for treating reacting systems

exhibiting chemical instability.8'49_52

In particular
stochastic methods are needed to analyze transitions between
multiple steady states in such systems. The Sch15gl model
has been widely used for studying phenomena associated with
multiple steady states. SchlESgl4 analyzed the model
deterministically and showed that for certain sets of parameter
values the model has two stable steady states and an unstable
steady state (section 2.1). Matheson 33 215 compared the
stationary distribution from a stochastic treatment with the
deterministic steady state results (section 2.2). Procaccia
and Ross33 have discussed the question of relative stability
of the two Stable steady states (section 3.4). Gillespie9
studied the time-dependent features of this model using first
passage times (section 3.4). Here the time-dependent probability
distributions are displayed for long times at which transitions
from one stable steady state to another are expected to occur.
The stochastic time evolution of the system is described
by a birth-and-death master equation (section 2.2). Since it
is a first-order linear equation, its solution is in principle

known in terms of the eigenvalues and the eigenvectors of the

83

84

transition matrix8 (section 2.3). A criterion for the validity
of the Kramers rule is given and results for the Schlégl model
are presented in this chapter.

Creel and R03353 have reported experimental observations
of transitions between multiple steady states of an NOZ/NZO

4
mixture illuminated with an argonéion laser.

4.1 THE SCHLOGL MODEL AND ITS STOCHASTIC FORMULATION

In the Schlagl model the chemical reactions

occur in a homogeneous Open system. The particle numbers of
species A and B (denoted by A and B) are assumed to be kept
constant by contact with external reservoirs or by appropriate
feeding into or removal from the reactor. Thus X is the only
variable, while A and B are external parameters. For
convenience we will also assume the numbers of A and B to be

equal. When the reaction parameter n defined by

n = k—E_ (4.1)

exceeds a critical value “c = 9, there is a range of B values

85

for which the system has three possible steady states,

X1 < X2 < X3, of which X1 and X3 are stable and X2 is unstable.
Outside this range of B and below the critical point, there is
only one stable steady state. The boundary separating these
two regions consists of the marginal stability points.

In the stochastic formulation X(t) is viewed as a random
variable taking integer values in {x:0<x<w} The state of the
system at a specified time t is then given by a probability
vector 2(t) whose x-th component is the probability that
X(t) = x. The dynamics of the system are described by the time
evolution of P(x,t) [ P(x,t) E Px(t) ] according to the

master equation,

a -- — —
fipm't) - a+(x 1)P(X 1,t) + a_(X+1)P(x+1,t)
- a(x)P(x,t), x > 0, (4.2)
where
C1
a+(x) = 7TBx(x-1) + c3B, (4.3)
C2
a_(x) = 7?x(x-1)(x-2) + c4x, (4.4)
and a(x) = a+(x) + a_(x). (4.5)

The parameters ci are related to the rate constants ki by

86

(4.6)

where ni is the molecularity of the i-th step in X, mi is

the total molecularity of the i-th step (e.g., n1 = 2, = 3),

m1
and V is the volume of the system. The stationary distribution
Ps(x) is obtained by setting the right-hand side of the master
equation equal to zero. Ps(x) has local maxima corresponding
to stable steady states and a minimum at the X-value for the

unstable steady state.

The master equation in matrix form is

EB = Q 2! (4.7)
where
Axy = a+(y)6x_1’y + a_(y)6x+1,y + a(y)6xy. (4.8)
Its solution is (section 2.3)
w A.t
P(x,t) = P (x) + 2 d.P.(x)e j , (4.9)
5 .=1 3 j
3
where

P.(x)P(x,0)
6. = 2 3 , (4.10)
x Ps(x)

 

87
and Aj and Ej are the nonzero eigenvalues and the correspon-

ding right eigenvectors of A respectively. In particular, if

P(x,0) = 6Xy, then

on. = Pj(Y) / PS(Y)o (4.11)

4.2. TRANSITION BETWEEN THE TWO STABLE STEADY STATES

For numerical studies the following values of the rate

constants used by Gillespie36 are selected:

c1 = 3 x 10'7 (1mo1ecu16‘3 time-1),
c2 = l X 10.4 (molecule -3 time-1),
and c3 = 1.5 X 10.3 (moleculen1 time-1);
then c4 = 1.5 (molecule-ltime-l) corresponds to the critical

point. To illustrate the transition between the two stable

5

steady states, consider c = 3.33333 and B = 1.01 x 10 as

4

an example. Figure 4.1 shows the stationary distribution
and the eigenvectors corresponding to the three eigenvalues

closest to zero. As seen in section 2.3, A0 = 0 and all other

eigenvalues are negative. Let the eigenvalues and eigenvectors

be arranged such that A0 > l1 > 12 ... . In the example

illustrated in Figure 4.1 the first two nonzero eigenvalues

Figure 4.1

88

Eigenvectors of the transition matrix A
for the Schlégl model, when the parameters
correspond to a point in the interior of
the multiple steady-state region.

(c = 3.33333 and B = 1.OIXI05).

4
The eigenvalues are:

a) 10 = 0,

b) )1 = -3.8404x10'8'
c) )2 = —1.0537,

d) A = -l.5619.

85(1)qu2

89

3.0 -

2.7-

[8'

[5)-

Q9 '-

05 "

 

 

 

 

 

0.3 P

 

 

1 1 1 A
0 DO 200 300 400 500 600 700

Figure 4.1.a The stationary mode ES.

1 L
800 900

L
IOOO

90

I

|.OI

OBI '-

0.6| -

O.4| -

 

 

- O.|9

T

P, (x): to2

-O.39

I

'059

I

I

'079

 

-' '9 L l L I l l
O O [00 200 300 400 500 600

I

Figure 4.1. b The transition mode g

 

\/

1
700

1.

l
800

900

IOOO

P2(x)x104

H9-

099-

079-

059'

039

(M9

 

\‘fi I

'00!

I

-QZ|

~04“-

 

-Q$|

1

I

-0£I

 

Figure 4.1.c

91

1 1 1 1 1 l
0 ICC 200 300400 500.600

X

The eigenvector g

2

 

l

1 1 1
700 IKE) 900 I000

corresponding to 12.

P3001: IO2

92

I36 '-

|.46 )-

LOG

0.66 '-

0.26 j

‘014 1-

 

 

‘094 1-

4.34

l

 

4.74

l

 

 

1 1 1 1 1 1 1 1 1
0 I00 200 $0 400 500 600 700 800 900 IOOO
K

Figure 4.1.d The eigenvector 23 corresponding to 13.

93

1 = .3,84o4x10—8 and A2 = -1.0537) differ by several

orders of magnitude. This means that the eigenmode g

(A
1
persists for a long time after all other modes have decayed.
Therefore, after a sufficiently long time, P(x,t) has

significant contributions only from as and 2 Moreover,

1.

21 has a positive peak corresponding to one of the peaks in

the stationary distribution £8 and a negative peak correspon-

ding to the other (see Figures 4.1.a and 4.1.b). As 21

decays, P(x,t) decreases in one steady-state region and

grows in the other. Let 31 be called the transition mode.
Figures 4.2 - 4.5 show the time evolution of the

probability distribution for various initial values of x.

A quasi-stationary distribution is established on a short

time scale. Then this distribution slowly relaxes to the

stationary distribution in accordance with the Kramers rule

(explained in section 4.3). In general the quasi-steady

distribution is determined by the initial distribution.

If the initial distribution is a delta function peaked at

an x-value close to one of the stable steady states, then

the quasi-stationary distribution is peaked almost

exclusively at that steady state. This behavior was noted

by Oppenheim 33 3131. When the initial x-value is close

to the unstable steady state, however, the quasi-steady distri-

bution has contributions from both stable regions. (Figure 4.4)
The direction in which the transition occurs depends

on the sign of a in equation (4.9), which in turn depends

1

P(x,t)x 10"

3.6 -
3.2 .
2.8 - W
2.4 -

2.0 1-

 

 

0.4 )-

 

Figure 4.2

 

94

L I

1
800 900 IOOO

Time evolution of probability distribution
when the parameters correspond to a point
in the interior of the multiple steady-state
region. c4 = 3.333333 and B = 1.01Xl05 (same

as in Figure 4.1). Initial distribution is

6(x-500).

 

 

95

 

 

 

1-1‘
3.61- "'0
mo7
1-0'
3-2-
2.81-
24-
N
9 20-
1
E:
a Lab
L2-
0.8-
t = 108 7
O i L l l l J l l
o no 2a: Km (“D «X» ax) 7pc am: 9a: mm:
x
Figure 4.3 Time evolution of probability distribution.
1.OIXl05.

c4 = 3.33333 and B =

Initial distribution is 6(x-50).

P(x,t)x 102

3.6 -

32 ..

2.4 -

 

 

 

Figure 4.4

28- nr’Z//”

96

 

200 300 400 500 600 700 800 900 IOOO

Time evolution of probability distribution.
c4 = 3.33333 and B = 1.01x105.
Initial distribution is 6(x-253), peaked at

the unstable steady-state X-value.

97

3.6 "

3.2 v-

2.4 -

2.0 -

P(x.t) x )02

 

"IO
A.) A?

l L
0 DO 200 300 400 500 600 700 800 900 IOOO

 

 

Figure 4.5 Time evolution of probability distribution.

c4 = 3.33333 and B = 1.01X105.

Initial distribution is 6(x-230).

98

on the sign of P1(x0) if the initial distribution is a delta
function at x0. Further, the two regions between which the

transition occurs are separated by the zero of £1 and not
by the minimum of 23' These points differ significantly

in some cases as will be seen below.

4.3. THE VALIDITY OF KRAMERS RELAXATION

A bistable system is said to undergo Kramers relaxation54

if the probability distribution equilibrates within each
region rapidly and then diffuses slowly from one region to
another thereby effecting a pure transition. Gardiner55 has
reformulated the Kramers rule based on a one—dimensional
Fokker-Planck equation. Similar results can be derived from

a birth-and-death master equation describing a bistable

system. For this purpose, let us define

X
M1x,t) = 2 P(y.t), (4.13)
y=0
N1(t) = M(X2-l,t), (4.14)
N2(t) = P(X2,t), (4.15)
N3(t) - E P(x,t), (4-16)

x=X +1

99

2
n1 = E Ps(x), (4.17)
x=0
n2 = PS(X2), (4.18)
and n3 = Z Ps(x). (4.19)
x=X2+1

From the master equation (4.2) and from the detailed balance

relation (2.23) satisfied by Ps(x), it follows that

a _ P(x+l,t) _ P(x t)
at M(x,t) - a+(x)PS(x) [ 527;:TT— 527:7- ]. (4.20)

Dividing by a+(x)PS(x) and summing over x from X1 to Xz-l,

we obtain

X -1

é [a (>01> (xn’l —3 M(x t) -
_ + s at ' -

P(szt) - P(let)

(4.21)

100

If Kramers relaxation is followed, then we can set

P(x,t)

P(x,t)

With these,

To calculate

M(x,t)

Ill

equation (4.21)

1a+(x)Ps(x)1‘

N1(t) x

 

N1(t)

 

Ps(x) n

 

becomes

1 3

the left-hand side:

X Ps(y) =

n y=0

1

N1(t) [1 - ¢(X)]

 

<
for x X2,

 

 

for x > X2. (4.22)
= N2(t) - N1(t)
n2 n1

(4.23)

N1(t) Xz'l
n1 - Z Ps(y)

n y=x+1

1

for X < X , (4.24)

101

 

where
r -1 x§-1
n P (y) for x < X
1 y=x+l S 2'
@(x) = 1
-1 X
n3 -2 Ps(y) for x > X2. (4.25)
[ y—X2+1

Substituting (4.24) in (4.23), we get

N2(t) N1(t)

 

 

 

d - -
K 3EN1(t) — n n , (4.26)
2 1
where
X2-1
l - 4(x)

K = Z . (4.27)

x=X1 a+(x)Ps(x)

Similarly by summing equation (4.20) from X2+1 to X3, we can

obtain

N2(t) N3(t)

d

u EENz (4.28)

(t)

 

 

n2 n3

where

102

X
3 1 - 4(x)

“ = .3. +1 auxwxr

 

(4.29)

Equations (4.26) and (4.28) permit a three-state interpretation

similar to Eyring's transition state theory.

 

 

 

(Kn1)'1 (un21’1
. 4} 1*“) .
Xl-region , X2 , X3-region.
‘ —1 ‘ -1

The probability distribution for the Schlégl model
satisfies equation (4.22) if the following conditions are met:
1) The eigenvalue corresponding to the transition mode must be
sufficiently smaller in magnitude than all other nonzero
eigenvalues. A quantitative criterion is given below.

2) The zero of the transition mode must match the minimum of
the stationary mode. Moreover, Ps(x) and P1(x) should be
In other

proportional in 0 < x < X and similarly for x > X

2’ 2'
words, the shapescfifthe transition mode and the stationary
mode should match in each region separately.

A transition between stable steady states is well defined
if there exists a time T when a considerable fraction of the
transition mode remains after all faster modes have decayed
almost completely. It is possible to find a relation
between A and A such that the coefficient of g has

1 2 l

decayed to a fraction f1 of its initial value and the

103

coefficient of £2 to a fraction f2 of its initial value ,

i.e., e = fl,

and e = f . (4.30)

The required relation is

_2 > log f / log f (4.31)

1.

Thus the coefficient of the transition mode will remain at 90%
(50%) of its initial value after the coefficient for the
faster mode has decayed to 1% of its initial value, if the
ratio of the eigenvalues exceeds 43 (6.6 respectively).

This condition differs from that given by Oppenheim g; 218 in

terms of the difference [A -1 Their criterion allows the

2 1|-
possibility that the transition mode itself is negligible
compared to the stationary mode by the time the next faster
mode has decayed substantially.

The first three nonzero eigenvalues have been plotted
versus B in various regions of the parameter space (Figures
4.6 - 4.8). Decay of the transition mode is very slow well
inside the multiple steady-state region, but it becomes

faster as either the critical point or a marginal stability

point is approached. Figure 4.9 shows the region in which

-5.0

-60

 

Figure 4.6

 

104

 

The first three nonzero eigenvalues of the
transition matrix A. c4 = 3.33333.
The dashed lines enclose the B-range in which

there are three steady states.

105

0.0 [-

-066

4.32

l

“1.98

'2.64 -

“3.30

'3.96

'4.62

I

’528

“5.94

 

-650 1 1 1 1 1 1 1 1 1 1
5.0 5.35 5.7 6.05 6.4 6.7 7.l Z45 7.8 8)!) 85

Bx I0"

Figure 4.7 The first three nonzero eigenvalues of the
transition matrix Q. c4 = 1.7 .

The dashed lines enclose the B-range in which

there are three steady states.

106

«240

“300

'360

‘420

 

 

 

BxKT‘

Figure 4.8 The first three nonzero eigenvalues of the

transition matrix A. c4 = 1.3 (no multiple

steady states).

107

Figure 4.9 Regions of the parameter space where the
conditions on the eigenvalues (equation
(4.31)) is satisfied. The solid lines

consist of the marginal stability points.

f2 = 0.01
f1 = 0.9 for the dotted line.
f = 0.5 for the dashed line.

1

108

 

 

l0-

..-o. x m

2.5

20

|.5

|.0

109

condition (4.31) on the eigenvalues is satisfied.

Figures 4.10 and 4.11 show the stationary distribution,
the first three eigenvectors with nonzero eigenvalues, and
the time-dependent distributions near a marginal stability
point. It is interesting to note that the stationary
distribution does not have a significant peak near the higher
steady state, but the transition mode does. Consequently,
the intermediate probability distributions are peaked at both
steady states. However, the relaxation to the lower steady
state is fast and there is no quasi-stationary distribution.

Figures 4.12 - 4.14 illustrate situations near the
critical point. In this region the zero of the transition
mode differs considerably from the unstable steady state
X-value. In Figure 4.14.b, for example, a transition occurs
between two regions of approximately equal size, eventhough
the stationary distribution is asymmetric.

Although the deterministic analysis sharply divides
the parameter space into single and multiple steady state
regions, the stochastic analysis reveals that in the multiple
steady state region close to the boundary, the behavior is
similar to that in the single steady state region, i.e.,
transition between multiple steady states and equilibration
in the vicinity of a single steady state are similar near
the boundary. Thus there is no abrupt change in the
qualitative behavior at the boundary separating the single

and multiple steady-state regions.

110

Figure 4.10 Eigenvectors of the transition matrix A
near a marginal stability point.
(c4 = 3.33333 and B = 9.6 x 104).
The eigenvalues are:

a) 10 = 0,
b) 11 = -1.6732 x 10'2,
c) 12 = —7.3902 x 10’1,

d) A -1.4365

111

4.0 r-

3.6

3.2 -

2.8-

 

 

 

l J
0 DO 200 300 400 500 600 700 800 900 000

Figure 4.10.a The stationary mode gs.

P, (1)1110'

7.|9 '-

569 "

4.|9 -

2.69 .

 

l.l9

'L

 

‘03)

'l.8l

I

'3.3|

I

'4.8|

I

I

‘ 6.31

 

- L l
1&0 I00 200

Figure 4.10.b

300

The

112

 

transition mode g

1.

113

1.44 -

H4-

084 -

0.54 -

0. 24

 

 

- 0.06

P2 (X)X '07

 

I

-0.36

-0.66 "

~096

'|.26 '

 

4.56 L L ‘

 

Figure 4.10.c The eigenvector 22 corresponding to 12.

114

zyzr n
2.22 -
1728
1.22 -
0.72 -
0.22 [J
~O.28

-0.78

 

 

0
P3 (11):: IO

1

-1.28

 

-1 .78

I

-2.28

‘I

 

- 2.78 L

 

L l l L L i l I L
0 I00 200 300 400 500 600 700 800 900 I000
X

Figure 4.10.d The eigenvector g3 corresponding to 13.

115

 

 

 

 

L////,t = 1000
3.6 r
L—t = 100

3.2 -

2.8 '-

24-
"9 2.0 -
1
fifin

1.2 b m

t = l tflOO
0.3 .- fl/
“IOOO
0.4 '-
o J l l l l I I
0 I00 200 300 400, 500 600 700 800 900 1000

I

Figure 4.11 Time evolution of the probability distribution

near a marginal stability point. c = 3.33333

4
B = 9.6x104. Initial distribution is 6(x-400).

4.4 -
4.0 r
3.6 -
3.2 r-

2.8

2.4

I

Ps (11):: (03

0.8 *-

0.4 1-

 

0.0

Figure 4.12

116

l 1 l l

1 1 1 1 1
(00 ISO 200 250 300 350 400 450 500
X

Figure 4.12.a

Eigenvectors of the transition matrix near

the critical point. c4 = 1.7 and B = 6.39X104.
The eigenvalues are a) 0, b) -4.9051X10-2,

c) -2.7271x10’1, and d) -5.3060x10'1.

117

 

 

 

‘68 L L

 

l 1 l l l l
0 50 l00 I50 200 250 300 350 400 450 500

X

Figure 4.12.b

P2001: IO3

 

118

 

 

-2

-3

-4

-5 l _L L l l L l l 1 ¥
0 so 100 150 200 250 300 350 400 450 500

X

Figure 4.12.c

4.73

3.73

2.73

L73

0.73

'027

P3 (11)): IO3

’22?

“3.27

-4.27

4.27L

T

I

I

I

 

'5.27

I50

200

119

l
250

I

Figure 4.12.d

1 L
300 350

1
400

 

L

450

500

P(x,t)x IO3

 

 

120

 

 

0 50

Figure 4.13

l L L I I L .. L I
ICC 150 200 250 300 350 400 450 500
X

Figure 4.13.a

Time evolution of probability distribution
near the critical point. c4 = 1.7 and
B = 6.39X104. Initial distributions are

a) 6(x-120), b) 6(x-300), and c) 6(x—182).

P(x,t)): 103

 

  

[00

121

 

1
250
X

1 1
(50 200

Figure 4.13.b

 

500

P(x,t)x 103

 

122

1 1 1 _1_
200' 250 300 350
I

Figure 4.13.c

400

 

 

1 1
450 500

P(x,t)x Io2

 

50

123

 

 

/, .

1 1 1 L 1 J
'00 I50 250 300 350 400 450 500
x

Figure 4.14.a

Figure 4.14 Time evolution of probability distribution near

the critical point and near a marginal stability
point. c4 = 1.7 and B = 6.45X104. Initial
distributions are a) 6(x-300), b) 6(x-100), and

c) 6(x-190).

801.”): IO3

 

50

I
I00

124

  
   

1 L 1 L L 8.. 1 1
(50 200 250 300 350 400 450 500
X

Figure 4.14.b

 

'125

'50 200 250

X

Figure 4.14.c

300

 

350 400

I
450

L

500

126

In summary, the asymptotic solutions of the master
equation for the Schlbgl model have been constructed from
the eigenvectors of the transition matrix. In the cases
studied, deviations from the Kramers rule are attributable
not only to the similarity of the decay times of various
modes, but also to distortions in the shape of the slowest

mode.

CHAPTER V

HYSTERESIS IN TRANSITIONS BETWEEN MULTIPLE NONEQUILIBRIUM

STEADY STATES OF THE SCHLOGL MODEL

Many models have been proposed and studied with the
aim of understanding regulatory phenomena such as chemical
oscillations and pattern formation in chemically reacting

open systems far from equilibrium.LEV-5'56"6O

One interesting
phenomenon associated with the presence of multiple steady
states in reacting systems far from equilibrium is hysteresis
in transitions between the steady states. Hysteresis can be
viewed as a special kind of oscillation caused by an external
variation of a system paramter. I have simulated hysteresis
in the Schlbgl model system both deterministically and
stochastically. Hysteresis has been observed by Creel and
Ross53 in the dissociation of N204 when irradiated with an

argon-ion laser. The steady state patterns of the chemical

system are similar to those of the model system.

5.1 SCHLOGL MODEL AND ITS STEADY STATES

In the Schlbgl model the chemical reactions

128

occur in a homogeneous open system; the numbers of molecules
of A and B (denoted by A and B) are assumed to be externally
controllable and they are constrained to be equal. Thus the
number of molecules of X (denoted by X) is the only system

variable. When the reaction parameter 0 defined by

exceeds a critical value me = 9, there is a range of B values
for which the system has three possible steady states,

X < X < X3, of which x and X are stable and X is unstable.

1 2 1 3 2
Outside this range of B and below the critical point there is
only one stable steady state. The boundary separating these
two regions consists of the marginal stability points.

In the stochastic formulation X(t) is viewed as a random
variable taking integer values in {x:0<x<w}. The state of the
system at a specified time t is then given by a probability
vector g(t) whose x-th component is the probability that
X(t) = x. The dynamics of the system are described by the

time evolution of P(x,t) [ P(x,t) - Px(t) ] according to the

master equation for the Schlagl models’9

5%p(x,t) = a+(x-l,t)P(x-1,t) + a_(x+1,t)P(X+1,t)

" a(xlt)P(XIt)l (Sol)

129

where

c
a+(x,t) = J B(t)x(x-1) + c3B(t), (5.2)
2
c:2
a_(x,t) = —— x(x-1)(x-2) + c4x, (5.3)
6
and a(x,t) = a+(x,t) + a_(x,t). (5.4)

This equation allows for time dependence of B. The parameters

C1 are related to the rate constants ki by

c. = k. n.! V (5.5)

where 111 is the molecularity of the i-th step in X, mi is the
total molecularity of the i-th step (e.g., n = 2 and m = 3),

1 1
and V is the volume of the system.10

5.2 DETERMINISTIC SIMULATION OF HYSTERESIS
For numerical studies we select the following values of

the rate constants used by Gillespie:36

c1 = 3 x 10.7 (molecule —3 time-1),
c2 = 1 x 10-4 (molecule ‘3 time-l),
- -3 -1 . -1
c — 1.5 x 10 (molecule time ); (5.6)

130

then c4 = 1.5 (molecule-n1 time-1) corresponds to the critical
point. For c4 values greater than this, the steady state

curve is S-shaped (see Figure 5.1). If we first prepare the
system at a single steady state in the lower branch and then
slowly increase B, X follows close to the steady state value
until near a marginal stability point and then it jumps to the
vicinity of the other steady state branch. If B is slowly
decreased, a downward jump occurs near the other marginal
stability point. This generates a hysteresis loop, shown as

the curves marked with arrows in Figure 5.1. The curves are

the solutions of the equations of evolution

C C
Q = J- B x (x-1) + c B - —2— X(X-l) (x-2) — c x
dt 3 4
2 6
9'2 = +
dt "B O (5.7)

Hysteresis occurs to some extent below the critical point
also, i.e., where there are no multiple steady states at all.
This is shown by the curves with arrows in Figure 5.2.

The area enclosed by the hysteresis loop is a measure
of the extent of hysteresis. In Figure 5.3, the area obtained
from the deterministic treatment (equations (5.7)) is plotted
versus the rate of change of B for various values of the
critical parameter, c4. Thfisincreases with an increase in 8

indicating the dynamic contribution to the hysteresis effect.

500r
450-
400L
350T
300-
250-

200-

 

L
5.8 5.93

Figure 5.1

131

I L I J L L I L

6.06 6J9 632 6.45 6.58 6.?! 6.84 69?
ENG”

Results from deterministic (curves with

arrows) and stochastic (noisy curves)

 

7.!

simulations of hysteresis in the multistable

region (c4 = 1.7). The S-shaped curve is the

set of steady states.fi 8 50.

3|5

280

245

2|0

I75

I40

I05

70-

35’

 

0.0 ‘
4.60 4.72

Figure 5.2

132

 

L

4.84 496 508 5.2 5.32 5.44 5.56 5.68 5.8

I 1 I l L J

8 1: IO"

Results from deterministic (curves with
arrows) and stochastic (noisy curves)
simulations of hysteresis in the single

steady-state region (c4 = 1.3). P =‘fiL

133

C4: LB

|.6

|.5

AREA x (0'5
4>

|.3

 

1 1 l l l ‘

_ 1
0110210304 5060708090

B

Area inside the hysteresis loop versus the rate

0

Figure 5.3

of change of B. (Deterministic results)

134

In the figure, the points corresponding to B = 0 were obtained
by integrating the area enclosed between the two branches of
stable steady states. This is the area obtained when vertical
transitions (transitions with dX/dB -% 1m) occur exactly at the
marginal stability points. It is interesting to note that the
area of the hysteresis loop extrapolates to this finite value
as 8 tends to zero in the multiple steady state region, and

it extrapolates to zero in the single steady state region.

The finite limit indicates a static contribution to the
hysteresis effect from the multi-state region, but below the

critical point the hysteresis is purely dynamic.

5.3 STOCHASTIC SIMULATION

Gillespie has deviced an algorithm to simulate chemical
reaction systems stochastically10 and has extended it to
include time dependent transition rates20 (section 2.3).
This algorithm first randomly selects the time of the next
reaction and then draws another random number to determine
which of the various reaction steps should be carried out.
If a reaction step hasoccurred.at time t, bringing the system
to state x, then the probability P1(T;X,t)dT that the next

reaction step will occur in the time interval (t+r, t+r+dr)

is given by

t+T
P1(T;X,t) = a(x,t+1) exp — I a(x,s)ds . (5.8)
t

135

The conditional probability that the next reaction is a
forward (or backward) step given that the next reaction occurs

after an interval T, is

a (x,t+r)
P2(tIT;x,t) = . (5.9)
a(x,t+T)

H-

The simulation algorithm starts by selecting two random numbers

r1 and r2 with uniform probability from the unit interval.

The distribution function corresponding to P1(T;X,t) is

T
F1(T;X,t) = J P1(s;x,t)ds
0
T
= I a(x,t+s) exp [ A(t) - A(t+s) ] ds
0
(5.10)
where
t
A(t) = I a(x,s)ds (5.11)

After performing the integral we obtain

T
Fl(r;x,t) = 1 - exp - J a(x,t+s)ds . (5.12)
0

136

We need to select 1 such that

i.e., { a(x,t+s)ds
0

- 1n (l-rl). (5.13)

In equation (5.13) we can replace l-r by r since both are

1 1
statistically equivalent. For the Schlégl model a_(x,t) is

independent of B and a+(x,t) depends linearly on B. When B

is varied linearly with time,

a+(x,t)
a(x,t+s) = a(x,t) + ——-———— 85. (5.14)

BIt)

Therefore,

a (x,t) 2
a(x,t)“: + —i——— 81—. (5.15)

B(t) 2

a(x,t+s)ds

OB—fifl

Combining equations (5.13) and (5.15), we obtain a quadratic

equation for T:

a+(x,t)
B—-———-— T + a(x,t)r + 1n r1 = 0. (5.16)
2B(t)

Solving the quadraric equation and selectig the solution that

137

is continuous in r and gives a nonnegative value for T, we

1

have the following relation between r1 and T:

 

 

I - a(x,t) + [ a(x,t)2 - B%%T a+(x,t) 1n r1 ];5
a+(x,t)8/ B(t)
if 8 > 0 or — 1n r1 < - Bét)[%a+(x,t)-a_(x,t)]
T =
1 B(t) _
a_ x,t (Te-Wk“ 1“ ‘1
L if B < 0 and - ln r1 > - EéEL[!5a+(x,t)-a_(xpt)]

 

(5.17)

From equation (2.65) we also have a relation between r2 and m:

+ 1 if r < a+(x,t+1) / a(x,t+r)

2

- 1 if r

2 > a+(x,t+r) / a(x,t+T). (5.18)

The time variable is advanced by T and the particle number
for X is incremented or decremented by one depending on
whether m is positive or negative. This procedure is then
repeated for the new state. The iteration produces a random
path of the system in time. This approach yields results
equivalent to those obtained by the master equation. Two
hundred sample paths have been simulated and the averaged
paths are shown in Figures 5.1 and 5.2.

In the deterministic simulation, the transition cannot

138

occur before the marginal stability point is reached (Appendix
A) ; but as seen in Figure 5.1 (noisy curves) , fluctuations

induce the transition within the multiple steady state region.
The areas obtained from the stochastic simulations As and

the correspoding deterministic areas Ad are plotted versus

B and c4 in Figures 5.4 and 5.5 respectively. The error bars
indicate expected deviations from the mean; they were
calculated from the standard deviations using the central

limit theorem}7 In general the area AS obtained from stochastic
simulation is less than that from the deterministic integration
(Ad), though the qualitative behaviour is the same in both
cases. Another important difference between the two kinds of
simulation is that as the rate of change of B tends to zero,

Ad tends to a finite value in the multiple: steady state
region, but AS tends to zero.

When B is changed at an infinitesimal rate, the
deterministic analysis predicts that the forward and backward
transitions occur at two different marginal stability points
whereas the stochastic analysis predicts a unique transition
point on the average in the multiple. steady state region.
This point was noticed and emphasized by Turner?1 However
this result should not be interpreted as a complete absence
of hysteresis effects in the stochastic simulation; rather,
our results show that hysteresis in the SchlBgl model is a
dynamic effect and that the static contribution predicted
by a deterministic analysis disappears when fluctuations

are taken into account.

139

AREA 11 10"

 

1 L L I

L
0 20 40 60 80 I00 I20 I40 I60 I80 200

 

Figure 5.4 Comparison of deterministic (upper line) and
stochastic simulation (lower line) results.

c4 = 1.7

140

AREAx )0"

 

Figure 5.5 Comparision of deterministic (upper line) and

stochastic simulation (lower line) results.

8 = 50.

CHAPTER VI
DETERMINISTIC STUDY OF MULTIPLE STEADY STATES IN COUPLED

FLOW TANK REACTORS

The question of relative stability of different steady
states arises when two or more of the steady states compatible
with fixed fluxes, fixed elemental composition of the system
and external conditions are asymptotically stable to small
concentration fluctuations. Resolutions have been proposed

on the basis of stochastic analysisBG, thermodynamic

30-32, and mixing experiments29 (section 3.4).

analysis

The purpose of this chapter is to describe the new
steady states and dynamical behavior of coupled systems when
each exhibits multiple steady states. It is also shown how
the outcome of mixing experiments with exchange of material
between two continuously stirred flow tank reactors can be
predicted from a deterministic analysis; the results obtained
differ from those predicted by Noyeszg, except for simple
cubic reaction rate laws. It is shown that deterministic
studies of coupled tanks do not yield information about the
relative stabilities of single tank steady states, even when
the outcome of a mixing experiment is independent of the

mixing method, i.e., independent of the changes made in the

coupled rate constant as a function of time.

141

142

6.1 MULTIPLE STEADY STATES IN A SINGLE CONTINUOUSLY STIRRED

FLOW-TANK REACTOR.

A model continuously stirred flow tank reactor (CSTR)
is characterized by two assumed operating conditions:

1. The volume of material in the tank reactor is fixed
and there is a negligible change in the density of material as
a consequence of reaction.

2. Mixing is assumed to be sufficiently rapid that the
material is uniform in composition throughout the tank.

The concentration of reactant i in the input is Rg,and
if product j is present in the input stream, its concentration
is P3. The reactant concentration in the output stream is
identical to that in the tank. The flow of material through

the tank is characterized by a flow rate constant k the

ol
inverse of the average residence time for material in the

tank; equivalently k is the ratio of the volume of material

0

entering the reactor per unit time to the total reactor volume.
In the following development we restrict attention to

the case where there is a single reactant R and a single

product P in the tank, as in the work by Noyeszg. With a single

and if the tank is

reactant and product, if R = R + P

T 0 0'
flushed with material from the input stream so that the total

reaction plus product concentration is equal to R initially,

T
then deterministically R + P = RT at all times.
In the tank, reactant R is converted into product P

by reaction according to the rate law

143

dP/dt = - dR/dt = V(R,P) = v(R,RT), (6.1)

where the last equality is valid because P can be replaced by
RT-R. With the addition of the rate of change of R due to flow
through the tank, the total rate of change of reactant

concentration is

dR/dt = kO(Ro-R) - v(R,RT). (6.2)

Three different polynomial forms for v(R,RT) and the
deterministic steady states associated with these models are

)

used as explicit examples in this chapter. The forms for v(R,R

T

are

1. Noyes' cubic model. For this model

V(R,P) = a(R-P/K) (1 - bR + ch), (6.3)
and the variables have values a = 8.97 x 10"4 5.1, ab =
5.985 x 10'2M'ls’1, ac = 1.00 x 10’2M'2s'1, K = 106 and
R + P = RT = 0.1M

2. Noyes' quartic model. For this model,

V(R,P) = a' (R-P/K) (1 - b'R + c'R3), (6.4)
with a' = 10’35'1, b' = 50 M'l, c' = 1.87 x 104 M'3, K = 106
and R = R = 0.1 M.

144

3. Two quintic forms Q1 and Q2, for which

v = (R-P/K) ( a + bP + cp2 + dP3 + ep4), (6.5)

K = 106 and RT = 0.2 M. The quintic expressions are aphysical
as reaction rates, but they are useful in illustrating the
difference between mathematical results valid without
restrictions on v(R,RT) and those steady-state and dynamical
results applicable when v(R,RT) is any cubic polynomial. For

models Qland Q2, the constant sets take the following values:

01:
a = 3.028222 x 10"9 ’1
b = 3.822292 x10"3 M'1 s‘1
c = 1.984543 x 10‘2 M'2 s’1
d = 6.649973 x 10'2 M'3 s"1
e = 0.999999 M'4 s'l,

Q2:
a = 2.966667 x 10"8 '1
b = 3.788482 x 10‘3 M'1 s'1
c = 1.969993 x 10"2 M'2 s-1
d = 6.499974 x 10'2 M'3 s'1

4 1

e = 0.999999 M-

Corresponding to a given value of RO

(less than or equal to RT),

145

there is a range of k0 values for which there are three real

positive solutions of dR/dt = 0 for Noyes' cubic model and
for the two quintics. The three steady state solutions are

designated Ra' RB' and Ry, with Ra < RY < R < R In Figure

B T'
6.1, the steady-state R concentrations are plotted as

functions of k0 for six different R0 values in the range from

0.03M to 0.1M, with RT = 0.1M for Noyes' cubic model (in Noyes'
work RT 5 R0). For the quintic models, attention is restricted

0’ k0 = 7.659631 x 10"4 s'1 for model 01
and k0 = 7.5924 x 10"4 s-1 for model 02' In general, for the

quartic model, there are four steady state solutions of the

to single values of k

deterministic kinetic equations, but the branch of solutions
with the smallest R values is unstable to small concentration
fluctuations; Noyes has plotted the remaining three solution
branches as functions of k0 in reference 29.

The stability of the steady-state solutions of equation
(6.2) may be determined by a linear stability analysis. If R
is a steady state solution of the deterministic kinetic

equation for reactant concentration, then the time-evolution

of a small concentration fluctuation 6R about R is governed by

I 6R, (6.6)
J

 

on the unidimensional subset of the full concentration space

defined by the condition R + P = R Therefore a single tank

T0

146

 

 

qu4' 1 I l I I I I I I’ I I I r‘# I I
- .1
0,40 8- R -0.100 -
0 0.075 -
’ 0.050 .
Q36 " 0.040 .
0.32,

 

 

 

0281

024*

KDR ‘
020'

OE?
02?
008-

004

        
     

  

, I

 

 

00-1.
00 02 CM-

I I I L

K) 12 I4 I5

06 08
4
101(0

Figure 6.1 steady states for the cubic model in a single
flow tank reactor. The reactant concentration
(R) in the tank at steady state is plotted as
a function of the flow rate constant (k0) for
various values of reactant concentration (R0)

in the input stream.

147

state is stable if k0 + (3v/8R)E > 0, marginally stable if

k0 + (av/3R)§ = 0, and unstable otherwise. For all of the
models (6.3 - 6.5), linear stability analysis shows that the
intermediate concentration y branch is unstable to small
concentration fluctuations, while the a and 8 branches are
stable. The points of coalescence of the y branch with the
a or B branch are marginal stability points (see Figure 6.1).
Our aim is to investigate what happens when two tanks,
one initially in the a state and the other in the 8 state,
are coupled by exchange of material. As part of this
investigation we find the possible steady states of the
coupled tank system and their dependence upon the coupling
rate constant kx. This analysis shows that the outcome of

mixing experiments is not an indicator of the relative

stability of the asymptotically stable steady states.

6.2 COUPLED FLOW-TANK REACTORS

Let us consider two tanks, operated with the same flow
0’ T’ and R0 = RT (i.e., each has
an input stream containing reactant only). Conversion of

rate constant k the same R
reactant to product is assumed to proceed according to the
same rate law in each tank. If the two tanks are coupled by a
machanical pumping of material from tank 1 to tank 2, while
at the same time, material is pumped from tank 2 into tank 1
with the same flow rate constant kx' then the rate of change

of reactant concentrations in the two tanks are

148.

)

de/dt k0 (RO - R1) - V(R1,RO) + kx'(R2 - R1

kO (R0 - R2)

and dRz/dt - V(R2,R0) + kx (R1 - R

2).

(6.7)
Thus, independent of the values of kx' the points (R1,R2) =
(Ra'Ra)' (RB,RB),and (Ry’RY) are always steady-state solutions
of the coupled tank equations.

Because the condition R + P = RT reduces the single-tank
problem to a single—variable problem, there is a potential
that governs the deterministic time-evolution of the single-

tank system, and there is also a potential in the coupled-

tank case. If V(R) is defined by the condition

R
V(R) = I V(X,Ro)dx, (6.8)

and if 0 is defined by

2
R + V(Rl)

¢(R1,R2;k0,kx) = - k R R + 1/2 k0 1

0 0 1

2
- koROR2 + 1/2 kOR2 + V(Rz)

2
+ 1/2 kX(R2-R1) , (6.9)

then 0 satisfies the conditions

deldt = - 30/8R1 (6.10)

149

and
dRZ/dt = - 80/3R2. (6.11)

For fixed kO and kx' then

and d0/dt = 0 only at steady states.

Figure 6.2 shows the possible steady-state concentrations
in tank 1 as a function of kx for Noyes' cubic model with
k0 = 1.5 x 10-55-1. For kx = 0 and for small kx' there are
nine possible steady-state pairs (R1,R2). These are designated
0a, dY, dB, YB' yy, Yd' 80' BY and BB in order of increasing
reactant concentration in the tank. The subscript indicates
the corresponding steady-state branch in tank 2; i.e., if
both de/dt = 0 and dRZ/dt = 0, and if tank 1 is in state 08,
then tank 2 is necessarily in state 80' The steady-state R

concentration in tank 2 may be computed from the steady-state

concentration in tank 1 since

2 1 0 - R1) ] / kx' (6.13)

0

The number of steady-state solutions of the coupled-tank

equations obviously depends upon the value of kx' It decreases

10

0.44

0.40

0.36

032

0.28

0-24

0.20

0.I6
0.I2

0.08

0.04
0.00

Figure 6.2

150

 

 

 

 

I I I I I I T I I l l l I I I I
- 4
q
—(
:4
.I
)— .11
h d
- d
)- .1
F -—1
I- -I

 

 

 

1114111111111

00 2.0 4.0 6.0 8.0 l0.0 l2.0 I40 I60
5 .
l0 kx

 

Steady states of coupled reactors system for
Noyes cubic model. The reactant concentration
in one of the tanks (R1) is plotted versus
the mixing rate constant(kx) for

5 -1

k0 = 1.5x10' 5

151
from nine to five and then to three as kx is increased.
Figures 6.3 and 6.4 show the possible steady-state concentra-

tions in tank 1 for Noyes model, as in Figure 6.2, but for

different values of the flow rate constant k0; k0 = 2.7 x 10"5
s"1 for Figure 6.3, and k0 = 2.516963 x 10.5 s.1 for Figure
6.4.

For the quintic case with the constant set Q1, the
behaviour of the stationary solutions of equations (6.7) is
qualitatively similar to that for Noyes cubic model. For small
kx’ nine steady states are found; with an increase in kx'
these coalesce pairwise to leave five states. With a further
increase in kx, only the three steady states for which R1 = R2
are found. With the constant set 02, in contrast, new steady
state solutions of equations (6.7) are found. Again, for small
kx' there are nine steady states; when kx is increased, new
steady state branches appear. In particular, over the kx
range from 3.688463 x 10'9 3’1 to 4.298794 x 10’9 5'1, there
are seventeen possible steady states for equations (6.7)
where v(R,Ro) has the quintic form Q2. With further increase
in kx' this number decreases to thirteen, then nine, five,
and finally three (all three with R1 = R2).

Linear stability analysis shows that the stability
properties of coupled tank differ from those for single
tanks. If (R1,R2) designates a steady state of the coupled

tank system, and if we consider a small concentration

fluctuation in each tank (6R1,6R2), subject to the

152

 

 

 

 

 

 

 

044 I I I I I I I | I I I 1 I I I
(14(37 4
0'36\ —
CX3EL- -
0.28 '— .1
)0 R1 .. J
C124I- -
020- _
(1H5’ 1
CIHZ- -
0.08 ..
r- -I
(Dffi4
000 J 1 1 1 1 l 1 1 1 1 I I 1 I 1 d
0.0 2.0 4.0 6.0 8.0 l0.0 I20 I40
105 k,

Figure 6.3 Steady states of coupled reactor system for

the cubic model. k0 = 2.7x10'5 s’l.

153

 

0.44IITITTITIIIII_FI

 

 

 

 

 

IDKKD' -
CM3ES- -
CX322- -
0.281- V -
KDF31 - _
024- -
ClZI)
CIKS- -
(IEB- -
(ICHBF- -
r- -
(DCM4. _
0.00 1 I L 1 1 l 1 1 1 .

 

00 2.0 4.0 6.0 8.0 10.0 12.0 (40'

105k,

Figure 6.4 Steady states of coupled reactor system for

the cubic model. k0 = 2.516963x10'5 s’l.

154

conditions R1 + P1 = RT and R2 + P2 = RT,

evolution of the fluctuation may be approximated by the

then the time-

linear equations

 

 

 

 

 

6R 6R
6%? 1 = “1.4. I] (6.14)
6R2 6R2
where
r 3v ‘
k0 + 8R - + kx - kx
R
1
g = (6.15)
3v
kx k0 + 8R - + kx
I R2 .

The stability of the steady state (R1,R2) to small fluctuations

depends upon the eigenvalues of the matrix Q, A and A If

1 2'

both A and 12 are positive, the state (R1,R2) is stable.

1
Marginal stability occures when det g = 0. At this point,
at least one of the eigenvalues of g vanishes, and in the
linear regime, a fluctuation of the character of the associated
eigenvector is neither damped nor amplified in time. The

condition for marginal stability of the (a '80) state is

B

 

155

Equivalently,

8v
(k0 + “R

 

 

(6.17)
If kx is less than the kX value at marginal stability, then the
inhomogeneous steady state (aB'Ba) is stable to small concent-
ration fluctuations.

Each steady state of the coupled-tank system is an
extremum point of 0. At a stable steady state (R1,R2) 0 has a
local minimum; for a state unstable to any fluctuation, 0 has
a local maximum, and for a state stable to certain perturbations
(6R1,6R2) but not to others, 0 has a saddle point. Contour
plots for the potential 0 for Noyes cubic and quartic models
are shown in Figures 6.5 - 6.10. In Figure 6.5 - 6.7, 0(R1,R2)
is plotted for Noyes cubic model with k0 = 2.517 X 10'5 s’l;
the change in potential with change in kx is illustrated by
the sequence of Figure 6.5 with kX = 0 (no coupling between

the tanks), Figure 6.6 with kx = 4.0 X 10.5 3.1, and Figure

6.7 with kx = 9.061 X 10-5 5.1. The kx value for Figure 6.7

is the point of marginal stability for the coupled-tank state
(08,8a). Figure 6.8 also shows results for the cubic model,
but in this case, k0 = 1.5 x 10’5 s"1 and kx = 3.9 x 10'5 s'l,
the marginal stability point of the inhomogeneous coupled—
tank steady state for this k0 value. Figures 6.9 and 6.10

show results for the quartic model obviously lacking the

100 R2

156

 

\\

 

 

Figure 6.5

      
 

100 R1

The potential function 4

with k0 = 2.517><1o'5 s'1

The extremum points of ¢

steady states (R1,R2) of

system.

for the cubic model

 

0

and k = 0.
x
correspond to the

the coupled tank

157 '

 

100 R2

 

 

 

 

 

Figure 6.6 The potential function 0 for the cubic model with

k0 = 2,517x10—Ss-1 and kX = 4.0X10-5s-1.

158

 

 

 

 

Figure 6.7 The potential 0 for the cubic model.

k0 = 2.517x10"5 s"1 and kx = 9.06lx10'5 s’l.

100 R

159

 

 

7”;- “- “ ‘z‘- “~ _____,—-¢-' 6' f ’[Ifl’:_’; —-‘\v ‘
4 ° 0 Liar—‘- m “ #755561», fiz'z/{M‘RQ‘X‘I

’.f’;fF-—-‘-“‘ “ ‘ c.._ v...- ‘. .4—

1’ .4}... s _ // /,I"//7l”/' “‘.\\‘\\J

 

 

‘\

,R ‘. , .‘ . ,

m. -—--- x .68:
\ . \7 Is \ '- “'\ V \ ‘9 ' .

Vfl‘x‘w x)»: ‘ \x. “‘ ..- \ ...—-—\\ -. "II
fl$§\\\\x‘\\-i .:‘\\ ‘ \\"\\\}\\\\\ \ WQ?‘ '\\\‘\\‘\‘.‘ .1
1 . 0 4-«2 ‘\x\“ ~\\‘“\..\.\\‘ 1 \‘t x\‘ \. m "
s. \ ‘ \ \ o \ \ ~\ \ ' ' ' \I
fl‘\\{‘\_\\\'{x{.‘\$\\ \\\‘\\\‘{\:’\\ :'\\\ \ fl ‘ '\ \‘1 14
PN‘.‘\\\“\\‘\.\. \. {k \‘\\\"\\\\¥;\‘\\\\'\\\ 5“ \ \\ \\\. . 3.1

.- “ - \l -_ .‘ ‘ I I x ‘5." v \ '. .1. "
‘\ \‘\.\‘ ‘\.\{\:\\\\1\. ‘\.\:\‘1\\\\‘\'\ \\ '\.\‘ x \ m‘\\\\ I, I '1‘" I

'4

fl ‘\ 1 1'1 “‘0“ \xx\'.\ .\.
I“ \ . . \ \ . \
0 ° 5 Ir’”\“3\\\\\\\¢§\\\\$\NW“- \\ ‘ \ \\

 

.' \\\\ \\I , -1,"1"'."1')' ’1'). I. '

) 1““vI'I“ )(IIIIIIIII 1 -2 , ..,.

a I I] > I I II I} . 1' ,t) ' I I 3 3, [I I l ./ . if; ,1" ,1" XIII/.2341.

0 . 0 I: -1331?! Ill/'IIII’IX‘I’ILQIIIII/ ’1 ”(VIII/16¢! {Ii/r’fa’”. {/ij . 'j'
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.

.'a
.1

 

100 R1

Figure 6.8 The potential 0 for the cubic model.

k = 1.5x10‘5 5‘1 and kx = 3.9x10'5 s

-1
0 o

160

    
 

   
  

     
 

 

-’%#g—Ifl :7, "/7 0;?”092’5-‘5—‘04.
-' _— "-o' x -.-z / x «J;- s. .
«ref/“‘3’; véfi/zwrvfoz/A-«Rm
f§/"p"_L-/L/ 1,11 ,‘f. . I/:/' I./ fix," A \ 2‘31

, - - . . a My

I ?;//,/l/;//l’”:” ,' / /’/l
(4717/ 511/741], / ,1" /
‘Iy (”f/(0,1,; / // ‘/ ,,
3.0 {IIAVHO/‘/ ¢._~H
1‘. . l
I

I .' - x - ,
I 1 ’1 .. \‘~ . \ wow
V/fl/ F1 .4 1' // / \ \\ \\\\\ \\\ /
,‘ If} I _ ‘\ .. \ '
I 1' "

 
   

  
 

 

II!

 

 

 

 

 

Figure 6.9 The potential 0 for the quartic model.

k0 = 3.374xlo's s‘1 and kx = 1.28x10'

5-1

In

161

.u .lla.) w'flylxtllfl11y
qllfl‘QAijllulrfl'M‘ulili ..J . u ..1 . . .

y ‘ (I it t .. /.‘..,/ I , .. . ,
gnu/...flwfljflhmué . .‘ . .: . . .;. ..m
a} l/lH.//J.///... l/Iurlrhoy. .. . . ‘ . . .. . ..II ..7
Ir/r //./..U/./ 21”: . , :yf
/. .0. x/ x //J / / ..f.. . r
. , , /,, :. ”z;
a // . 12...?
.h _ ./ ../ .40
x x . .
.I(l\ \ \ . x / 1,

(I-.. \\\.\\ \\\ \\\\.\\ \ a,
l‘ ‘\ \ \ .
(\ ~.\ \\ \ \\ \\ \ \ a
/ \. \ xx x t .\ 1.11;],
11...! \\ x \\ \ \ll / 4 , ,1/
//i|\. \ \ 1/ .r 2.5..
. \ \ . '1’.)

                  

   

'1
f}

ff};
, , ,

_;—

“—L—r;
‘3’—

:11
I
I
I
l
\

   

   

 

__ __ I, {In
. I .
_ .F... .. / . . \ .111
, x z . . .\ .r
I I, X / \. t .r/
- V
... ,//// rr .1 34/91}!!! \ x. t. x r.
0.”; oz /. i ..l. x. ,/
.5. z , I, ..r. f; ,1
:1. z , (r I-.!| , I ,4
,_x..,..u.,/ 1 // ...!1. 115:1! ._ J.
In“ .I.I/.l .lr. IlJlIf / J
5””? .r//. I'lf/IJHslulll 1|»!fIJ 1! . _
a I .I.’ : I ll! IIIIIII| ’1' ll . i y —
if; /:H;/ I. I}! .a a; x _ . \
II.- n 4. "Illrl J) .
,vwr , ’1‘. .../huh. V: .11.] IfoHIIf 1|...wa I: . . ..
ufifi%flV%WWHWWHHHhU .., . . .bx,:.1 .6
«27,11, (Ii/Ell: ...." I‘ll! [I’ll II. I‘ . . .. . l/II . . t .. ‘ x
.1: .../4/4..JMWW1..{MJU.I.IINH1{W-AHI . ‘ .. ..h t . H h .. .(r . . I. - x \L
_ J # .
o 5 o 5 5
Q o i O I O I I O
4 3 3 2 N 2 1 1 0 0

-4 -1
S .

1XIO

2.5
and k = l.
x

0

2.
100 R1

= 3.5><1o'5 5'1

1:0
ko

0.5

 

Figure 6.10 The potential ¢ for the cubic model.

010

162

symmetry present for the cubic model. Figure 6.9 shows the

5 -1 4 -1

potential for k = 3.374 x 10' s and kx = 1.28 x 10’ s ,

0
again the marginal stability point of the inhomogeneous
coupled-tank state; Figure 6.10 shows results for k0 =

3.5 x 10"5 s'1 and kx = 1.1 x 10'4 5'1.

Next the effect of changes in the coupling constant kx
are considered. First we consider gradual increases in the
coupling constant, beginning with the tanks uncoupled (kx=0)
and in different steady states. If kX is increased
sufficiently slowly, the coupled tanks will always be found
at time t in the inhomogeneous steady state appropriate for
kx(t)' until the marginal stability point is reached. The
equation of motion for the system point in the (R1,R2) plane

is derived as follows: Let (§1{§2) be a steady state solution

of the coupled-tank equations

0 = R (R0 - R

O - V(R1,RO) + kX (R - R )

1’ 2 1

0 = k (RO - R

O - v(R2,R0) + kx (R1 — R2), (6.18)

2)

and consider the changes in R1 and §2 produced by an

infinitesimal change in kx to kX + ka. The new steady state

(R1+6R1, R2+6R2) satisfies
0 = k0 (R0 - R1 - 6R1) - v(R1 + 6R1, R0) + (kX + 6kx)
X (R2 + 5R2 - Rl - 6R1)

163

0 = k (R0 - R

0 - 6R2) - V(R2+6R2, R0)

2

+ (kx + kaHR1 + 6R1 - R2 - 6R2).
(6.19)
Away from marginal stability points of the coupled-tank
system, the changes in steady-state concentrations in the

tanks can be expanded as series in powers of 6kx. To first

order in ékx, 5R1 and 6R2 satisfy

_ _ - 22 _ - _ —
o - kOGRl 3R!§ 5R1 + kx(6R2 6R1) + kx(R2 R1)
1
o - — k 6R - fl! 6R + k (6R - 6R ) + k (E — fi )-
' o 2 an E 2 x 1 2 x 1 2 '
2
(6.20)
i.e.,
(SR "13 51$
- g 1 = 1 _2 akx, (6.21)
6R2 Rz-Rl

 

 

 

 

where the same matrix g determines the damping of concentra-
tion fluctuations. The equation has a solution for 6R1 and
SR

2 linear in 5kx, provided that det g ¢ 0.

If kX is increased sufficiently slowly that the coupled-
tank system initially in the state (R2,‘Rg) satisfies the
steady-state conditions for coupled tanks with coupling

constant kx(t) at all times t, then the values of'fi1 and £2

164

satisfy the differential equation

(6.22)

 

 

If kX = 0 at t = 0, the initial conditions are R? = Ra’ R8'
0

orliyand R = Ra’ R or Ry. The matrix 5 depends upon k

2 k

B' 0' x,

R1, and R2 and it is singular when the state (R1, R2) is

marginally stable.

In addition to the case of adiabatic mixing treated
above, it is also interesting to consider the limit of
instantaneous mixing, with kx approaching infinity. In this

limit, only solutions with R1 = R2 can be found for the

coupled-tank steady-state equations. When the single—tanks

exhibit three steady states Ra < RY < R with Ru and R

8' B
stable to small concentration fluctuations and RY unstable,

tanks with reactant concentrations R1 and R2 prior to coupling

lie in the domain of attraction of the (Ra’Ra) state if

R1 + R2 8,

state if R1 + R2 > 2 Ry. This result applies to all systems

with a single reactant and a single product independent of

< 2 Ry, and in the domain of attraction of (R R8)

the form of V(R,P), provided that the system shows the steady-
state pattern specified above.

Next, we consider the position of separatrices between
the domains of attraction of different stable steady states,
and show how these move as the value of kx is changed. When

kX = 0, there are four stable steady states for the coupled

165

tanks, (Ra’Ra)’ (Ra'RB)' (RB,Ra),and (RB,RB). The domains of
attraction of these steady states are separated by the
perpendicular lines R1 = RY and R2 = RY. As kX is increased,
these separatices shift. If no new inhomogeneous steady states
appear with increasing kx' these separatrices coalesce when
the (aB'Ba) state reaches marginal stability; the state (a8,8a)
lies on the degenerate separatrix at that point. In the limit
as kx goes to infinity, there is a single separatrix for the
domains of attraction of the (Ra'Ra) and (RB'RB) states, and
as discussed above, it has the equation R1 + R2 = 2 Ry.

The cubic case shows several special features. If
v(R,RT) is a cubic with multiple steady states in a single

flow tank, then for some flow rate constant k the three

0!

- R = R - R . For this value of
B Y Y a

the potential ¢(R1 R k
I

steady states satisfy R

k kx) is symmetric with respect

2; o'
= 2 RY for any value of kx (see Appendix

0!
to the line R1 + R2

F). If R+ and R_ denote the solutions of the quadratic

equation
k + £2 = 0 (6.23)
0 8R R ' .
t
then R+ + R_ = 2 Ry, (R+,R_) is a possible inhomogeneous

steady state solution of the coupled tank equations (see

Appendix F). If c is defined by the relation

k0(Ro-R) - v(R,RO) = C (R - Ra)(R - Ry)(R - RB)

(6.24)

166

for this selected ko value, then there is a single separatrix
between the domains of attraction of the steady states (Ra’Ra)

and (RB'RB) when
- R ) . (6.25)

Because of the symmetry of the potential, in this particular
case, the separatrix remains fixed as kx is increased. Hence
the point (R+,R_) lies on the separatrix for every value of
kX greater than or equal to - %(RB-RY)2’

For the cubic mechanism with any other choice of k0,
the marginal stability point of the inhomogeneous state lies
on the separatrix when kx takes on the marginal stability
value, but with an increase in kx' the separatrix sweeps
past this point.

To illustrate the contrast between the cubic case and

M8

the general case, let kx denote the kx value at which the

inhomogeneous (a,8) state is marginally stable, and designate
the R concentrations at marginal stability by‘E and E. The

values of E and E obviously depend upon k . For a particular

0
* __
k0 value, denoted k0, (a,8) lies on the separatrix of the

domains of attraction of the (Ra’Rd) and (RB'RB) states, and

it remains on the separatrix when the coupling constant is

MS

increased to kx + dkx. For anotherk. denoted kg.(3.§) lies

0!
on the separatrix of the domains of attraction of the (RQ,Ra)

and (RB'RB) states in the limit as kX goes to infinity.

167

* .
In the cubic case k0 = kg,

not hold. For the quartic model and the first quintic form,

but in general this equality does

it is possible to find marginal stability points (3.6) that

lie in the domain of attraction of the (Ra’Ra) state if kx

S

is increased suddenly to k2 + Akx, but that lie in the

8

domain of attraction of (R RB) if kx is increased to k: +

at
2Akx. Further, if new inhomogeneous steady states appear

with an increase in kx' as for the second quintic form, the
state (a8’8a) still lies on a degenerate separatrix at its
marginal stability point, but this is a separatrix between a

, ’b
homogeneous state and a new inhomogeneous (a 6d) state,

BI
rather than the (Ra'Ra) - (RB’RB) separatrix.

The analysis above essentially explains the outcome
of all possible mixing experiments conducted by changing
the value of kx in time. For example, equation (6.22)
describes the behaviour of coupled systems that are mixed

"adiabatically" , so that the system moves along a path of

steady-state points in the R R2 plane. When kx reaches kgs,

ll

*
the matrix Q is singular. Then unless k0 = k0,
mal increase in kx leaves the system in the domain of

an infinitesi-

attraction of a new steady state, which may be homogeneous
or inhomogeneous.

In the limit as kx approaches infinity, the results of
instantaneous mixing experiments are obtained. Any difference
in outcomes of mixing experiments with different functions
kx(t) obviously result from the kX dependence of the potential

(see Figures 6.5 - 6.10). The cubic case is unusual, because

-168

the outcome of instantaneous mixing experiments and gradual

mixing experiments is necessarily the same.

6.3. COMPARISION WITH THE NOYES ANALYSIS
In this section, the above results for coupled-tank
experiments are compared with the results of Noyes. These
results are also contrasted with the proposed relative
stability criterion. Noyes has proposed that dynamic equi-
stability at constant R0 can be determined experimentally
by mixing experiments, beginning with one tank in the Ra
state and the other (initially uncoupled) tank in the R

8
state. Then the reactors will go to the (a',B') state.

Initially k0 + %% > 0 for both the a' and 8' states. As kX
increases it decreases for a' and B'. If k0 + %% = 0 for a
a!
given kx value while k0 + %% ¢ 0, then according to Noyes
BI

a further increase in kX will drive the combined system to
settle in the 8 state and therefore the 8 state is more stable

than the a state. If k + 3v reaches zero first then the

o 372'
8'

a state is said to be more stable. As shown above and as noted

by Noyes, the actual condition for marginal stability of the

coupled tanks is given in equation (6.17). Therefore, if

k + £2 = 0, but k + 23 > 0, the coupled tank state

is stable. This raises the possibility of finding counter-
examples to Noyes prediction of the outcomes of mixing
experiments. Two mathematical counterexamples are considered

below. The first is provided by the quintic form Ql where

169

k0 + 3% = 0 when kX = 9.373 X 10-10 s-1 but an increase in
GI

kx by a factor of N40 is required to reach the marginal
stability point and when kx is increased slightly beyond kgs,
the system initially in the (a,B) state settles into the
(Ra'Ra) state, contrary to Noyes prediction. In the second

counterexample provided by Q2, a small increase in kx beyond

. . 8v
the pOint at which k0 + 5R a.

inhomogeneous tank state to the marginal stability point;

= 0 does suffice to drive the

but it does not settle into a homogeneous steady state.
Instead, it is driven to a new inhomogeneous state. When with
further increase in kx the new inhomogeneous state goes
unstable, the system again settles into the (Ra'Ra) state.
It is important to recognize that coupling between the
tanks alters the nature of the steady states in the single
tanks to the extent that information about the relative
stabilities of single-tank steady states is not available
from mixing experiments, even if all mixing experiments have
the same outcome. This point can be illustrated in several
ways. For example, the equation for the change in R concent-

ration in tank 1 when coupled to tank 2 can be written in

 

the form
dR k R + k R
-—1 = (k0 + kX) 0 O X 2 - R1 " V(RIIRT).
dt k0 + kX

(6.26)

170

This is identical to a single-tank equation with k0 replaced
by k0 (see equation (6.28)), R replaced by R0, and RT left
unchanged. Similarly

dR2 _ _.

g—— = k0(R0 - R2) - v(R2,RT), (6.27)

t

with

k0 = k0 + kX

RO - (kORo + kaZ) / (k0 + kx)
and R0 = (koRo + kXRl) / (k0 + kx). (6.28)

- - o n o o - _'
R0 and R0 satisfy the identities R0 < R0 and R0 < R0. At the

steady state (R1,R2) of the coupled tank system, the concent-

ration R in tank 1 satisfies the equation

R ) = O. (6.29)

Hence a single-tank system with parameters k R0, and RT has

In this

0'
a steady state with the same R concentration R1.
sense, the coupled tank steady states may be mapped onto

single-tank steady states. However, the single-tank parameter

R0 differs for the two coupled tanks, and both ED and R0

differ from the values for the original single-tank

171

experiments. Additionally, the stability properties for the
coupled tank states differ from those for single-tank states,
even when the adjustment for the change in effective flow
rates and input-stream reactant concentrations is made. This
is apparent from a plot of the coupled tank steady states on
single-tank plots, as in Figure 6.11; the coupled-tank marginal
stability points do not coincide with the single-tank marginal
stability points.

Finally, it is significant to notice that the eigen-
vector of the matrix g with eigenvalue zero has nonvanishing
components in both R

and R except for the case when one

1 2'
tank is at marginal stability prior to coupling. Thus, in the
general case, the stability of the c0upled tanks cannot be
localized in a single tank.

The analysis above has shown several new features for
coupled multiple steady state systems. First, an effective
potential and the shifts in location of the steady states
have been derived for model coupled folw tank reactors.
Second, it has been shown that the stability properties of
the coupled tank systems differ in general from single tank
stability properties; coupling stabilizes the tanks to small
concentration fluctuations and makes it possible to operate

coupled tanks in a region of the concentration space where

stable operation of single tanks is impossible.

172

 

 

 

 

 

 

(141; I l I I | I I l I I I I I I I I
040 - :
0.3 6 " r“ ‘
0.32 — :
A—fi .1
0.28 " ‘
l0 R - -
0.24 " "
0.20 " :
016 " ‘
0.12.. :
)— ...
0.08 - ..
0.04 " "
000 .. l l 1 l l 1 I l l I I l I I -

0.0 0.2 0.4 0.6 0.8 (.0 (.2 (.4 (.6

1043?.)

Figure 6.11 Comparision of steady state patterns of the

single and coupled flow reactors.

CHAPTER VII
FUTURE WORK AND DEVELOPMENT

Transitions between multiple steady states and.hysteresis
were studied using the birth-and-death master equation. This
does not include local fluctuations. A complete mesoscopic
description is provided by the multivariate master equation.
It is possible that the multivariate master equation leads to
slightly different results for these problems from those
presented here. Schlagl's coexistence condition4 is obtained
by including the effect of diffusion but neglecting
fluctuations in the particle number, whereas the Nicolis

condition30

is obtained by including number fluctuations and
neglecting spatial inhomogeneities. Grassberger38 has recently
argued that there is no sharp coexistence point when both
kinds of fluctuations are taken into account. If this is the
case, a study using the multivariate master equation may
indicate a static contribution to hysteresis. Hence such a
study seems worthwhile.

Since the eigenvector method and the simulation method
involve extensive computations, it is desirable to have a set
of evolution equations for the first few moments from the

birth-and-death master equation. By multiplying the master

equation by xn and summing, we get

173

174

n
d r1 _. n n—r _ r
E<X > _ Z [If] <x [a+(x) + ( l) a_(x)]> .
r=l
Since a+(x) and a_(x) are cubic in x, the equation for <xn>

. n+2
involves <x

>. A systematic approximation might be
developed by a truncation scheme.
Finally, it would be interesting to study the stochastic

and time-dependent behaviour of the coupled tank reactor

system.

APPENDICES

APPENDIX A

Proof that the transitions generating hysteresis cannot
occur before a marginal stability point is reached according

to the deterministic rate equation:

Eliminating t from equation (5.7), we obtain

——' = i F(XIB)/Bl (A°1)

where

C C
F(X,B) = —l 3x2 - —3 x3 + c3B - c4X.
2 6

The initial condition is specified by (XS,B) where XS is the
X-value at a stable steady state corresponding to B. Since
F(X,B) has a bounded derivative in a finite interval
(Lifshitz condition), the solution is unique. i.e., single-

valued (and continuous). Since 8 > O, (A.l) gives

dX

dB 0 at any steady state. (A.2)

The set of steady states are described by

F(XS(B)'B) = 0.

175

176

Considering the derivative of F along the steady state curve,
we find

a_F‘ .121) d__xs . o
33 x=x a x=x dB
S S
Therefore,
SEE = - 3: / 2: (A 4)
dB BB x=xs 3X x=xs
Also
3F) _ 2
331x=x - klxs + k3 > 0.

To find the sign 015°8-E

8 , we expand F at a steady state:

'x=x
(x-xs) + ... (A.5)

when B is time-independent. Thus

8F . .
3X < 0 if XS is stable, (A.6)

and hence from (A.4)

177

dX
—E > o (A.7)

dB

at all stable steady states.

Suppose there exists a solution which crosses the middle
(unstable) branch of the steady state curve. In order to
satisfy the initial condition the solution has to pass through
a stable steady state. This is possible only if its slope
exceeds that of the steady state curve at the point of
intersection. Since the slope of the stable branch is always
positive according to equation (A.7), this contradicts

equation (A.2). Thus the supposed solution does not exist.

APPENDIX B
Method and program to find the eigenvalues of the transition

matrix

The infinite state-space of the stochastic process is
made finite by setting a+(N) = O for a large N. Then the
characteristic polynomial of the tridiagonal matrix A is

given by fN(x) where

fk(x) = - [a(k) + x]fk_1(x) - a+(k-l)a_(k)fk_2(x)
for k = 1,2,...,N;
f_1(x) = 1 and f0(x) = — [a(O) + x].

Let Fk(x) = fk(x) / g(k) where g(k) is chosen for numerical

convenience (i.e., to avoid overflow) as
-k
g(k) = 2 a(k)a(k-l)...a(1)a(0).
Then

Fk(x) = - 2 [1 + x/a(k)]Fk_l(x)

a+(k-1) a_(k)
- 4 Fk_2(x).
a(k-l) a(k)

 

It is known that there are n roots of fk(x) = 0 in the

178

179

interval (-w,a) if fk(a) # O and the sequence f_l(a), f0(d),
..., fk(a) changes sign n times.62 Since g(k) > 0 for all k,
fk(a) can be replaced by Fk(a) in the above statement. This
lets one bracket the desired root suitably for solving the

above equation using any root-finding routine.

000000

000000

180'

PROGRAM EIGENVL
Subroutine ZEROIN should be supplied for loading

LOGICAL L(997:1001)

DIMENSION P(997:1001), R(998:1001)
EXTERNAL CHARPOL

COMMON /CB/ C(4), B, K

DATA C /0.3E-6, 0.1E-3, 0.15E-2, 1.7/

OPEN (60, FILE
OPEN (61, FILE

'INPUT', STATUS = 'OLD')
'OUTPUT', STATUS = 'UNKNOWN')

998 READ (60,*,END=999) C(4), BLOW, BHIGH, BSTEP
Do this until input records are exhausted

DO 40 B = BLOW, BHIGH, BSTEP

L(997) = .FALSE.
L(998) = .FALSE.
L(999) = .FALSE.
L(1000) = .FALSE.
L(1001) = .FALSE.

Hope that an eigenvalue will bite in the range between
-500 and -10. This need not be one of the few we need.
Since we will also know its position in the spectrum,
it can be used as a bait to catch the ones we want.

XA = - 500.0
X = - 10.0
XB = 0.0

First we find four ranges of x enclosed by P(i) and P(i-l)
such that each interval contains exactly one eigenvalue.
Note that K is common with CHARPOL and that Y is ignored
since we are only interested in X and K.

10 CONTINUE
Y = CHARPOL(X)
IF (K .GT. 1001) THEN
PRINT*, 'SIGN CHANGE OCCURRED ', K,

$ ' TIMES. B = ', B, ' X = ', X
STOP 'BAD'
ELSE IF (K .GE. 997) THEN
L(K) = .TRUE.
P(K) = X
END IF

IF ( .NOT. L(997) ) THEN
IF (K .LT. 997) XA = X
IF (K .GT. 997) X3 = X
X = (XA+XB) / 2.0
GO TO 10

0000

000

0 00000

181

ELSE IF ( .NOT. L(998) ) THEN
X = P(997) / 2.0
IF ( L(1000) ) X = P(997) + (P(1000)-P(997))/3.0
IF ( L( 999) ) X = P(997) + (P( 999)-P(997))/2.0
GO TO 10
ELSE IF ( .NOT. L(999) ) THEN
X = P(998) / 2.0
IF ( L(1000) ) X = P(998) + (P1000)-P(998))/2.0
GO TO 10
ELSE IF ( .NOT. L(1000) ) THEN
X = P(999) / 2.0
GO T010
ELSE IF ( .NOT. L(1001) ) THEN
X = 1.0
GO TO 10
END IF

Now the eigenvalues can be calculated to desired accuracy
using any root finding routine.

DO 30 I = 998, 1001
G1 = P(I-l)
GZ = P(I)
ABSERR = 1.0E-15
RELERR = 1.03-15
CALL ZEROIN ( CHARPOL, G1, G2, ABSERR, RELERR,

$ IFLAG)
IF (IFLAG .GT. 3) THEN
PRINT*, 'ERROR FROM ZEROIN. IFLAG = ',
$ IFLAG, ' B = ', B, ' G1 = ', G1,
$ ' 62 = ', G2
STOP 'BAD'
END IF
R(I) = G1
30 CONTINUE

Write results

WRITE (61,4) C(4), B, R

4 FORMAT ( '0', 4X, F7.5, 4X, 1PE10.4, 3(4X,E20.13),
$ 4X, E10.3)

40 CONTINUE
GO TO 998

Go back to see if there is another input record

999 CONTINUE

END
FNCTION CHARPOL (X)

This function routine evaluates the characteristic
polynomial of the transition matrix and counts the
number of sign changes during the evaluation

COMMON /CB/ C(4), B, KOUNT

C3B = C(3) * B
ClBBY2 = C(l) * B / 2 0
CZBY6 = C(Z) / 6 0
FKMl = 1.0
KOUNT = 0
AK = C3B
FK = - 1.0 - X/AK
IF ( FK . LT. 0.0 ) KOUNT = l
APK = C3B
DO 10 K = 1, 1000
FKMZ = FKMl
FKMl = FK
APKMl = APK
AKMl = AK
APK = ClBBY2 * K * (K-l) + C3B
AMK = CZBY6 * K * (K-l) * (K-Z) + C(4) * K
IF ( K .EQ. 1000 ) APK = 0.0
AK = APK + AMK
FK = - 2.0 * ( 1.0 + X/AK ) * FKMl
$ - 4.0 * APKMl/AKMl * AMK/AK * FKMZ

IF ( SIGN(1.0,FK) .NE. SIGN(1.0,FKM1) ) KOUNT = KOUNT+1
10 CONTINUE
CHARPOL = FK
IF ( CHARPOL .EQ. 0.0 ) KOUNT = 995
RETURN
END

APPENDIX C
Programs to construct the eigenvectors of the transition

matrix

The eigenvectors are constructed interactively because
optimum values are to be selected for several parameters at
intermediate stages.

An arbitrary value of Pj(0) is selected, values of Pj(x)
are calculated recursively for x > O, and finally gj is norma-
zed. The value for Pj(0) must be large enough not to cause
underflow of Pj(x) for some x and small enough not to cause
overflow of the normalization constant. The normalization
constant is obtained by summing Pj(x)2/Ps(x). This ratio as
calculated is not at all accurate for large values of x for
which both Pj(x) and Ps(x) are almost zero. Hence every
individual case has a specific cut-off value of x.

Every step involved in the computation of an eigen-
vector is performed by invoking an apprOpriate command; after
each step the results may be viewed, and if the results are
not satisfactory the step can be repeated with different
values for the parameters. The commands are:

INIT for initialization
DATA for changing the following values:
J — index of the current eigenvector
being computed. J = 1, 2, or 3.
IMAX - index of component beyond which

terms need not appear in sum while

183

P0

PJ

TERM

SHOW

SUM

NORM

SAVE

PLOT

TIME

STOP

HELP

184

calculating the normalization constant.
IMIN - index of component before which terms
need not appear in sum (often IMIN = O)
IMID - optional starting point for computing
Es (For x < IMID, 25(x) will be
calculated by backward iteration.)
for calculating unnormalized Es
for calculating unnormalized Bj
for calculating Pj(x)2/Ps(x) for all x (If the
result will cause overflow it is replaced by
- 1.0)
for a display of various vectors (EDIT with no
output) and the current values of the
parameters

for calculating the normalization constant for Ej

for normalizing gj orgs
for saving an eigenvector for plotting
for plotting all four eigenvectors

for plotting time-dependent asymptotic solutions

Since many components of the normalized Es are almost

zero, it is recommended that as be normalized after all other

eigenvectors have been normalized. However the normalizing

program performs correctly irrespective of whether the

normalized or unnormalized Es is used for normalizing gj.

185

In a typical run the commands may be issued in the
following order:

INIT (prompts for c and E values)

4
P0 (prompts for Ps(0))
For j = 1, 2, 3
DATA (set J = j)
PJ
TERM
SHOW
DATA (set IMAX)
DATA (set IMIN = 0)
SUM
NORM (J)
SAVE (J)
‘NORM (0)
SAVE (0)
PLOT

TIME

STOP

mmmmmmmm-m-mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm

186

A procedure to compute the eigenvectors and time-
dependent asymptotic distributions interactively.
On-line instructuions are expected to be self-explanatory.

ON ERROR THEN GOTO CMD

ON CONTROL_Y THEN GOTO CMD
ASSIGN SYS$COMMAND SYS$INPUT
ASSIGN SYS$OUTPUT OUTPUT

INQ = "INQUIRE/NOPUNCTUATION"
WJ = "WRITE JOB"

WS = "WRITE SYS$OUTPUT"

l

l

!

CMD:

INQ COM "COMMAND ?"

IF COM.EQS."STOP" THEN EXIT

IF COM.EQS."HELP" THEN GOTO HEL
IF COM.EQS."TIME" THEN GOTO TIM

There is a file EIGx.EXE for every valid response x,
except when x = STOP, HELP, or TIME

RUN EIG'COM'

IF COM.EQS."SHOW" THEN EDIT/READONLY EIGOUT.DAT

IF COM.EQS.”PLOT" THEN PRINT/NOFEED/DELETE VECPLOT.DAT
GOTO CMD

IM:

A command procedure is prepared and submitted to a batch
queue.

.- .- O- .- a C- O- .-

DELETE EIGTIME.COM;*

OPEN/WRITE JOB EIGTIME.COM

WJ "$ASSIGN SYS$INPUT INPUT"

WJ "SASSIGN SYS$OUTPUT OUTPUT"

WJ "$LINK/EXE=PGPLOT EIGFUN+[Z]PGPLOT+VPLOT+MACPLOT"
WJ "$RUN PGPLOT"

WJ "$DECK"
I

! Conversation with user
I

WS '3 "

INQ NP "HOW MANY PLOTS ? ( 1 ) "
INQ Y "ALL IN ONE PAGE ?"
INQ T1 "TIME FOR FIRST PLOT "

INQ T2 "TIME FOR LAST PLOT

mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm

INQ
INQ
INQ
IF DEL

187

XMAX "MAXIMUM X = "
PMAX "MAXIMUM P(x,t) EXPECTED = "
DELTA "IS INITIAL DISTRIBUTION A DELTA FUNCTION ? "

TA THEN INQ INPOS "PEAKED AT "

IF DELTA THEN GOTO SKIP

INP
W5
W8
W8
W8
W3
W8
W5
W8
INQ
INQ
INQ
KIP:

Prep
func
it i

o—o—o—o—o—U)

IF Y T
R = NP

08 = - 1

II II

"Then other acceptable initial distributions are"
"restricted to the span of the first four eigenvectors"
II II

"Enter components of the initial distribution along“
"these eigenvectors"

" A(1) = 1.0"
A2 "A(2) = "
A3 "A(3) = "
A4 "A(4) = "

are a data deck compatible with program PGPLOT and
tion PLFUN. The function reads data the first time
5 used in a job.

HEN WJ "l"
+ 1

IF .NOT.Y THEN WJ R

IF Y THEN NN = NP +

IF .NO
I

IP = 0

L1
IF
IF
IF
IF
WJ
IF

L2:

IP

X
IF IP.
WJ II s
WJ II s
WJ II $
CLOSE
WS II II
W5 "A

PJH

T.Y THEN NN

x = "YES"

IP.EQ.NP THEN N 100
IP.NE.NP THEN N IP
X THEN WJ " ' ', 1, 'X', l, 'P(X)', 4, "

X THEN WJ NN,", 0.0, ", XMAX, ", 0.0", ", ", PMAX

n oTRUE- I"! N! "IIIIIIIII"

IP.NE.O THEN GOTO L2

WJ II II’Tl' II’II’ T2, II'II' INPOS, II'II’ NP, II’II
IF .NOT.DELTA THEN WJ " ",A2,",",A3,",",A4,","

= IP + 1

= .NOT.Y

LE.NP THEN GOTO L1

EOD"

PRINT/NOFEED/DELETE PLOT.DAT"
DELETE PGPLOT.*;*"

JOB

batch job is ready to be submitted. Say NO if you"

WS "have made a mistake."

(hm-mmmmmmmmmmmmmmmmmmmmmmmmmmmmm

"SUBMIT?

188

IF OK THEN SUBMIT/QUEUE=FASTQ/NOPRINT EIGTIME.COM

"Valid responses are:"

W S II II

INQ OK

GOTO CMD
HEL:

W S II II

WS

W S II II

WS " INIT
WS " DATA
WS "

W S II

WS "

W S II

W S II

w S II

W 8 II

W 8 II

WS " P0
WS " PJ
WS " TERM
WS " SHOW
WS " SUM
WS " NORM
WS " SAVE
WS " PLOT
WS " TIME
WS " STOP"
WS " HELP"
W S II II

GOTO CMD

for initialization"
for changing the following values"
J - index of the current eigenvector"
being computed. J = 1,2,3."
IMAX - index of component beyond which"
terms need not appear in sum"
IMIN - index of component before which"
terms need not appear in sum"
IMID - optional starting point when"
computing P0"
for calculating unnormalized P0"
for calculating unnormalized PJ"
for calculating PJ**2/P0 for each component"
for a display of various vectors"
for calculating the normalization constant"
for normalizing PJ"
for saving current vector for plotting"
for plotting the normalized eigenvectors"
for plotting time dependent solutions"

GOO

GOO

000

189

PROGRAM EIGINIT

PARAMETER (NX=1500, NE=3)

IMPLICIT REAL*8 (A-H,O-Z)

DIMENSION APLUS(0:NX), AMINUS(0:NX), EV(0:NE), C(4),
$ P0(0:NX), PJ(0:NX), TERM(0:NX)

REAL SAVE(0:NX,O:NE)

DATA SAVE /6004*-1.0/, C /3D-7,1D-4, 1.5D-3, 1.7/

$ P0 /1501*-1D0/, PJ /1501*-1DO/, TERM /1501*-1D0/

WRITE (*,*) ' INIT'
WRITE (*,1) ' C4 = '
1 FORMAT (/'$', A6)
READ (*,*) CIN
WRITE (*,1) ' B = '
READ (*,*) BIN

Search through the eigenvalue listing

EV(O) = 0.0
OPEN (15, FILE='EIGVALS', FORM='FORMATTED',
$ STATUS = 'OLD')

WRITE (*,*) 'Searching for eigenvalues'
DO WHILE (C(4).NE.CIN .OR. B.NE.BIN)
READ (15,2,END=20) C(4), B, EV(3), EV(Z), EV(l)

2 FORMAT ( 5X, E7.0, 4X, E10.0, 3(4X,E20.0))
END DO
CLOSE (15)
J = 0
IMIN = -1
IMAX = -1
IMID = -1

Compute aplus and aminus and store in a file.

ClBBY2 = C(l) * B * 0.500
CZBY6 = C(2) / 6D0
C38 = C(3) * B
DO 10 I = o, NX
x = DBLE(I)

APLUS(I) = ClBBY2 * X * (x-lDO) + C33
AMINUS(I) = CZBYG * x * (x-an) * (X-ZDO) + C(4)*X
10 CONTINUE
OPEN (10,FILE='EIGDATA', STATUS='UNKNOWN',
$ FORM='UNFORMATTED')
WRITE (10) c, B, EV, APLUS, AMINUS, J, IMIN, IMAX, IMID

Erase numbers from previous work in the following files
OPEN (20, FILE='EIGPO', STATUS='UNKNOWN',

$ FORM='UNFORMATTED')
SUMO = -1D0

190

WRITE (20) P0,SUMO

CLOSE (20)

OPEN (3o, FILE='EIGPJ', STATUS='UNKNOWN',

s FORM='UNFORMATTED')

SUMJ = -1D0

PNORMJ = -100

WRITE (30) PJ, SUMJ, PNORMJ

CLOSE(30)

OPEN (4o, FILE='EIGTERM', STATUS='UNKNOWN',
$ FORM='UNFORMATTED')

WRITE (40) TERM

CLOSE (40)

OPEN (60, FILE='EIGSAVE', STATUS='UNKNOWN',
s FORM='UNFORMATTED')

WRITE (60) C(4), B, EV, SAVE

CLOSE (60)

STOP ' '

Normal termination at this point

20 CLOSE (15)
STOP 'Eigenvalues not available'
END

191

PROGRAM EIGDATA
PARAMETER (NX=1500, NE=3)
IMPLICIT REAL*8 (A-H, O-Z)

DIMENSION C(4), EV(O:NE), APLUS(0:NX), AMINUS(0:NX),

$ PJ(0:NX), TERM(0:NX)
CHARACTER VAR*4
DATA PJ / 1501*-1D0 /, SUMJ / -1DO /,
$ TERM / 1501*-1D0 /

WRITE (*,*) '

OPEN (10, FILE='EIGDATA', STATUS='OLD'
$ FORM='UNFORMATTED')

READ (10) C, B, EV, APLUS, AMINUS, J,
CLOSE (10)

WRITE(*,1) ' VARIABLE ? '
FORMAT ('$', All)

READ (*,2) VAR

FORMAT (A10)

WRITE (*,*) ' '

WRITE (*,3) VAR, ' = '
FORMAT ('S', A4, A3)

READ (*,*) IVAR

IF (VAR.EQ.'J' .OR. VAR.EQ.'j') THEN

J = IVAR

IMIN = —1

IMAX = -1

IMID = -1

OPEN (30, FILE='EIGPJ' , STATUS='
s , FORM='UNFORMATTED')

WRITE (30) PJ, SUMJ, PNORMJ

CLOSE (30)

OPEN (4o, FILE='EIGTERM', STATUS='
$ FORM='UNFORMATTED')

WRITE (40) TERM

CLOSE (40)
END IF

IF (VAR.EQ.'IMIN' .OR. VAR.EQ.'imin')
IF (VAR.EQ.'IMAX' .OR. VAR.EQ.'imaX')
IF (VAR.EQ.'IMID' .OR. VAR.EQ.'imid')

PNORMJ / -1DO /

DATA'

I

IMIN! IMAX,

UNKNOWN',

UNKNOWN',

IVAR
IVAR
IVAR

IMIN
IMAX
IMID

OPEN (10, FILE='EIGDATA', STATUS='UNKNOWN',

$ FORM='UNFORMATTED')

WRITE (10) C, B, EV, APLUS, AMINUS, J,
CLOSE (10)

END

IMIN, IMAX,

IMID

IMID

192

PROGRAM EIGPO

PARAMETER (NX=1500, NE=3)

IMPLICIT REAL*8 (A-H, O-Z)

DIMENSION P0(0:NX), APLUS(0:NX), AMINUS(0:NX), EV(0:NE),
$ C(4)

CHARACTER RES

WRITE( *,*) ' P0.
OPEN (10, FILE='EIGDATA', STATUS='OLD',
$ FORM='UNFORMATTED')

READ (10) C, B, EV, APLUS, AMINUS, J, IMIN, IMAX, IMID
CLOSE (10)
WRITE (*,1) ' STARTING FROM x = 0 ?'
1 FORMAT ('$', A21)
READ (*,2) RES
2 FORMAT (A1)
WRITE (*,1) 'STARTING CONSTANT = '
READ (*,*) CONST
IF ( RES.EQ.'Y' .OR. RES.EQ.'y' ) THEN
P0(O) = CONST
SUMO = P0(0)
DO 10 I = 1, NX
P0(I) = APLUS(I-l) / AMINUS(I) * P0(I-1)
SUMO = SUMO + P0(I)
10 CONTINUE
ELSE
IF (IMID.EQ.-1) STOP 'IMID not initialized'
P0(IMID) = CONST
SUMO = P0(IMID)
DO 20 I = IMID-1, O, -1
P0(I) = AMINUS(I+1) / APLUS(I) * PO(I+1)
SUMO = SUMO + P0(I)
20 CONTINUE
DO 30 I = IMID+1, NX
P0(I) = APLUS(I-l) / AMINUS(I) * P0(I-1)
SUMO = SUMO + P0(I)

30 CONTINUE
END IF
OPEN (20, FILE='EIGPO', STATUS='UNKNOWN',
$ FORM='UNFORMATTED')
WRITE (20) PO, SUMO
CLOSE (20)

END

193

PROGRAM EIGPJ

PARAMETER (NX=1500, NE=3)

IMPLICIT REAL*8 (A-H, O-Z)

DIMENSION PJ(0:NX), APLUS(0:NX), AMINUS(0:NX), EV(O:NE),
$ C(4)

WRITE(*,*) ' PJ'

OPEN (10, FILE='EIGDATA', STATUS='OLD',

$ FORM='UNFORMATTED')

READ (10)C, B, EV, APLUS, AMINUS, J, IMIN, IMAX, IMID
CLOSE (10)

OPEN (30, FILE='EIGPJ', STATUS='OLD',

$ FORM='UNFORMATTED')

READ (30) PJ, SUMJ, PNORMJ

CLOSE (30)

WRITE (*,1) 'PJ(O) = '

1 FORMAT ('$', A21)
READ (*,*) PJ(O)
PJ(1) = ( APLUS(O) + AMINUS(O) + EV(J) ) * PJ(0)/AMINUS(1)
DO 10 I = 2, NX

PJ(I) = ( ( APLUS(I-l) + AMINUS(I-1)) * PJ(I-l)
$ - APLUS(I-Z) * PJ(I-2) + EV(J) * PJ(I-l) )
$ / AMINUS(I)
10 CONTINUE

OPEN (30, FILE='EIGPJ', STATUS='UNKNOWN',
$ FORM='UNFORMATTED')

WRITE (30) PJ, SUMJ, PNORMJ

CLOSE (30)

END

194

PROGRAM EIGTERM

PARAMETER (NX=1500, NE=3)

IMPLICIT REAL*8 (A-H, O-Z)

DIMENSION P0(0:NX), PJ(0:NX), TERM(0:NX)

WRITE(*,*) ' TERM'
OPEN (20, FILE='EIGPO', STATUS='OLD',

$ FORM='UNFORMATTED')

READ (20) P0, SUMO

CLOSE (20)

OPEN (30, FILE='EIGPJ', STATUS='OLD',

$ FORM='UNFORMATTED')

READ (30) PJ, SUMJ, PNORMJ

CLOSE (30)

DO 10 I = 0, NX
IF ( DABS(PJ(I)) .GT. 1D19*DSQRT(P0(I)) ) THEN
TERM(I) = -1DO
ELSE
TERM(I) = PJ(I) / P0(I) * PJ(I)
END IF
10 CONTINUE

PNORMJ = DSQRT (SUMO)

OPEN (30, FILE='EIGPJ', STATUS='UNKNOWN',

$ FORM='UNFORMATTED')

WRITE (30) PJ, SUMJ, PNORMJ

CLOSE (30)

OPEN (40, FILE='EIGTERM', STATUS='UNKNOWN',
$ FORM='UNFORMATTED')

WRITE (40) TERM

CLOSE (40)

END

mum-m

2

mmmmmmmmmmm

195

PROGRAM‘EIGSHOW

PARAMETER

(NX=1500, NE=3)

IMPLICIT REAL*8 (A-H, O-Z)

DIMENSION C(4), EV(0:NE), APLUS(O:NX), AMINUS(0:NX),
P0(0:NX), PJ(0:NX), TERM(0:NX)

WRITE (*,*) ' SHOW'

WRITE (*,*) ' 1 : Data'

WRITE (*,*) ' 2 : PO'

WRITE (*,*) ' 3 : P0, Pj'

WRITE (*,*) ' 4 : P0, Pj, Term'

WRITE (*,1) 'SELE T A NUMBER : '

FORMAT (/ '5', A18)

READ (*,*) MODE

OPEN (10, FILE='EIGDATA', STATUS='OLD',
FORM='UNFORMATTED')

OPEN (20, FILE='EIGPO', STATUS='OLD',
FORM='UNFORMATTED')

OPEN (30, FILE='EIGPJ', STATUS='OLD',
FORM='UNFORMATTED')

OPEN (40, FILE='EIGTERM', STATUS='OLD',
FORM='UNFORMATTED')

READ (10) C, B, EV, APLUS, AMINUS, J, IMIN, IMAX,

READ (20) PO, SUMO
READ (30) PJ, SUMJ, PNORMJ
READ (40) TERM

CLOSE (10)
CLOSE (20)
CLOSE (30)
CLOSE (40)

OPEN (50, FILE='EIGOUT', STATUS='UNKNOWN',

FORM='UNFORMATTED')

IF (MODE.EQ.1) THEN

WRITE (*,2) C(4), B, EV, J, IMIN, IMAX, IMID
FORMAT (5X, 'C4 = ', F5.3 /
5X, '3 = ', 1PE11.4 /
5X, 'EV(0)= ', 0PF5.3 /
5X, 'EV(1)= ', 1PE11.4 /
5X, 'EV(2)= ', Ell.4 /
5X, 'EV(3)= ', Ell.4 /
5X, 'J = ', Il /
5X, 'IMIN = ', I4 /
5X, 'IMAX = ', I4 /
5X, 'IMID = ', I4 /
15X, 'The next line is irrelevant. ',
'Please type QUIT'l)
ELSE

IF (MODE.EQ.2) THEN
WRITE ( 50,3) (I, P0(I), I=O,NX )

$

END
END

196

FORMAT ( 5x, I4, 5x, 1PE11.4 /)
ELSE IF (MODE.EQ.3) THEN
WRITE (50,4) ( I, P0(I), PJ(I), I=O,NX)
FORMAT ( 5x, I4, 5x, 1PE11.4, 5x, Ell.4 / )
ELSE IF (MODE.EQ.4) THEN
WRITE (50,5) (I, P0(I), PJ(I), TERM(I), I=O,NX)
FORMAT ( 5x, I4, 5x, 1PE11.4, 5x, Ell.4, 5x,

END IF

WRITE (*,*)
WRITE (*,*)
WRITE (*,*)

WRITE (*,*)
IF

Ell.4 / )

I I

' Type C for viewing. '

' Type Control-Y (or Control-Z follow',
'ed by QUIT) to terminate SHOW'

I I

10

UNH

$

$

$

$

$

197

PROGRAM EIGSUM

PARAMETER (NX=1500, NE=3)

IMPLICIT REAL*8 (A-H, O-Z)

DIMENSION C(4), EV(O:NE), APLUS(O:NX), AMINUS(0:NX),
TERM(0:NX), PJ(0:NX)

WRITE (*,*) ' SUM'

OPEN (10, FILE='EIGDATA', STATUS='OLD',
FORM='UNFORMATTED')

READ (10) C, B, EV, APLUS, AMINUS, J, IMIN, IMAX, IMID

CLOSE (10)

OPEN (40, FILE='EIGTERM', STATUS='OLD',
FORM='UNFORMATTED')

READ (40) TERM

CLOSE (40)

OPEN (30, FILE='EIGPJ', STATUS='OLD',
FORM='UNFORMATTED')

READ (30) PJ, SUMJ, PNORMJ

CLOSE (30)

IF (IMIN.EQ.-1) STOP 'IMIN is not initialized'
IF (IMAX.EQ.-l) STOP 'IMAX is not initialized'

SUMJ = ODO

DO 10 I = IMIN, IMAX
IF (TERM(I).EQ.-1)STOP 'term negative'
SUMJ = SUMJ + TERM(I)

CONTINUE

WRITE (*,1) SUMJ
IF (IMIN.NE.O) WRITE (*,2) TERM(IMIN-l)
IF (IMAX.NE.NX) WRITE (*,3) TERM(IMAX+1)

FORMAT ( '0', 5X, 'Sum = ', 1PE11.4 )
FORMAT ( '0', 5X, 'previous term = ', 1PE11.4 )
FORMAT ( '0', 5X, 'next term = ', 1PE11.4 )

OPEN (30, FILE='EIGPJ', STATUS='UNKNOWN',
FORM='UNFORMATTED' )

WRITE (30) PJ, SUMJ, PNORMJ

CLOSE (30)

END

10

20

30

$

198

PROGRAM EIGNORM

PARAMETER (NX=1500, NE=3)
IMPLICIT REAL*8 (A-H, O-Z)
CHARACTER JA

DIMENSION P0(0:NX), PJ(0:NX)

WRITE (*,*) ' NORM'
CONTINUE
WRITE (*,1) ' 0 OR J ?'
FORMAT ('$', A10)
READ (5,2) JA
FORMAT (A1)
IF (JA.EQ.'0') THEN
OPEN (20, FILE='EIGPO', STATUS='OLD',
FORM='FORMATTED')
READ (20) PO, SUMO
CLOSE (20)
IF (SUMO.EQ.-1DO) STOP
' Normalization constant not available'
DO 20 I = O, Nx
P0(I) = P0(I) / SUMO
CONTINUE
SUMO = 1DO
OPEN (20, FILE='EIGPO', STATUS='UNKNOWN',
FORM='UNFORMATTED')
WRITE (20) P0, SUMO
CLOSE (20)
ELSE IF (JA.EQ.'J' .OR. JA.EQ.'j') THEN
OPEN (30, FILE='EIGPJ', STATUS='OLD',

FORM='UNFORMATTED')
READ (30) PJ, SUMJ, PNORMJ
CLOSE (30)
IF (PNORMJ.EQ.-1DO) STOP
'Normalization constant not available'
PNORMJ = PNORMJ * DSQRT (SUMJ)
DO 30 I = O, NX
PJ(I) = PJ(I) / PNORMJ
CONTINUE
PNORMJ = 1DO
SUMJ = 1DO
OPEN (30, FILE='EIGPJ', STATUS='UNKNOWN',
FORM='UNFORMATTED')
WRITE (30) PJ, SUMJ, PNORMJ
CLOSE (30)
ELSE
GO TO 10
END IF

END

199

PROGRAM EIGSAVE

PARAMETER (NX=1500,NE=3)

REAL*8 P(O:NX), SUM, PNORM, C(4), C4, B, EV(O:NE),
$ APLUS(O:NX), AMINUS(0:NX)

REAL SAVE(O:NX,O:NE)

CHARACTER JA

WRITE (*,*) ' SAVE'
OPEN (60, FILE='EIGSAVE', STATUS='OLD',

$ FORM='UNFORMATTED')

READ (60) C4, B, EV, SAVE

CLOSE (60)

WRITE (*,1) '0 OR J 2'
1 FORMAT ('$', A10)
READ (*,2) JA
2 FORMAT (A1)
IF (JA.EQ.'0') THEN
OPEN (20, FILE='EIGPO', STATUS='OLD',
s FORM='UNFORMATTED')
READ (20) P, SUM
IF (SUM.NE.1D0) STOP 'Not normalized'
DO 10 I = 0, Nx
SAVE(I,O) = SNGL (P(I))
10 CONTINUE
ELSE IF (JA.EQ.'J' .OR. JA.EQ.'j') THEN
OPEN (10, FILE='EIGDATA', STATUS='OLD',
$ FORM='UNFORMATTED')
READ (10) C, B, Ev, APLUS, AMINUS, J, IMIN, IMAX,
$ IMID

CLOSE (10)

OPEN (30, FILE='EIGPJ', STATUS='OLD',
s FORM='UNFORMATTED')

READ (30) P, SUM, PNORM

CLOSE (30)

IF (PNORM.NE.1DO) STOP 'Not normalized'
DO 20 I = 0, NX
SAVE(I,J) = SNGL (P(I))

20 CONTINUE
END IF
OPEN (60, FILE='EIGSAVE', STATUS='UNKNOWN',
$ FORM='UNFORMATTED')
WRITE (60) C4, B, EV, SAVE
CLOSE (60)

END

0000

000

000

200

PROGRAM EIGPLOT

Files [Z]VPLOT.OBJ and [Z]MACPLOT.OBJ should be provided
for LINKing.

PARAMETER (NX=1500, NE=3)
REAL*8 C4, B, EV(O:NE)
REAL P(O:NX,0:NE)

LOGICAL DONE(O:NE)
CHARACTER PLBL*10, ANS*3

WRITE (*,*) ' PLOT'
WRITE (*,1)
$ 'MAXIMUM VALUE OF x TO APPEAR ON PLOT AXIS = '

1 FORMAT ('$', A45)

READ (*,*) LIMIT

OPEN (60, FILE='EIGSAVE', STATUS='OLD,
$ FORM='UNFORMATTED')

READ (60) C4, B, EV, P

CLOSE (60)

DO 10 J = 0, NE
DONE(J) = P(0,J).NE.-1.0
IF (.NOT.DONE(J)) THEN
WRITE (*,2) J
2 FORMAT (/1x, 'P(', 11, ') is not available')
WRITE (*,3) 'OK ?'
3 FORMAT ('$', A4)
READ (*,4) ANS
4 FORMAT (A3)
IF (ANS(:1).NE.'Y' .AND. ANS(:1).NE.'y') STOP
$ 'Ignore the following message:'
END IF
10 CONTINUE

OPEN (61, FILE='OUTPUT', STATUS='UNKNOWN',
$ FORM='FORMATTED')
OPEN (71, FILE='VECPLOT', STATUS='UNKNOWN')

Each eigenvector is drawn separately

DO 60 J = 0, NE
IF (DONE(J)) THEN

Y-axis limit is determined below

PMAX = 0.0
DO 20 I = O, LIMIT
AP = ABS(P(I,J))
IF (AP.GT.PMAX) PMAX = AP
20 CONTINUE

 

000 O

000

000

(300

Scales

201

WRITE (71,5) C4, B, EV(J)

FORMAT ('1', 30X, 'C4 = ', F4.2, 5X, 'B = ',
1PE11.4, 5X, 'EIGENVALUE = ', Ell.4 )

CALL PLOTS (9.5, 12.0, 71)

for the two axes are calculated below

XSCALE = FLOAT (LIMIT) /11.0
IF ( J.EQ.O ) THEN

YL = 0.0
YH = PMAX
ELSE
YL = - PMAX
YH = PMAX
END IF
YSCALE = (YH - YL) / 8.0

The axes are drawn

CALL PLOT (0.5, 1.0, -3)
CALL PLOT (8.0, 0.0, 2)
CALL PLOT (8.0, 12.0, 2)

The X-axis is marked

30

Y-axis

‘40

DEL LIMIT / 10
AX 0.0
DO 30 AN = 0.0, FLOAT(LIMIT), DEL
CALL PLOT ( 8.0, AX, 3)
CALL PLOT ( 8.1, AX, 2)
CALL NUMBER (8.2, AX-0.1, 0.1, AN, 90.0, -1)
AX = AX + DEL / XSCALE
CONTINUE
CALL SYMBOL ( 8.6, 5.45, 0.1, 'X', 90.0, 1)

is marked

IEXP = LOG10(YH)

DEL = FLOAT(INIT((YH-YL)/10.0**(IEXP-1)))*
l0.0**(IEXP-2)

AY = 8.0

DO 40 AN = YL/10.0**IEXP, YH/10.0**IEXP,

DEL/10.0**IEXP
CALL PLOT (AY, 0.0, 3)
CALL PLOT (AY, -0.1, 2)
CALL NUMBER (AY, -0.7, 0.1, AN, 90.0, 3)
AY = AY - DEL/YSCALE
CONTINUE
IF (IEXP.LT.-9) THEN
PLBL = 'P(X)*10“' // CHAR (-IEXP/10+48) //
CHAR ( MOD(-IEXP,10) + 48 )
NLBL = 10

 

000

GOO

000 O

202

ELSE
PLBL

'P(X)*10“' // CHAR (~IEXP+48)

NLBL = 9
END IF

Let the user know what is happening at the moment

6
$

WRITE (61,6) J
FORMAT ( // '0', 5X, 'INITIALIZATION FOR PLOT',
I1, 'COMPLETED' )

An eigenvector is now drawn

50

Inform

7
END

XIN 0.0
YIN 8.0 - (P(0,J)-YL) / YSCALE
CALL PLOT ( YIN, XIN, 3)

DO 50 I = 1, LIMIT
XIN = FLOAT (I) / XSCALE
YIN = 8.0 - (P(I,J)-YL) / YSCALE

CALL PLOT ( YIN, XIN, 2 )
CONTINUE

CALL PLOT (0.0, 0.0, 999)
the user
WRITE (61,7) J

FORMAT ( '0', 5X, 'PLOT ', Il, ' COMPLETED')
IF

60 CONTINUE

END

20

30

203

FUNCTION PLFUN (IND, X)

PARAMETER (NX=1500, NE=3)

REAL*8 C4, B, EV

LOGICAL FIRST

DIMENSION A(O:NE), EV(O:NE), P(O:NX,0:NE)
SAVE P, FIRST T1, TSTEP

DATA FIRST /.TRUE./

IF (FIRST) THEN
FIRST = .FALSE.
READ (*,*) T1, T2, INPOS, N
TSTEP = (T2-T1) / (N—l)

OPEN (50, FILE='EIGSAVE', STATUS='OLD',
FORM='UNFORMATTED')

READ (50) C4, B, EV, P

CLOSE (50)

A(O) = 1.0
IF (INPOS.LT.0) THEN

READ (*.*) ( A(J). J=1.NE )
ELSE

DO 20 J = 1, NE
A(J) = P(INPOS,J) / P(INPOS,O)
CONTINUE
END IF
END IF

PLFUN = 0.0
T = T1 + IND * TSTEP
DO 30 J = 0, NE
I = INT (X)
PX = (I+1-X)*P(I,J) + (X-I)*P(I+1,J)

PLFUN = PLFUN + A(J)*PX * EXP (SNGL(EV(J)) * T )

CONTINUE

RETURN
END

APPENDIX D

Deterministic Simulation of hysteresis

0000

0000

00000

PROGRAM DETERM

Files ZEROIN8.0BJ and [Z]ADAMSPC.OBJ should be provided
for LINKing

PARAMETER (N=1)

IMPLICIT REAL*8 (A-H, O-Z)

EXTERNAL AUX, F1

DIMENSION B(O:1000), XF(0:1000), XB(O:1000), C(4),
$ X(N,8), SAVE(N,11), YMAX(N), ERROR(N)
LOGICAL FORW

COMMON /Z/ BZ

COMMON /BC/ BETA, C

DATA C /3.0D-7, 1.0D-4, 1.5D-3, 1.7D0 /

OPEN (60, FILE='INPUT', FORM='FORMATTED', STATUS='OLD')
OPEN (61, FILE='OUTPUT', FORM='FORMATTED',
$ STATUS='UNKNOWN')

READ (60,*) C(4), B1, BZ, BETA
WRITE (61,*) C(4), B1, B2, BETA

1 FORMAT ('1', 10X, 'C(4) = ', F5.3 /
$ '0', 10X, 'B1 = ', 1PE11.4 /
$ '0', 10X, 'B2 = ', Ell.4 /
$ '0', 10X, 'BETA = ',_E11.4 )

The B-interval is divided into 1000 subintervals and data
are collected at the dividing points.

BSTEP = (BZ - B1) / 1000.0
B(O) = Bl
DO 10 I = l, 1000
B(I) = B(I-l) + BSTEP
10 CONTINUE

Two integrations are needed in the same interval. First in
the forward direction, (i.e., with B increasing; the loop
control variables are adjusted accordingly

BZ B(O)
ISTART = l
IEND = 1000
IDIFF = 1
H = 1.0D-12
Set initial condition (B = BZ, X = X(1,1) ) depending on

204

205

C whether forward or backward variation of B.

C
20 CONTINUE
G1 = 10.0
G2 = 1500.0
CALL ZEROIN8 (F1, G1, G2, 1.0E-8, 1.0E-8, IFLAG)
IF (IFLAG.GT.2) THEN
WRITE (61,2) IFLAG, BZ, G1, G2
2 FORMAT ( 'OERROR FROM ZEROIN'
$ '0', ' IFLAG = ', I1 /
$ '0', ' B = ', 1PE11.4 /
$ '0', ' G1 = ', Ell.4 /
$ '0', ' 62 = ', Ell.4 )
STOP 'ERROR'
END IF
C
IF (FORW) THEN
XF(O) = G1
ELSE
XB(1000) = Cl
END IF
X(1,1) = G1
HMIN = 1.0E-18
EPS = 1.0E-8
JSTART = 0
C
DO 30 I = 1, N
YMAX(I) = 1.0
30 CONTINUE
C
C Integration loop
C
DO 40 I = ISTART, IEND, IDIFF
CALL ADAMS (AUX, N, BZ, X, SAVE, H, HMIN, B(I),
$ EPS, YMAX, ERROR, KFLAG, JSTART)
IF (KFLAG.LT.0) THEN
WRITE (61,3) KFLAG, FORW, BZ, B(I)
3 FORMAT ('0', 5x, 'ERROR FROM ADAMS' /
$ '0', 10X, 'KFLAG = ', I2 /
$ '0', 10X, 'FORW = ', L10 /
$ '0', 10X, '32 = ', E10.4 /
$ '0', 10X, 'B(I) = ', E10.4 )
STOP 'ERROR'
END IF
C
C Collect data after each step

C

IF (FORW) THEN
XF(I) = X(l,1)

ELSE

END IF

4O CONTINUE

0 000

000

0000

206

Integration is complete

IF (FORW) THEN
prepare for the backward sweep
FORW = .FALSE.

B2 = B(lOOO)
ISTART = 999
IEND = 0
IDIFF = -1
H = - H
BETA = - BETA
GO TO 20

END IF

Write data for plotting

OPEN (71, FILE='DETPLOT', STATUS='UNKNOWN',

$ FORM='FORMATTED')

WRITE (71,4) (B(I), XF(I), XB(I), I=0,1000)
4 FORMAT (1001(3(5X,E10.4)/))

Area inside the loop is computed by simplest (rectangular)
method

AREA = 0.0
DO 50 I = 0, 1000
AREA = AREA + BSTEP*(XB(I) - XF(I))
50 CONTINUE
WRITE (61,6) AREA
6 FORMAT ( '0', 10X, 'AREA = ', 1PE11.4)
END

0000

0000

207

SUBROUTINE AUX (B, X, XD, N)
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION X(N), XD(N), C(4)
COMMON /BC/ BETA, C

This routine computes dB/dX. It is used indirectly through
the library subroutine ADAMS

XD(l) = F1 (X(1)) / BETA
RETURN
END

REAL*8 FUNCTION F1(X)
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION C(4)

COMMON /Z/ B

COMMON /BC/ BETA, C

This function computes dX/dt. It is called by AUX and
indirectly used by program DETERM through ZEROIN8.

F1 = ( 0.5*C(l)*X*(X-1.0)+C(3) ) * E

$ -C(2)/6.0 *X*(X-1.0)*(X-2.0) - C(4)*X
RETURN

END

APPENDIX E

Stochastic simulation of hysteresis

PROGRAM GILSIMT
C File ZEROIN8.0BJ Should be supplied for LINKing.

DOUBLE PRECISION x, B, APLUS, AMINUS, AMBYAP, QRAN,

$ FACTOR, BETA, C1, c2, c3, c4, CZBY6,
$ C1BY2

EXTERNAL F1

PARAMETER (IDIM=2001)

DIMENSION XSUM(0:IDIM-1,2), BGRID(O:IDIM-1), A(4),

$ XNOW(0:IDIM-1), AR(2)

LOGICAL FORW, GET, PUT, PLOT

COMMON /AA/ Bz, C(4)

C
OPEN (60, FILE='INPUT', STATUS='OLD')
OPEN (61, FILE='OUTPUT', STATUS='UNKNOWN')
C
READ (60,*) Bl, BZ, C, BETA, NTIMES, ISEED, FORW, GET,
$ ' PUT, PLOT
C1 = C(1)
C2 = C(2)
C3 = C(3)
C4 = C(4)
WRITE (61,1) B1, BZ, C, BETA, NTIMES, ISEED, GET, PUT,
$ PLOT
1 FORMAT ( '1', 1P, 10X, 'B1 = ', Ell.4 /
$ '0', 10X, '82 = ', Ell.4 /
$ '0', 10X, 'C = ', 4(E11.4,3X) /
$ ‘0', 10X, 'BETA = ', Ell.4 ///
$ '0', 10X, 'NTIMES = ', IlO /
$ '0', 10X, 'ISEED = ', IlO /
$ '0', 10X, 'GET = ', L10 /
$ '0', 10X, 'PUT = ', L10 /
s '0', 10x, 'PLOT - = ', L10 /// )
C
C The B-interval is divided into IDIM-l subintervals and
C data will be collected at the dividing points.
C
BGRID(O) = B1
BSTEP = (B2-Bl) / FLOAT(IDIM-l)
DO 10 I = 1, IDIM-l
BGRID(I) = BGRID(I-l) + BSTEP
10 CONTINUE
C
C This program can be run either to refine existing data
C statistically, or to acquire new data. In the former case

208

209

C data should be read and in the latter case initial

C
C

0000

conditi

IF (

20
ELSE

3O

Startin
steady

on should be set.

GET) THEN
OPEN (70, FILE='GRPREV', STATUS='OLD',
FORM='UNFORMATTED')
READ (70) NPREV, XSUM, ASUM, AZSUM
CLOSE (70)
DO 20 I = 1, IDIM-l
XSUM(I,1) = XSUM(I,1) * NPREV
XSUM(I-1,2) = XSUM(I-1,2) * NPREV
CONTINUE

DO 30 J = 1, 2
DO 30 I = 0, IDIM-l
XSUM(I,J) = 0.0

CONTINUE
ASUM = 0.0
AZSUM = 0.0
NPREV = 0

g values of X are Obtained by solving the cubic
state equation.

G1 = 1.0

G2 = 1500.0
B2 = BGRID(O)
ABSERR = 0.0
RELERR = 1.0E-5

Note that 82 appears in common block AA.
CALL ZEROIN (F1, G1, 62, ABSERR, RELERR, IFLAG)
IF (IFLAG.GT.2) THEN

WRITE (61,2) G1, G2, IFLAG

FORMAT ('OERROR FROM ZEROIN. G1 = ', F8.4,
' G2 = ', F8.4, ' IFLAG = ', Il)

STOP 'ERROR'

END IF

XSUM(0,1) = G1

Gl = 1.0

G2 = 1500.0

B2 = BGRID(IDIM-l)

Note that BZ appears in common block AA
CALL ZEROIN (F1, G1, G2, ABSERR, RELERR, IFLAG)
IF (IFLAG.GT.2) THEN

WRITE (61,2) G1, G2, IFLAG

STOP 'ERROR'

END IF
C1BY2 = C1 / 2.0
C2BY6 = C2 / 6.0

(3()0(3()0 ()0

0

(3(30 ()0

210

The outermost lOOp begins *************************

DO 70 N = 1, NTIMES

The same B-interval is covered twice. J = 1 means B is
increasing. J = 2 means B is decreasing. Variables
controlling the innermost D0 are adjusted in accordance

with J.

DO 65 J = 1, 2

BETA = ABS (BETA)

IF (J.EQ.1) THEN
B = BGRID(O)
X = XSUM(O,1)
XNOW(O) = X

ISTART = O
IEND = IDIM - 1
IDIFF = 1
ELSE
BETA = - BETA

B = BGRID(IDIM-l)
X = XSUM(IDIM-1,2)
XSUM(IDIM-l) = X
ISTART = IDIM-1
IEND = O
IDIFF = -1

END IF

APLUS = ( C1BY2 * X * (X-1.0) + C(3) ) * B
AMINUS = ( C2BY6 * (X-1.0) * (X—2.0) + C(4) ) * X
These must have some values to begin with.

DO 60 I = ISTART+IDIFF, IEND, IDIFF

************************ The simulation section begins

50

mm

CONTINUE
QRAN = 1DO - 2D0 * BETA / APLUS / B *

ALOG (1.0-RAN(ISEED) )
AMBYAP = AMINUS / APLUS

IF ( ( J.EQ.1 .OR. QRAN .GT. -2D0*AMBYAP )
.AND. DABS(QRAN/AMBYAP+2DO) .GT.
AMBYAP*lD-5 ) THEN

FACTOR = - AMBYAP + DSQRT ( QRAN +

(2DO+AMBYAP) * AMBYAP )

ELSE

FACTOR = 1D0 + 0.5D0 * QRAN / AMBYAP
END IF
B = B * FACTOR
APLUS = APLUS * FACTOR
IF (1.0/(l.0-RAN(ISEED)).GE.1.0+AMBYAP) THEN

APLUS = APLUS + C1*B*X

AMINUS = AMINUS + 3DO * C2BY6 * X * (X-lDO) + C4

(10(UFJO

()0(3()0

(10(1 0(3()0(3

(3(30

211

X = X + lDO
ELSE

X = X - 1DO

APLUS = APLUS - C1*B*X

AMINUS = AMINUS - 3D0*C2BY6*X*(X-1DO) - C4
END IF
IF (B.LT.BGRID(I) .EQV. J.EQ.1) GO TO 50

************************ simulation section ends.
Collect data after each step

XSUM(I,J) XSUM(I,J) + X
XNOW(I) =

60 CONTINUE

X

A one-way trajectory has been completed. Find area under
this trajectory and update sums involved in average and
standard deviation.

AR(J) = 0.0
DO 62 K = 1, IDIM-1
AR(J) = AR(J) + 0.5*(XNOW(K)+XNOW(K-1))*BSTEP
62 CONTINUE
65 CONTINUE
AREA AR(2) - AR(1)

ASUM ASUM + AREA
AZSUM = AZSUM + AREA*AREA
7O CONTINUE

The outermost loop ends ******************************
Sums are converted to averages.

NTOTAL = NPREV + NTIMES
DO 80 I = 1, IDIM-1
XSUM(I,1) = XSUM(I,1) / NTOTAL
XSUM(I-1,2) = XSUM(I-1,2) / NTOTAL
8O CONTINUE

Save results if desired.

IF (PUT) THEN

OPEN (71, FILE='GRNOW', STATUS='NEW',
$ FORM='UNFORMATTED')

WRITE (71) NTOTAL, XSUM, ASUM, AZSUM
END IF

Write a few lines on the output.
AVEG ASUM / NTOTAL

SDEV SQRT (AZSUM/(NTOTAL-l) - ASUM/NTOTAL*ASUM/
$ (NTOTAL-1))

0000

00000

212

WRITE (61,3) AVEG, SDEV, SDEV/AVEG

3 FORMAT ( '0', 'AREA = ', 1PE11.4 //
z :3? :SDEV = 1' FEE-z x
, RATIO= PE .
$ '0', '********;**************I )

Write formatted output if desired. These can be later
read by a plotting program

IF (PLOT) THEN
OPEN (72, FILE='PLOTFL', STATUS='NEW',
$ FORM='FORMATTED')
WRITE (72,4) (BGRID(I), XSUM(I,1), XSUM(I,2),
S I = 0, IDIM-1 )
4 FORMAT ( 5X, 1PE11.5, 5X, E11.5, 5X, E11.5 )
END IF

END

FUNCTION F1(X)

The cubic equation used for setting initial conditions.
Note that program GILSIMT uses this indirectly through
subroutine ZEROIN.

COMMON /AA/ B, C(4)

F1 = C(1) * B * X * (X-1.0) / 2.0 - C(2) * X * (X-1.0)
$ * (X-2.0) / 6.0 + C(3) * B - C(4) * X

RETURN

END

APPENDIX F
Special features of the cubic mechanism in a coupled flow

tank reactor

Let v(R1,RO) be a cubic rate law, supporting three
steady states in a flow tank. Let k0 be chosen so that

RB - RY = RY - Ra. Let R+ and R_ define the points satisfying

3V _

Then

i) ¢(Rl,R2;k0,kx) is symmetric with respect to the line
R1 + R2 = 2 Ry,

ii) R+ + R_ = 2 RY'

iii) (R+,R_) is a solution of the coupled tank equations
(6.7),

and iv) the state (R+,R_) is marginally stable.

Proof:

The potential for the coupled tank satisfies

8O de
- —— = -——-= k (R -R ) - V(R ,R ) + k (R -R )
3R1 dt 0 0 1 1 0 x 2 1
BO dR2
and - -—— = ——— = k0(RO-R2) - v(R2,R0) + kX(R1-R2). (F.2)
8R2 dt

Since Ra' R and RY are the single-tank steady states,

8!
213

214

k0 (RO-R) - V(R,R0) = C (R-Ra) (R-RB) (R-RY)
=C[R3-(R+R+R)R2+(RR+RR+RR)R
a B Y a B a Y B Y
- RaRBRY]' (F.3)

The point obtained by reflecting (R1,R2) on the line

R1 + R2 = 2 RY Is (ZRY-RZ' 2Ry-R1). Let us compare 0(R1,R2)

and 2R -R 2R -R .
¢( Y Y 1)

2i

¢(R1,R2) - ¢(2Ry-R2, 2RY-R1) = f(Rl) + f(RZ), (F.4)

where

- - _ 1 4 _ _ 4 l 3 -
f(R) — c { 4 [R (2RY R) 1 + 3 (Ra+RB+Ry)[R
3 1 2 2
(ZRY-R) ] - 5 (RGRB+RaRy+RBRY) [R - (2Ry-R) ]
+ RaRBRy[R - (2Ry-R)]} (F.S)

Expanding f(R), it is found that the coefficients of all

powers of R vanish when RB-RY = RY-Ra. Therefore,

¢(R1,R2) = ¢(2Ry-R2, 2Ry-R1). (F.6)

This proves (i).

From (F.1) and (F.2) it follows that R+ satisfy

2 -
3 R1 - 2 (Ra + R8 + Ry)Ri + (RGRB + RQRY + RBRy) — 0.

(F.7)

215

Therefore,

_ Z
R + R - 3(Ra + R8

8
ii) is Proved.

since R + R = 2R .
a Y

Solving (F.7), we Obtain

1
= + — -
Rt RY _ /3(RY RGRB) .
Let R - R = R ’ R = A.
B Y Y a
l
= + I
Then R:t RY - 5 (.3 .

(R+,R_) is a coupled steady state if

k0(R0-R+) - V(R+,R0) + kx(R_-R+)

and

ko(R0—R_) — V(R_,R0) + kX(R+-R_)

(F.12) can also be written as

c (R+-Ra)(R+-RB)(R+-Ry) + kX(R_-R+)

 

2
i.e., - -_—: CA - k

+ R ) = 2R
Y Y

(F.8)

(F.9)

(F.10)

(F.11)

(F.12)

(F.13)

(F.14)

216

Therefore, (R+,R_) is a steady state when
= —E <
k 3 A . (c 0)

Thus (iii) is proven.

The condition for marginal stability is

(

3v] 3v

k0 + BRlR k0 + OR
a I

 

 

 

 

(R+,R_) satisfies this condition by construction

(F.15)

(F.16)

(F.1).

REFERENCES

REFERENCES

1. G. Nicolis and I. Prigogine, "Self-organization in
Nonequilibrium Systems" (John-Wiley, New York, 1977).

2. A. A. Andronov, A. A. Vitt, and S. E. Khaikin, "TheOry of
Oscillations" (Pergamon, Oxford, 1966).

3. F. Scthgl, Z. Phys. 248, 446 (1971).

4. F. Scthgl, Z. Phys. 253, 147 (1972).

5. I. Matheson, D. F. Walls, and C. W. Gardiner, J. Stat.
Phys. 12, 21 (1975).

6. D. A. Sanchez, "Ordinary Differential Equation and
Stability Theory: An Introduction" (Dover, New York, 1968).

7. J. L. Doob, "Stochastic Processes" (John-Wiley, New York,
1953).

8. I. Oppenheim, K. E. Shuler, and G. H. Weiss, Adv. Mol.
Relax. Proc. 1, 13 (1968).

9. D. T. Gillespie, Physica 255, 69 (1979).

10. D. T. Gillespie, J. Comp. Phys. 22, 403 (1976).

11. H. A. Kramers, Physica, Z, 284 (1940); J. E. Moyal, J. Roy.
Stat. Soc. 82, 150 (1949).

12. A. D. Fokker, Ann. Phys. 43, 810 (1914); M. Planck,
Sitzber. Preuss. Acad. Wissensch. 1917, p.324.

13. N. G. van Kampen, Adv. Chem. Phys. 34, 245 (1976).

14. R. Kubo, K. Matsuo, and K. Kitahara, J. Stat. Phys.
2, 51 (1973).

15. D. T. Gillespie, J. Chem. Phys. 12, 5363 (1980).

217

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

27.

28.

29.

30.

218

K. Matsuo, K. Lindenberg, and K. E. Shuler, J. Stat. Phys.
22, 65 (1978).

L. Arnold, "Stochastic Differential Equations: Theory
and Applications", (John-Wiley, New York, 1974).

S. Chandrasekhar, Rev. Mod. Phys. 22, 1 (1943).

P. G. Hoel, S. C. Port and C. J. Stone, "Introduction
to Stochastic Processes" (Houghton-Mifflin, Boston, 1972).
D. T. Gillespie, J. Comp.Phys. 22, 395 (1978).

L. Onsager and S. Machlup, Phys. Rev. 22, 1505 (1953);
S. Machlup and L. Onsager, Phys. Rev. 22, 1512 (1953).
K. L. C. Hunt and J. Ross, J. Chem. Phys. 12, 976

(1981) and references therein.

For example, U. Weiss, Z. Phys. B 22, 429 (1978);

J. D. Bekenstein and L. Parker, Phys. Rev. 222, 2850
(1981); H. Hara, Z. Phys. B 22, 159 (1981).

A. Nitzan, Phys. Rev. A 21, 1513 (1978).

C. Escher and J. Ross, J. Chem. Phys. 12, 3773 (1983).
F. SCthgl, C. Escher, and R. S. Berry, Phys. Rev. A
21, 2698 (1983).

D. FOrster, "Hydrodynamic fluctuations, Broken Symmetry
and Correlation Functions" (W. A. Benjamin Inc. Reading,
Mass., 1975) Ch.7.

K. Matsuo, J. Stat. Phys. 22, 169 (1977).

R. M. Noyes, J. Chem. Phys. 12, 5144 (1979).

G. Nicolis and J. W. Turner, Ann. N. Y. Acad. Sci.

316, 251 (1979).

31.

32.

33.

34.

35.

36.

37.

38.

39.

40.

41.

42.

43.

44.

45.

46.

219

I. Oppenheim, K. E. Shuler, and G. H. Weiss, Physica
222, 191 (1977).

I. Procaccia and S. Mukamel, J. Stat. Phys. 22, 633
(1978).

I. Procaccia and J. Ross, J. Chem. Phys. 21, 5558 (1978).
R. D. Levine J. Chem. Phys. 22, 3302 (1976).

I. Procaccia and R. D. Levine, J. Chem. Phys., 22,

3357 (1976).

D. T. Gillespie, Physica 222, 69 (1979).

V. Seshadri, B. J. West, and K. Lindenberg, J. Chem.
Phys. 12,1146 (1980).

P. Grassberger, Z. Phys. B 21, 365 (1982).

M. Malek-Mansour and G. Nicolis, J. Stat. Phys. 22,

197 (1975).

D. T. Gillespie, J. Stat. Phys. 22, 311 (1977).

D. T. Gillespie, J. Phys. Chem. 22, 2340 (1977).

W. Horsthemke, M. Malek-Mansour, and B. Hayez. J. Stat.
Phys. 22, 201 (1977).

G. Nicolis and M. Malek-Mansour, J. Stat. Phys. 22,

495 (1980).

H. E. Stanley, "Introduction to Phase Transition and
Critical Phenomena" (Oxford University Press, New York, 1971).
G. Dewel, D. Walgraef, and.Ph Borckmans, Z. Phys. B. 22J
235 (1977).

D. Walgraef, G. Dewel, and.Pu Borckmans,‘Adv. Chem. Phys.

Q, 311 (1982).

220

47. S. K. Ma, "Modern Theory of Critical Phenomena"
(W. A. Benjamin, Inc., Reading, Mass., 1976).
48. F. Scthgl, Z. Phys. B 2, 301 (1975).

49. H. K. Janssen, Z. Phys. 2_2, 67 (1974).

50. C. W. Gardiner, K. J. McNeil, D. F. Walls, and
I. S. Matheson, J. Stat. Phys. 22, 309 (1976).

51. S. Grossmann and R. Schranner, Z. Phys. B 22, 325
(1978).

52. M. Mangel, J. Chem. Phys. 22, 3697 (1978).

53. C. L. Creel and J. Ross, J. Chem. Phys. 22,3779 (1976).

54. H. A. Kramers, Physica 1, 284 (1940).

55. C. W. Gardiner, J. Stat. Phys. 22, 157 (1983).

56. A. Lotka, Proc. Nat. Acad. Sci. USA. 2, 420 (1920).

57. I. Prigogine and R. Lefever, J. Chem. Phys. 22,1695
(1968).

58. J. Tyson, J. Chem. Phys. 22,3919 (1073).

59. R. J. Field and R. M. Noyes, J. Chem. Phys. 22, 1877
(1974).

60. RI Field, J. Chem. Phys. 22, 2289 (1975).

61. J. S. Turner, Bull. Math. Biology, 22, 205 (1974);
Adv. Chem. Phys. 22, 63 (1975).

62. S. J. Hammarling, " Latent Roots and Latent Vectors"

(University of Toronto Press, Toronto, 1970). Ch. 5.