J. Phys. Chem. 1993,97, 10245-10247
10245
Kinetic Instability in Reaction-Diffusion Systems: The Effect of Time Delay7 S. R. Inamdar and B. D. Kdkami' Chemical Engineering Division, National Chemical Laboratory, Pune-41 I 008, India Received: April 9, 1993; In Final Form: June 25, 1993'
The effects of time delay on the transition curve marking the onset of kinetic instability in a class of reactiondiffusion systems are analytically evaluated by linearizing the system equations. The analysis in general shows that the transition curves become multivalued and indicates considerable shrinkage of the parameter plane where such transitions can occur. Also, the effect is more severe for higher extents of time delay.
Introduction Chemical reactions operating far from equilibrium can show a variety of phenomena that includes multistationarity periodic oscillations, quasiperiodic oscillations, and chaos. The origin of such behavior lies in the nonlinearities present in the system and typically involves autocatalysis, cross catalysis, and other types of feedback loops. In mathematical terms these systems are described by a set of ordinary or partial differential equations, and their solution behavior for variations in parameters has been extensively investigated. In many instances systems may involve time delays in their operative mechanism that consequentlyleads to a set of delay differential equations. In simple terms the rates of change of the variables now depend on the past history and properties of temporal evolution of such systems have led to some useful conclusions and conjectures. A number of earlier investigations demonstrate the utility of differential delay equations in treating the problem of chemical and biochemical systems. They includethe investigation of Rappl showing the occurrence of oscillations in biochemical systems arising due to end product inhibition with a time delay, the onset of respiratory disorders investigated by Mackey and Glass,2and the simulation of DNA replication in T3/T7 bacteriophage in the presence of time delay^.^ The delayed feedback is known4 to stabilize unstable states and can generate new bifurcations in an illuminated thermochemical system. Weiner et aLs studied the experimental chemical oscillator that used a delayed feedback on itself to show that the period of controlled oscillationsincreases or decreases as a function of delay time and also showed the existence of an overlapping region of birhythmicity. Chevalier et a1.6 showed the existence of chaos for the minimal bromate oscillator when nonlinear delayed feedback was employed. Epstein' studied the influence of linear feedback on the reduced Oregonator and found that for an appropriate choice of the delay time, the reduced model behaved much like the full rigorous system. Schell and Rosss provided an exhaustive study detailing the effects of time lags on the temporal evolution of homogeneous chemical reactions and studied the effects of reinjections of the trajectories near the saddle and periodic orbits. They observed the presence of chaos and hyperchaos. More recently, Roesky et aL9 experimentally studied the B-Z reaction with a built-in time delay between the input and output and observed a variety of periodic states when the time delay and the coupling strengths were increased. Chaotic states are produced in the presence of delays which are normally not observed in the free running oscillator for the same residence time. Badola et aI.'O also studied the influence of time delay on autocatalator and showed the existence of bursting solutions. In the present investigation, we consider the effects of time
* Author to whom all correspondence should be addressed. NCL communication no. 5759. .Abstract published in Advance ACS Abstracts, September 1, 1993. f
0022-3654/93/2097-10245$04.00/0
delays on the transition of the threshold of kinetic instability for a reaction-diffusion system. The governing set of equations is linearized making it possible to analytically evaluate the effects of time delays on such transitions. The results are thus limited to the linearized version and not the full nonlinear problem. In a latter section, we shall present the general theory and exemplify it by considering the specific example of the reaction-diffusion system described by Encillator.1' General Theory and Linear Stability Analysis A dynamical system described by a set of first-order ordinary differential equations can be written as x = f(x,L) (1) where x = (xlrx2, ..., xn) is a vector spanned in a n-dimensional space of Rn. The introduction of one or more time delays in the governing system dynamics gives rise to the following system equations -
x = f(X,CL,X1-+vT19Tp. .9Tr) (2) To simplify our task, we shall consider a two-dimensionalsystem and write it in a operator form as
dx
- = kA,x(t-r,) dt
+ N(t,x) +
where the following definitions are made: number of time delays associated with the system dynamics the return value of variable at time ( t + T I ) the constant matrices (or Jacobian matrices) associated with the system dynamics the nonlinear vector remaining after extraction of the linear differential operator but not
r x(t - Ti)
A, N(t,x)
XW- li),t)
containing the time delay variable terms
the nonlinear vector containing only time delay terms after extraction of the linear
operator the time delays present in the system dynamics Following Shimanov,l2 the characteristic equation for a twodimensional system can be written in a matrix form as TI
( i = 1,2,
... r )
A(X) =
For thecaseof thesystem withonly one timedelay, theeigenvalue 0 1993 American Chemical Society
Inamdar and Kulkarni
10246 The Journal of Physical Chemistry, Vol. 97,No. 39, 1993
a2x ax = D,- + xo- Rx + (1 - R)x at a3
- Da, x exp(ay) (12)
- -
In the limiting case of T or T 0, we recover the normal second-order polynomial for the original system without delay. In the general case, expanding the exponential term into a series form, we see that, eq 5 leads to an infinite-order polynomial. The first few terms of this polynomial turn out to be QJ
Equation 6 clearly indicates that, with the introduction of time delay, the system of specified finite dimension becomes a system with infinite dimensions. The equation is suggestive of the fact that the system can exhibit more rich bifurcation cascades under the effect of time delay. We would now like to focus our attention upon the problem of obtaining the transition curve equations in the presence of time delay and for this reason rearrange eq 5 as,
where x and y denote the chemical species, Dal and Da2 the Damkohler numbers, D1 and 0 2 the respective diffusivities; a is the exponential autocatalytic parameter, R the recycle ratio with T units of time delay incurred while reentering the reactor. Equations 12 and 13 provide a description of an exponentially autocatalyzed reaction4iffusion system, and the basic model had been extensivelyanalyzed earlier.” The model schemeshows diverse behavior ranging from a simple stationary point to oscillations and chaotic behavior known to be exhibited by a number of chemical and biochemical systems. The homogeneous solution of the reaction-diffusion system would remain the same even after introduction of time delay. The solution can be given as,
+ + ~ , ( -l 2R)
x0 yo ys =
(2R - 1
+ Da,)
Defining deviations from steady state as x = u + xs,y = u we obtain the evolution equations as
u, = D2u,
+ (Da,eayJ)u - ( a x p a l e a y J- R - Da,)v + (1 - Rblf-,
Defining the first transition bifurcation point as the appearance of a pair of pure imaginary eigenvalues and substituting the eigenvalues as the roots of the characteristic equations leads to two expressions for real and imaginary parts equated to zero. Further simplification gives us
A,
+ A, COS(W7) + A, COS(2OT) + B2w sin(or) - w2 (9)
- A , sin(wr) - A, sin(2w~)+ B,w
+ y,,
+ h ( ~ (17) )
where h(u,u) = aDaleay#(uv). The linearized evolution equations in matrix operator form can be written as D,V - (R + Daleay1) DaleayJ
-axpuleayJ D2V + a x p a l e a y 8- R - Da,
1.
+ B2w COS(WT)= 0 (10)
Defining constants kl, k2, k3, kd, and e as
Defining z = COS(WT), we have
k, = -(R
sin(w7) = (1 - z2)l/’
+ DaleuyJ)
k, = -(axpale&) k, = (Daleaya)
Elimination of the fundamental frequency of oscillation w results in two polynomials of sixth and eighth order, which can be simplified by eliminating z. In fact, such a polynomial reduction can also be effected numerically. The algebraic expression after elimination of w and z = COS(UT) can be used to construct a bifurcation diagram. The Eacillstor Delayed
The reaction-diffusion equations with a single time delay element introduced in the dynamics can be given as
k, = (axpaleaY*- R - Da,) f=(I-R) (19) the characteristic polynomial for the reaction4iffusion system can be written as det([
;:lmzT2
+ k1
k2 -D2m2?r2+ k,
1-
The Journal of Physical Chemistry, Vol. 97, No. 39, I993
Kinetic Instability in Reaction-Diffusion Systems
10247
leading to the eigenvalue problem, -D,mZ?r2
+ k, - ce-"
-X
k, -D,m2?r2
+ k, - te"'
-X
1.
We assume a solution for the linearized equations, giving the eigenfunctions as14
3.50
2.50
Now, considering the transcendental form in eq 5 for the critical transition point defining a boundary for linearly stable solutions, we will obtain the expressions for the constant expressions. They are given as
1
+ k,)(-D,mZ?r2 + k,) - k,k, A, = -e;"(-D,m2?r2 - D2m2?r2+ k, + k,)
A, = (-D,m2a2
1.50
0.00
+
+ k,)
B, = 2cel' The eigenvectors c1 and c2 in eq 22 are obtained as c1 =
(2)'/2 (1
c2 =
+
(2)'/26, (1
+6 y / 2
where
6, =
+
D,mz?r2 ee-"
+ X - k,
k2 We have then obtained the expressions for eigenfunctions and eigenvectors. We now proceed to obtain the expression for the critical wave number me. Equations 9 and 10 are differentiated with respect to T , and the resulting z'equations are equated to obtain a cubic in z as (4A;)z'
+ (4A2A,)z2 + [(2A,(A, - A , ) + A:) + B:w4]z + [A2(A1- A 3 ) + B,B2w4] = 0 (26)
The critical wave number can be obtained by taking a derivative of this expression with respect to mc2and equating it to zero. This results in
+ + D2)?r2(k,D, + k,D,)( 1 - td") +
m~?r44(3ce"'?r2D,D2(D,+ D,)] m:?rZ((D,
2D,D2?r2(k1+ k, - 2 4 2 ) ) + (Ql where the expressions 521 and
522
SI, = 4z3A3, + ee-"(4A3z2)(D,
= (k4D,
(27)
are given by
+ D,)?r2 + 2A,'z +
(2A2ee-")(D, Q2
+ Q,] = 0
+ D,)?r2 + zB:w4
0.02
0.04
0.06
0.08
0.10
Da 1 Figure 1. Bifurcation diagram showing a locus of type I instability points in a delayed reactiondiffusion system for the case of the Encillator model. (1, T = 0.3; 2, T = 0.4; 3, T = 0.5; 4, T = 0.6). Other parameters are a = 20.0, xo = 0.50, R = 0.6, yo = (1 - XO).
A, = (tZe-ZT') B, = m2r2(D, D,) - (k,
4 3 2
Results and Discussion Now that we have obtained most of the properties of the reaction4iffusion system, for the purpose of carrying out numerical experiments,we may construct a bifurcation diagram giving loci of type I instability points. Putting the wave number m to zero, and using eqs 9 and 10 and 14 and 15, we obtain such a bifurcation diagram in the Dal-Da2 plane. Equations 9 and 10lead to a eighth and six order polynomial and require successive polynomial reduction to eliminate z. The results are presented in Figure 1 and clearly show the multivalued nature, that an infinite number of solutions are possible, and that the loops in the diagram keep appearing as we go on reducing the size of the parameter window and reduce the step size and interval for search of solutions to the extent of the machine accuracy level. The bundle of trajectories marked in the figure shows the effect of increasingthe value of the time delay. The transition curve keeps shrinkingas the delay magnitude increasesfor all solutions found. For reasons of comparison, the single transition curve for the case of no time delay, i.e., T = 0, is also presented. In general, introduction of time delay causes the parameter space (in the Dal-Da2 plane) to shrink. The unique transition curve marking the onset of kinetic instability becomes multivalued, in tune with the infinite order that the system assumes in the presence of time delay, and higher values of time delay cause increased shrinkage of the parameter plane where such solutions can be found. References and Notes (1) (2) (3) (4)
Rapp, P. Faraday Symp. Chem. Soc. 1974, 9, 215. Mackey, M. C.; Glass, L. Science 1977, 197, 287. Buchholtz, F.; Schneider, F. W. Biophys. Chem. 1987, 26, 171. Zimmermann, E. C.; Schell, M.; Ross, J. J. Chem. Phys. 1984,81,
.+*.I
IJLI.
(28)
+ k,D2)?r2(2A,z- k, - k,)+ ee-"?r2(D, + D,) X (2A,z + k,k4 - tZe-2T'] + B204a2(D,+ D,) (29)
This method for obtaining the expression for the critical value of wave number m, is analogous to that in the case of plane wave solutions.ls Substituting X = f i w in eqs 27-29 and separating the real and imaginary parts lead to the required expression. Putting the variablesrecycle ratio (R) equal to unity and the time delay ( T ) to zero, for the case of a simple zero eigenvalue, we see that the results match the expressions reported elsewhereI4in the absence of recycle and time delay.
(5) Weiner, J.; Schneider, F. W.; Bar-Eli, K. J . Phys. Chem. 1989, 93, 2704. (6) Chevalier, T.; Freuend, A.; Ross, J. J. Phys. Chem. 1991, 95, 308. (7) Epstein, I. R.; Luo, Y. J . Chem. Phys. 1991, 95, 244. (8) Schell, M.; Ross, J. J . Chem. Phys. 1986, 85, 6489. (9) Roesky, P. W.; Doumbouya, S. I.; Schneider, F. W. J . Phys. Chem. -1993. - - -,97. - . , -198. - -. (10) Badola, P.; Rajani, P.; Ravi Kumar, V.; Kulkarni, B. D. J. Phys. Chem. 1991,95, 2939. ( 1 1 ) Inamdar, S. R.; Rajani, P.; Kulkarni, B. D. J. Phys. Chem. 1991,95, 3422. (12) Shimanov, S. N. J . Appl. Math. Mech. 19!49,23,1198 (Translation of Soviet Journal: PMM 1959, 23 (5), 836). (13) Inamdar, S. R. Ph.D. Thesis, University of Poona, 1990. (14) Inamdar, S. R.; Kulkarni, B. D. J . Phys. A: Math. Gen. 1990,23, L461. (1 5 ) Kuramoto, Y. Chemical Oscillations, Waves, and Turbulence; Springer-Verlag: Berlin, 1984.
