- Email: [email protected]
Homoclinic orbits in a simple chemical model of an autocatalytic reaction
Physica D 32 (1988) 461-470 North-Holland, Amsterdam
HOMOCLINIC
ORBITS IN A SIMPLE CHEMICAL
AUTOCATALYTIC Christian
KAAS-PETERSEN
Department
of Applied Muthemutics
Stephen Department
MODEL
OF AN
REACTION
and Centre for Non-Linear
Studies, University of Leeds, Lee&
LS2 9JT. U.K.
K. SCOTT of Physical Chemistry and Centre for Non-Lmear
Studies, University of Leeds, Leeds LS? UJT, U.K.
Received 5 February 1988 Revised 6 June 1988 Communicated by A. V. Holden
In a one-parameter family of ODES a periodic solution may disappear via the formation of a homoclinic orbit as the parameter c, say, passes through a critical value c*. The period of the solution becomes infinite, growing as In / c ~ c* / as c -+ c* when the stationary solution at which the homoclinic orbit begins and ends is hyperbolic. If the stationary solution is non-hyperbolic, i.e. a bifurcation point, then the period tends to infinity more quickly than this logarithmic form, growing as ( 1c - c’ I)- ‘I*. We show that although bifurcation point and homoclinic orbit formation are both codimension-1 phenomena, a homoclinic orbit related to a non-hyperbolic stationary solution is essentially a codimension-1 situation. Both types of homoclinic orbit have been found in a simple model for an autocatalytic chemical reaction scheme.
1. Introduction
in lRn) exist for c < c*, say, and no closed loops
A homoclinic orbit solution Q(t) for a dynamical system described by an autonomous non-linear ODE f =f( x), is characterized by G(t) + x0 for
exist for c > c* as illustrated in fig. 1. The periodic solution may be stable or unstable - the periodic solution is stable if trace [Df( x0, c*)] < 0, otherwise unstable cf. [ll, p. 2921. We describe here a
f + + 00 and jectory tends
Q(t) + x0 for t + - 00, i.e. the trato the same point x0 in the state
space Iw” as time increases or decreases indefinitely [ll]. For this to occur, x,, must be a stationary solution of the governing ODES, f(xa) = 0. In geometrical terms the trajectory a(t) and the point x0 form a closed loop in the state space. (x0 is necessarily a saddle point and when, as will be the case here, Df(xo, c) has only real eigenvalues, a homoclinic orbit is also known as a saddle connection or separatrix loop [17].) Homoclinic orbits may appear at distinct parameter values of one-parameter families of ODES parametrized by c: 1 =f( x; c). If the homoclinic orbit exists for c = c*, then periodic solutions (also closed loops
0167-2789/88/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
computational bit solutions autonomous,
method for locating homoclinic orin such a one-parameter family of non-linear ODES.
For c < c* there exists a say stable periodic solution with period T. See fig. 1. As c increases towards c*, the period T tends to infinity whilst the amplitude of the solution remains bounded. We show below that T goes to infinity as In ]c - c* 1, when x,,(c*) is a hyperbolic saddle point, cf. [5], but T tends to infinity faster than this logarithmic form when x0( c*) is a non-hyperbolic saddle point (a SNIPER bifurcation [2]). The periodic solution. disappears for c > c*. This transition is sometimes referred to as a homoclinic explosion or a blue-sky catastrophe [22].
B.V.
We present orbit
Q(t),
a method
c*, in the situation linearized
for locating
and the corresponding about
ble manifold, real eigenvalue
a homoclinic
parameter
value
where the vector field f(x; x0 has a one-dimensional
i.e. the matrix
Df(xo;
c)
unsta-
c) has one
X,(c) > 0 and all other eigenvalues
have real parts less than zero, i.e. x0 is a hyperbolic saddle point. The method is based on the observation
that
the homoclinic
out from the stationary unstable one-dimensional manifold, labelled IV” parameter
values
orbit
must
start
solution xg along the manifold. This unstable in fig. 1, exists for all
c close to c*.
The task of finding homoclinic orbits has been approached in other ways. One approach employs path following techniques to track the periodic solution as a function of c until the period becomes sufficiently large [4]; other approaches are based on the formulation of a bifurcation equation via the LyapounovvSchmidt reduction [12]. Homoclinic orbits have been studied extensively in the Lorenz equations [20, 211. The present method should also be adaptable to that scheme, but our motivation has arisen from studies of a dynamical system of chemical reaction equations. The model we investigate involves an autocatalytic reaction in a well-stirred flow reactor (the chemical engineers
2. Theory
Fig. 1. Phase plane diagrams for three situations: (a) c -c c.*a stable periodic solution and corresponding manifolds of stationary solution: (b) c = c*, the two manifolds coincide and a homoclinic orbit forms; (c) c > c*, the periodic solution has disappeared but the manifolds still exist with changed relative positions.
CSTR)
under
isothermal
conditions.
of homoclinic orbits
In this section we demonstrate three features. First, that the period T of a periodic solution approaches infinity logarithmically when the periodic solution approaches a homoclinic orbit whose stationary solution x,, is a hyperbolic saddle point with real eigenvalues. Next we derive the corresponding growth of the period if x0 is a nonhyperbolic saddle point. Then we explain why homoclinic orbits with non-hyperbolic saddle points are codimension-1 events. The key idea is that the periodic solution comes close to a stationary solution x0 such that (i) the system spends most of the time along the trajec-
tory near to x0 and (ii) the vector field near x0 is well approximated by its linearization. We thus split the vector field in a linear part valid near x0 and a nonlinear part that remains almost fixed as the parameter c is varied. Similar reasoning was used in [21]. Let
i=f(x;c),
XER”,
CER
be a system of ODES which has a stationary solution x0 = x0(c). We assume x0 is a saddle point and, in fact, further that Df(x,; c) has only real eigenvalues which are ordered as X ,I I X, ~, I . . I X, < 0 < hi. Associated with these eigenvalues we have n eigenvectors pi,. , u,,. Near x0 the vector field is well represented by its linear approximation Df(x,; c)(x - xc,). If x belongs to the span { 02,. . . , q,} it will approach x0 from the uZ direction, as X, has the smallest magnitude. With the above ordering of the eigenvalues all other directions become unimportant and hence we concentrate on the system with n = 2. To be consistent with fig. 1 we assume hi > -X, so that the periodic solution is stable. Let us also assume that co-ordinates are chosen such that x,) is at the origin and ui and v2 are the axes. The linearized vector field can then be written 2.2,= A,u,, k, = h,u,.
C. Kaas-Petersen
and S. K. Scott/Homo&Cc
463
orbits
ua
:11_I: -_- ---
I
c-s.011
1
7
t _.__..-_!
-..
/
I
Fig. 2. Trajectory saddle point.
(s,O)
Ui
“2
of linear system in the neighbourhood
of a
Consider fig. 2, where a box B is shown through which passes a trajectory of the linearized flow. The trajectory starting at the point (s, E) on one side of B at time t = 0 emerges from B at the point (77, s) at time 7:
Thus the time taken to travel through the box is r = (l/X,)ln(.s/~). This describes the flow near to the stationary solution x0. Now we return to the periodic solution. Let c < c* (see fig. l), such that the dynamical system has a periodic solution coexisting with the stationary solution. (There must exist another stationary solution inside the closed loop.) When we are near to homoclinic behaviour (c close to c*) the trajectory spends a lot of time near to x0. The period of the solution is T which we will express in the form T = a + 7, where r is the time spent moving through a small box near to x0 and (Y is the time spent moving along the rest of the loop (i.e. to return to the box). Substituting from above for r we obtain T = a + (l/A,)
cS,O)
Fig. 3. Homoclinic orbit attached to non-hyperbolic saddle point for c = c*. For c < c* there is no stationary solution in the box B; for c > c* there are two stationary solutions.
I-
4
IB
In (S/a),
where s is the box size and is a number fixed by us. We wish to relate E to the parameter c: E is the distance from the periodic solution to the stable manifold Ws measured along the line M. At c = c* we have E = 0: for c < c*, E > 0. Let us introduce
the distance d, as the distance from W” is positive for c < along M. Then d, d, < E. Returning again to c = c* we have = 0. For c > c*, E is no longer defined as no periodic solution: d, is, however, still and is negative. We assume that then typically d, = dh( c*)( c terms, with d&(c*) # 0, [ll, p. E is approximated by a similar only for c 5 c*). Thus we obtain T=(Y+filnIc-c*I
to W” c* and d,(c*)
there is defined
we may linearize, c*) + higher order 2921. The distance expression (valid
asc-*c*,
where /? involves A,, s and dh(c*). A similar formula has been derived for the case in which X, and A, are complex conjugate eigenvalues with negative real parts [5]. When the stationary solution is a non-hyperbolic saddle point, as shown in fig. 3, the flow through the box B will have the form til=U:-(c-c*), l.i,=
-Au,.
For c < c* there are no stationary solutions and the time taken to traverse the box B from the left side to the right is found by solving the ur equation. The period of the orbit is then dominated by the time for this traverse (plus a constant for the period along the remainder of the orbit loop). We get T=
To+
2 JIc _ c*l tan-’
&c*,
as c + c*.
C. Km.-Petersen
464
and S. K. Scott / Homochc
orbits
\
x2 '
y(t)
\ \ "I
y(t') y(0)
Fig. 4. The unstable manifolds W” on the outside of the two stable manifolds WIS and W”’ allows homoclinic orbits attached to a bifurcation point.
x0'
i As tan-’ E converges to 7r/2 for E + 0, the period grows like Jc - c* 1-I/*, which is faster than logarithmic. This behaviour is confirmed by numerical results. The occurrence of a homoclinic orbit with a non-hyperbolic saddle point is a codimension-1 phenomenon. The explanation for this is as follows (see fig. 4). With c > c* there are no periodic solutions, and the trajectory W” leaving the unstable stationary point is attracted to the stable stationary solution. In its approach, W” passes on the outside of the stable manifold W,S. The relative orientations of the stable and unstable manifolds shown in this picture will remain the same for all values of c up to c*, where the saddle node bifurcation point is found, and there the homoclinic orbit is formed. For c < c* the saddle point has disappeared and we are left with the limit cycle. The important feature here is that the trajectory leaving the unstable saddle point is sent round into the region of attraction of the stable stationary
solution.
Fig. 5. $(I) is fold emanating imation of the before shooting
;P
\
method
In this section we describe a method mine homoclinic orbits of the non-linear cal system of ODES
to deterdynami-
XER”,
A, < 0 -c A,, the stationary state is a saddle point. Corresponding to the eigenvalue h, there exists an eigenvector q. The line (one dimensional manifold) x0 + hu, is a linear approximation to the unstable manifold W”. Let h, be a small number and let us obtain the trajectory #(t) of the ODES starting at x0 + h,v, at time t = 0. If h, is small enough, 4 (t) will be a good approximation to the unstable manifold W “. After some time t*, +(t) will again approach the stationary solution x0 (see fig. 5). The distance d = \I$( t) - x01( becomes smaller until d reaches a minimum at a certain to the closest approach. time t, corresponding After this time d will grow again: then either G(t) opposite
to q
(shown
in fig. 5)
so that the velocity vector f( $ (t), c) and q are almost antiparallel or 4(t) moves in the q direction. In the first case let the function G take the value - 1: in the second case let G = + 1. The function G is parametrized by c with
CEK!. G=G(c)=
Let (x,; c) be a stationary solution. The stability of this solution is found by examining the eigen-
x1
values X, of the n x n matrix Df(x,; c). Suppose the eigenvalues are such that X, I X, I < . . I
+1
i=f(x;c),
/ \
a trajectory approximating the unstable manifrom xg. The vector o, spans the linear approxmanifold. $(I) comes closest to x0 at time t,, off and away from G(O).
takes a direction 3. Computational
Yhrl'
I
- 1
\
if # ( t ) moves coparallel with ur on return, if q(t)
moves antiparallel
with ur on return.
465
C. Kaas-Petersen and SK. Scott/ Homoclinic orhris
We may now use the method
of false position
to
bracket the parameter value of c* for which the homoclinic orbit occurs and at which G changes
coefficients
for all species
metrical
profile,
become
[3, 181
and an assumed
and the dimensionless
sym-
equations
sign. Here we shall just give some practical the
computations.
solution
near
the direction
First
we compute
to the homoclinic
tive error control of lo-* and have taken h,= 10e6. The actual value of h, effects the time taken for the integration, i.e. t,, but has relatively little influence on the result for c*. We did not compute the minimal distance d accurately, but simply obtained the points of q(t) equidistant in time, with unit time step between successive points. This minimal distance could however be combined with the appropriate sign in defining G.
4. Model of an autocatalytic chemical system isothermal, A + B + C.
The first step A + B is taken to be an autocatalytic process. Reaction proceeds in a well-stirred tank into which a fresh supply of the reactants A and B flows continuously and from which there is a matching outflow of a mixture of A, B and C. The concentrations of the two independent species A and B are governed by a set of mass-balance equations which are non-linear ODES. In a suitable dimensionless form these may be written as [3,7,81 da/dt
=
t,(l
-
a)
dfi/dt=t,(&,-/3)+(up2+
-
aP2
-
KJY,
This gives
are eigenvec-
tors so care must be taken when normalizing ol. We have solved our particular ODES using a rela-
We consider here the prototype non-linear chemical reaction scheme
i&q’& = D a2a/ax2 - LX/~’-
a periodic
orbit.
of ui: both oi and -vi
details of
K”(Y,
K,(Y-K2P.
If we consider the same reaction scheme but with no stirring and the supply and removal replaced by diffusional processes across boundaries then we obtain a set of reaction-diffusion equations. In the simplest form we may take equal diffusion
with
aqax=ap/ax=o
at
(~=l
atx=l.
and
p=&
x=0,
The stationary state and Hopf bifurcation responses of this model have been discussed extensively elsewhere [15, 181. These two systems admit homoclinic orbits. In both systems homoclinic orbits have been located using the method of the previous section. For a typical example we specify the values of K",~~ and &. The dimensionless residence time t, (or diffusion coefficient D) plays the role of the parameter c, and the critical value t,* is determined accurately. In this case the stable limit cycle oscillations exist for a range of t,? -c t, and these are then computed. We have plotted the period of these oscillations against In )t, - t:) and a satisfactory linear relationship is observed as predicted. In fig. 6(a) we plot the stationary solution cy,, against t,. The stationary locus has two folds and hence two regions of multiple solutions. The pattern displayed is termed a mushroom (this is more obvious when 1 - (Y,, is used as the ordinate). In between the regions of multiple solutions there is a Hopf bifurcation at t: along one branch from which emerges a family of stable periodic orbits. The amplitude of the limit cycle increases as t, decreases below t,“. The periodic solution clearly terminates via a homoclinic orbit at the point t; as the limit cycle “collides” with the middle branch of saddle points, in the way described above. If one of the other parameters, & say, is now changed slightly there is a quantitative change in the diagram, see fig. 6(b). The point along the middle branch at which the homoclinic orbit is
formed moves closer to the upper non-hyperbolic saddle point which occurs at fy. For some particular &,, see fig. 6(c), the homoclinic orbit tion coincides exactly with the bifurcation l* = tr. find
that
bifurcation
We note that if we continue the
homoclinic
point
persists,
orbit
%s
formapoint
to vary & we
attached
i.e. for fixed
K~
to the and
K,
there is a range of & over which this phase-space arrangement exists.
5. Path following for homoclinic orbits lb) By plotting the values of t, at which turning points (ty) or Hopf bifurcations (ty) or homoclinic orbits (t:) for
0.05 and K, = (The plane has been extended to negative of the dimensionless concentration ratio &: unrealistic physically but it is convenient to this region mathematically.) Points A and degenerate bifurcation points where two K* =
0.0005. values this is include F are
eigenvalues are equal to zero. The locus of homoclinic orbits emerges from the degenerate point A (where the amplitude of the orbit is zero) and moves smoothly across the parameter plane until it reaches the curve of bifurcation point at the point locus
B (shown in greater detail in fig. 7(b)). The of homoclinic orbits does not cease at this
point but follows the curve from B to F. The stability of the periodic solution bifurcating along the Hopf curve (t:) has been evaluated routinely, using BIFOR2 [13]. Between points A and D and from E to F the periodic solution is unstable, but otherwise it is stable. The significance of this is that near D the homoclinic locus lies to the right (higher t,) of the Hopf point and hence we can infer the existence of a turning point or saddle-node bifurcation in the path of periodic solutions. Fig. 8 shows some accurate examples of stationary and periodic solutions loci for systems with PO varying from -0.01 to 0.2. These show multiple stationary states, multiple Hopf bifurca-
(cl
Fig. 6. Schematic representation of stationary and periodic solutions as function of 1, for various /IO identifying parameter values for Hopf bifurcation t,“, saddle-node bifurcation r,‘ and homoclinic orbit I:,
tions, multiplicity of periodic solutions and homoclinic orbits formed from both stable and unstable limit cycles.
6. Discussion and conclusions We have presented a new method for locating homoclinic orbits which can be applied to sets of non-linear ODES and PDEs. This approach differs from previous algorithms. For instance in [l, 41 periodic solutions are analyzed as fixed point of a Poincare map. Specifically, in [4] one finds a periodic solution with a specified (large) period and in
C. Km.-Perersen
nnd S.K. Scott/
467
Homoclinic orbits
1
-0
0
1
2
tr
i 0
3
Fig. 7. Bifurcation diagram in fr-& parameter plane for K, = 0.0005: (a) the closed curve with two cusps, AGF solutions. The curve ADEF (H is not on this curve) of elsewhere unstable. The locus A B F of homoclinic orbits.
02
-61
tr
03
the simple autocatalytic reaction in the well-stirred case with K> = 0.05 and HA of bifurcation points, within this region there are three stationary Hopf bifurcation points: on section DE the periodic solutions are stable, (b) Blow-up of region close to point B.
[l], which can only be applied to systems of two coupled equations, the return plane in the Poincare map is taken as the line dividing the angle spanned by the two eigenvectors ur and 0,. Periodic solutions only exist on one side of the homoclinic orbit and so this latter method does not permit the efficient bracketing possible with our method. We have not developed the path following routine for tracing the homoclinic curve in twoparameter systems, but this should be possible as should extension of the method to the location of heteroclinic orbits. Our systems have been dissipative rather than conservative. For Hamiltonian systems perturbed by a small time-dependent periodic forcing term the Melnikov technique is applicable. Our method is essentially based on a one-dimensional unstable manifold. For systems such as the Lorenz equations the trajectories shoot
in and
spiral
out
near
the non-zero
stationary
states. To examine this we either need to modify the method here or to reverse time in the equations. Evidence
for homoclinic
orbits
at non-hyper-
bolic saddle points in real chemical systems has been reported by Gray et al. [9, lo] in the hydrogen-oxygen and carbon monoxide-oxygen reactions. Both of these systems display an explosion limit in well-stirred flow reactors separating steady reaction (low reaction rate) from ignition (high stationary state reaction rate). On crossing this curve of bifurcation points (by raising the reaction temperature) the system jumps into large amplitude oscillations. Further increase in ambient temperature reveals a classic Hopf bifurcation extinction of the limit cycle, which reappears as the temperature is reduced again, without hysteresis.
C. Kaa.-Petersen
and S.K. Scott / Homochic
no=02
10
orbits
1 0 -1
1
no=0
4
0.5
OOC~
0
-
-T
01
02
+,
03
OOL-0
I
02
01
03
tr
(b)
-.. f-R no=005
,/I
1.0
00’0 02
.
I
d
_*--
/
---
,/-
/
/’
^
:
1
I I \ \
00
I--
O2
01
t,
n
.03
Cc) Fig. 8. Accurately computed stationary and periodic solution diagrams for the system various p,,: full lines - stable stationary solutions; dashed lines unstable stationary periodic solutions; open circles- amplitude of unstable periodic solutions.
01
02
t,
03
Cd) in fig. 5 with K~ = 0.05 and K,, = 0.0005 for solutions: full circles - amplitude of stable
C. Kaar-Petersen
and S. K. Scott/
Homoclimc
00 .Ol
02
tr
03
orbits
469
01
02
t,
03
(0
(e) Fig.
8. Continued
Further decrease in temperature causes a lengthening in period, with an apparent faster than logarithmic growth. The large amplitude oscillations are extinguished at exactly the same ambient temperature as that at which they appear on the upward traverse, so there is no hysteresis which would normally be associated with homoclinic orbit formation. Acknowledgements We thank Dr. C. Sparrow for providing us with helpful insight into the origin of homoclinic orbits with saddle points. We are grateful to Professors J.K. Hale and R. Aris for suggesting this problem and to SERC for support (grant GR/D 56990).
References [l] R. At-is and I.G. Kevrekidis, private communication (1987). [2] K. Bar-Eli and R.M. Noyes, Relevance of a two-variable
Oregonator to stable and unstable steady states and limit cycles, to thresholds of excitability and to Hopf vs. SNIPER bifurcations, J. Chem. Phys. 86 (1987) 1927-1937. [3] J. Brindley. C. Kaas-Petersen. J.H. Merkin and SK. Scott, in preparation. [4] E. Doedel, AUTO: Software for continuation and bifurcation problems in ordinary differential equations, Caltech (1986). [5] P. Glendinmng and C. Sparrow, Local and global behaviour near homoclinic orbits, J. Stat. Phys. 35 (1984) 645-697. [6] B.F. Gray and S.K. Scott, Multistability and sustained oscillations in isothermal open systems, J. Chem. Sot. Faraday Trans. 1, 81 (1985) 1563-1567. [7] P. Gray and S.K. Scott, Autocatalytic reactions in the isothermal continuous stirred-tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38 (1983) 29-43. [8] P. Gray and SK. Scott, Autocatalytic reactions in the isothermal continuous stirred-tank reactor: oscillations and instabilities, Chem. Eng. Sci. 39 (1984) 1087-1097. [9] P. Gray, J.F. Griffiths and SK. Scott, Proc. R. Sot. London A 394 (1984) 243-258. [lo] P. Gray, J.F. Griffiths and SK. Scott, Proc. R. Sot. London A397 (1985) 21-44. [Ill J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields (Springer, New York, Berlin, Heidelberg, 1983).
470
C. Km.-Petersen
und S.K. Scott/
[12] J.K. Hale, Introduction to dynamic bifurcation, in: L. Salvadori, ed. Bifurcation Theory and Applications, Springer Lecture Notes in Mathematics, no. 1057 (1984). [13] B.D. Hassard, N.D. Kazarinoff and Y.-H. Wan, Theory and Applications of Hopf Bifurcation (Cambridge Univ. Press, 1981). [14] C. Kaas-Petersen, PATH - User’s guide, Centre for Nonlinear Studies, University of Leeds (1987). [15] C. Kaas-Petersen and S.K. Scott, Stationary-state and Hopf bifurcation patterns in isothermal reaction-diffusion systems, Chem. Eng. Sci. 43 (1988) 391-392. [16] S.R. Kay, SK. Scott and P.G. Lignola, The application of singularity theory to isothermal autocatalytic reactions, Proc. R. Sot. London A409 (1987) 433-448. [17] S. Schecter, The saddle-node separatrix-loop bifurcation,
Homoclinic orhts
SIAM J. Math. Anal. 18 (19X7) 1142-1156. [18] S.K. Scott, Isolas, mushrooms and oscillations In isothcrmal autocatalytic reaction-dilfusion equations, Chem. Eng. Sci. 42 (1987) 307-316. [19] SK. Scott and W.W. Farr, Dynamic tine structure in the cubic autocatalator, submitted to Chem. Eng. Sci. [20] C. Sparrow, The Lorenz Equations: Bifurcations, Chaos and Strange Attractors (Springer, New York, Heidelberg. Berlin, 1982). [21] C. Sparrow, The Lorem equations, in: A. Holden, ed. Chaos (Manchester Univ. Press, 19X6). [22] H.B. Stewart and J.M.T. Thompson, Towards a clasaitication of generic bifurcations in dissipative dynamical systems, Dynamics and Stability of Systems 1 (19X6) X7-96.
HOMOCLINIC
ORBITS IN A SIMPLE CHEMICAL
AUTOCATALYTIC Christian
KAAS-PETERSEN
Department
of Applied Muthemutics
Stephen Department
MODEL
OF AN
REACTION
and Centre for Non-Linear
Studies, University of Leeds, Lee&
LS2 9JT. U.K.
K. SCOTT of Physical Chemistry and Centre for Non-Lmear
Studies, University of Leeds, Leeds LS? UJT, U.K.
Received 5 February 1988 Revised 6 June 1988 Communicated by A. V. Holden
In a one-parameter family of ODES a periodic solution may disappear via the formation of a homoclinic orbit as the parameter c, say, passes through a critical value c*. The period of the solution becomes infinite, growing as In / c ~ c* / as c -+ c* when the stationary solution at which the homoclinic orbit begins and ends is hyperbolic. If the stationary solution is non-hyperbolic, i.e. a bifurcation point, then the period tends to infinity more quickly than this logarithmic form, growing as ( 1c - c’ I)- ‘I*. We show that although bifurcation point and homoclinic orbit formation are both codimension-1 phenomena, a homoclinic orbit related to a non-hyperbolic stationary solution is essentially a codimension-1 situation. Both types of homoclinic orbit have been found in a simple model for an autocatalytic chemical reaction scheme.
1. Introduction
in lRn) exist for c < c*, say, and no closed loops
A homoclinic orbit solution Q(t) for a dynamical system described by an autonomous non-linear ODE f =f( x), is characterized by G(t) + x0 for
exist for c > c* as illustrated in fig. 1. The periodic solution may be stable or unstable - the periodic solution is stable if trace [Df( x0, c*)] < 0, otherwise unstable cf. [ll, p. 2921. We describe here a
f + + 00 and jectory tends
Q(t) + x0 for t + - 00, i.e. the trato the same point x0 in the state
space Iw” as time increases or decreases indefinitely [ll]. For this to occur, x,, must be a stationary solution of the governing ODES, f(xa) = 0. In geometrical terms the trajectory a(t) and the point x0 form a closed loop in the state space. (x0 is necessarily a saddle point and when, as will be the case here, Df(xo, c) has only real eigenvalues, a homoclinic orbit is also known as a saddle connection or separatrix loop [17].) Homoclinic orbits may appear at distinct parameter values of one-parameter families of ODES parametrized by c: 1 =f( x; c). If the homoclinic orbit exists for c = c*, then periodic solutions (also closed loops
0167-2789/88/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
computational bit solutions autonomous,
method for locating homoclinic orin such a one-parameter family of non-linear ODES.
For c < c* there exists a say stable periodic solution with period T. See fig. 1. As c increases towards c*, the period T tends to infinity whilst the amplitude of the solution remains bounded. We show below that T goes to infinity as In ]c - c* 1, when x,,(c*) is a hyperbolic saddle point, cf. [5], but T tends to infinity faster than this logarithmic form when x0( c*) is a non-hyperbolic saddle point (a SNIPER bifurcation [2]). The periodic solution. disappears for c > c*. This transition is sometimes referred to as a homoclinic explosion or a blue-sky catastrophe [22].
B.V.
We present orbit
Q(t),
a method
c*, in the situation linearized
for locating
and the corresponding about
ble manifold, real eigenvalue
a homoclinic
parameter
value
where the vector field f(x; x0 has a one-dimensional
i.e. the matrix
Df(xo;
c)
unsta-
c) has one
X,(c) > 0 and all other eigenvalues
have real parts less than zero, i.e. x0 is a hyperbolic saddle point. The method is based on the observation
that
the homoclinic
out from the stationary unstable one-dimensional manifold, labelled IV” parameter
values
orbit
must
start
solution xg along the manifold. This unstable in fig. 1, exists for all
c close to c*.
The task of finding homoclinic orbits has been approached in other ways. One approach employs path following techniques to track the periodic solution as a function of c until the period becomes sufficiently large [4]; other approaches are based on the formulation of a bifurcation equation via the LyapounovvSchmidt reduction [12]. Homoclinic orbits have been studied extensively in the Lorenz equations [20, 211. The present method should also be adaptable to that scheme, but our motivation has arisen from studies of a dynamical system of chemical reaction equations. The model we investigate involves an autocatalytic reaction in a well-stirred flow reactor (the chemical engineers
2. Theory
Fig. 1. Phase plane diagrams for three situations: (a) c -c c.*a stable periodic solution and corresponding manifolds of stationary solution: (b) c = c*, the two manifolds coincide and a homoclinic orbit forms; (c) c > c*, the periodic solution has disappeared but the manifolds still exist with changed relative positions.
CSTR)
under
isothermal
conditions.
of homoclinic orbits
In this section we demonstrate three features. First, that the period T of a periodic solution approaches infinity logarithmically when the periodic solution approaches a homoclinic orbit whose stationary solution x,, is a hyperbolic saddle point with real eigenvalues. Next we derive the corresponding growth of the period if x0 is a nonhyperbolic saddle point. Then we explain why homoclinic orbits with non-hyperbolic saddle points are codimension-1 events. The key idea is that the periodic solution comes close to a stationary solution x0 such that (i) the system spends most of the time along the trajec-
tory near to x0 and (ii) the vector field near x0 is well approximated by its linearization. We thus split the vector field in a linear part valid near x0 and a nonlinear part that remains almost fixed as the parameter c is varied. Similar reasoning was used in [21]. Let
i=f(x;c),
XER”,
CER
be a system of ODES which has a stationary solution x0 = x0(c). We assume x0 is a saddle point and, in fact, further that Df(x,; c) has only real eigenvalues which are ordered as X ,I I X, ~, I . . I X, < 0 < hi. Associated with these eigenvalues we have n eigenvectors pi,. , u,,. Near x0 the vector field is well represented by its linear approximation Df(x,; c)(x - xc,). If x belongs to the span { 02,. . . , q,} it will approach x0 from the uZ direction, as X, has the smallest magnitude. With the above ordering of the eigenvalues all other directions become unimportant and hence we concentrate on the system with n = 2. To be consistent with fig. 1 we assume hi > -X, so that the periodic solution is stable. Let us also assume that co-ordinates are chosen such that x,) is at the origin and ui and v2 are the axes. The linearized vector field can then be written 2.2,= A,u,, k, = h,u,.
C. Kaas-Petersen
and S. K. Scott/Homo&Cc
463
orbits
ua
:11_I: -_- ---
I
c-s.011
1
7
t _.__..-_!
-..
/
I
Fig. 2. Trajectory saddle point.
(s,O)
Ui
“2
of linear system in the neighbourhood
of a
Consider fig. 2, where a box B is shown through which passes a trajectory of the linearized flow. The trajectory starting at the point (s, E) on one side of B at time t = 0 emerges from B at the point (77, s) at time 7:
Thus the time taken to travel through the box is r = (l/X,)ln(.s/~). This describes the flow near to the stationary solution x0. Now we return to the periodic solution. Let c < c* (see fig. l), such that the dynamical system has a periodic solution coexisting with the stationary solution. (There must exist another stationary solution inside the closed loop.) When we are near to homoclinic behaviour (c close to c*) the trajectory spends a lot of time near to x0. The period of the solution is T which we will express in the form T = a + 7, where r is the time spent moving through a small box near to x0 and (Y is the time spent moving along the rest of the loop (i.e. to return to the box). Substituting from above for r we obtain T = a + (l/A,)
cS,O)
Fig. 3. Homoclinic orbit attached to non-hyperbolic saddle point for c = c*. For c < c* there is no stationary solution in the box B; for c > c* there are two stationary solutions.
I-
4
IB
In (S/a),
where s is the box size and is a number fixed by us. We wish to relate E to the parameter c: E is the distance from the periodic solution to the stable manifold Ws measured along the line M. At c = c* we have E = 0: for c < c*, E > 0. Let us introduce
the distance d, as the distance from W” is positive for c < along M. Then d, d, < E. Returning again to c = c* we have = 0. For c > c*, E is no longer defined as no periodic solution: d, is, however, still and is negative. We assume that then typically d, = dh( c*)( c terms, with d&(c*) # 0, [ll, p. E is approximated by a similar only for c 5 c*). Thus we obtain T=(Y+filnIc-c*I
to W” c* and d,(c*)
there is defined
we may linearize, c*) + higher order 2921. The distance expression (valid
asc-*c*,
where /? involves A,, s and dh(c*). A similar formula has been derived for the case in which X, and A, are complex conjugate eigenvalues with negative real parts [5]. When the stationary solution is a non-hyperbolic saddle point, as shown in fig. 3, the flow through the box B will have the form til=U:-(c-c*), l.i,=
-Au,.
For c < c* there are no stationary solutions and the time taken to traverse the box B from the left side to the right is found by solving the ur equation. The period of the orbit is then dominated by the time for this traverse (plus a constant for the period along the remainder of the orbit loop). We get T=
To+
2 JIc _ c*l tan-’
&c*,
as c + c*.
C. Km.-Petersen
464
and S. K. Scott / Homochc
orbits
\
x2 '
y(t)
\ \ "I
y(t') y(0)
Fig. 4. The unstable manifolds W” on the outside of the two stable manifolds WIS and W”’ allows homoclinic orbits attached to a bifurcation point.
x0'
i As tan-’ E converges to 7r/2 for E + 0, the period grows like Jc - c* 1-I/*, which is faster than logarithmic. This behaviour is confirmed by numerical results. The occurrence of a homoclinic orbit with a non-hyperbolic saddle point is a codimension-1 phenomenon. The explanation for this is as follows (see fig. 4). With c > c* there are no periodic solutions, and the trajectory W” leaving the unstable stationary point is attracted to the stable stationary solution. In its approach, W” passes on the outside of the stable manifold W,S. The relative orientations of the stable and unstable manifolds shown in this picture will remain the same for all values of c up to c*, where the saddle node bifurcation point is found, and there the homoclinic orbit is formed. For c < c* the saddle point has disappeared and we are left with the limit cycle. The important feature here is that the trajectory leaving the unstable saddle point is sent round into the region of attraction of the stable stationary
solution.
Fig. 5. $(I) is fold emanating imation of the before shooting
;P
\
method
In this section we describe a method mine homoclinic orbits of the non-linear cal system of ODES
to deterdynami-
XER”,
A, < 0 -c A,, the stationary state is a saddle point. Corresponding to the eigenvalue h, there exists an eigenvector q. The line (one dimensional manifold) x0 + hu, is a linear approximation to the unstable manifold W”. Let h, be a small number and let us obtain the trajectory #(t) of the ODES starting at x0 + h,v, at time t = 0. If h, is small enough, 4 (t) will be a good approximation to the unstable manifold W “. After some time t*, +(t) will again approach the stationary solution x0 (see fig. 5). The distance d = \I$( t) - x01( becomes smaller until d reaches a minimum at a certain to the closest approach. time t, corresponding After this time d will grow again: then either G(t) opposite
to q
(shown
in fig. 5)
so that the velocity vector f( $ (t), c) and q are almost antiparallel or 4(t) moves in the q direction. In the first case let the function G take the value - 1: in the second case let G = + 1. The function G is parametrized by c with
CEK!. G=G(c)=
Let (x,; c) be a stationary solution. The stability of this solution is found by examining the eigen-
x1
values X, of the n x n matrix Df(x,; c). Suppose the eigenvalues are such that X, I X, I < . . I
+1
i=f(x;c),
/ \
a trajectory approximating the unstable manifrom xg. The vector o, spans the linear approxmanifold. $(I) comes closest to x0 at time t,, off and away from G(O).
takes a direction 3. Computational
Yhrl'
I
- 1
\
if # ( t ) moves coparallel with ur on return, if q(t)
moves antiparallel
with ur on return.
465
C. Kaas-Petersen and SK. Scott/ Homoclinic orhris
We may now use the method
of false position
to
bracket the parameter value of c* for which the homoclinic orbit occurs and at which G changes
coefficients
for all species
metrical
profile,
become
[3, 181
and an assumed
and the dimensionless
sym-
equations
sign. Here we shall just give some practical the
computations.
solution
near
the direction
First
we compute
to the homoclinic
tive error control of lo-* and have taken h,= 10e6. The actual value of h, effects the time taken for the integration, i.e. t,, but has relatively little influence on the result for c*. We did not compute the minimal distance d accurately, but simply obtained the points of q(t) equidistant in time, with unit time step between successive points. This minimal distance could however be combined with the appropriate sign in defining G.
4. Model of an autocatalytic chemical system isothermal, A + B + C.
The first step A + B is taken to be an autocatalytic process. Reaction proceeds in a well-stirred tank into which a fresh supply of the reactants A and B flows continuously and from which there is a matching outflow of a mixture of A, B and C. The concentrations of the two independent species A and B are governed by a set of mass-balance equations which are non-linear ODES. In a suitable dimensionless form these may be written as [3,7,81 da/dt
=
t,(l
-
a)
dfi/dt=t,(&,-/3)+(up2+
-
aP2
-
KJY,
This gives
are eigenvec-
tors so care must be taken when normalizing ol. We have solved our particular ODES using a rela-
We consider here the prototype non-linear chemical reaction scheme
i&q’& = D a2a/ax2 - LX/~’-
a periodic
orbit.
of ui: both oi and -vi
details of
K”(Y,
K,(Y-K2P.
If we consider the same reaction scheme but with no stirring and the supply and removal replaced by diffusional processes across boundaries then we obtain a set of reaction-diffusion equations. In the simplest form we may take equal diffusion
with
aqax=ap/ax=o
at
(~=l
atx=l.
and
p=&
x=0,
The stationary state and Hopf bifurcation responses of this model have been discussed extensively elsewhere [15, 181. These two systems admit homoclinic orbits. In both systems homoclinic orbits have been located using the method of the previous section. For a typical example we specify the values of K",~~ and &. The dimensionless residence time t, (or diffusion coefficient D) plays the role of the parameter c, and the critical value t,* is determined accurately. In this case the stable limit cycle oscillations exist for a range of t,? -c t, and these are then computed. We have plotted the period of these oscillations against In )t, - t:) and a satisfactory linear relationship is observed as predicted. In fig. 6(a) we plot the stationary solution cy,, against t,. The stationary locus has two folds and hence two regions of multiple solutions. The pattern displayed is termed a mushroom (this is more obvious when 1 - (Y,, is used as the ordinate). In between the regions of multiple solutions there is a Hopf bifurcation at t: along one branch from which emerges a family of stable periodic orbits. The amplitude of the limit cycle increases as t, decreases below t,“. The periodic solution clearly terminates via a homoclinic orbit at the point t; as the limit cycle “collides” with the middle branch of saddle points, in the way described above. If one of the other parameters, & say, is now changed slightly there is a quantitative change in the diagram, see fig. 6(b). The point along the middle branch at which the homoclinic orbit is
formed moves closer to the upper non-hyperbolic saddle point which occurs at fy. For some particular &,, see fig. 6(c), the homoclinic orbit tion coincides exactly with the bifurcation l* = tr. find
that
bifurcation
We note that if we continue the
homoclinic
point
persists,
orbit
%s
formapoint
to vary & we
attached
i.e. for fixed
K~
to the and
K,
there is a range of & over which this phase-space arrangement exists.
5. Path following for homoclinic orbits lb) By plotting the values of t, at which turning points (ty) or Hopf bifurcations (ty) or homoclinic orbits (t:) for
0.05 and K, = (The plane has been extended to negative of the dimensionless concentration ratio &: unrealistic physically but it is convenient to this region mathematically.) Points A and degenerate bifurcation points where two K* =
0.0005. values this is include F are
eigenvalues are equal to zero. The locus of homoclinic orbits emerges from the degenerate point A (where the amplitude of the orbit is zero) and moves smoothly across the parameter plane until it reaches the curve of bifurcation point at the point locus
B (shown in greater detail in fig. 7(b)). The of homoclinic orbits does not cease at this
point but follows the curve from B to F. The stability of the periodic solution bifurcating along the Hopf curve (t:) has been evaluated routinely, using BIFOR2 [13]. Between points A and D and from E to F the periodic solution is unstable, but otherwise it is stable. The significance of this is that near D the homoclinic locus lies to the right (higher t,) of the Hopf point and hence we can infer the existence of a turning point or saddle-node bifurcation in the path of periodic solutions. Fig. 8 shows some accurate examples of stationary and periodic solutions loci for systems with PO varying from -0.01 to 0.2. These show multiple stationary states, multiple Hopf bifurca-
(cl
Fig. 6. Schematic representation of stationary and periodic solutions as function of 1, for various /IO identifying parameter values for Hopf bifurcation t,“, saddle-node bifurcation r,‘ and homoclinic orbit I:,
tions, multiplicity of periodic solutions and homoclinic orbits formed from both stable and unstable limit cycles.
6. Discussion and conclusions We have presented a new method for locating homoclinic orbits which can be applied to sets of non-linear ODES and PDEs. This approach differs from previous algorithms. For instance in [l, 41 periodic solutions are analyzed as fixed point of a Poincare map. Specifically, in [4] one finds a periodic solution with a specified (large) period and in
C. Km.-Perersen
nnd S.K. Scott/
467
Homoclinic orbits
1
-0
0
1
2
tr
i 0
3
Fig. 7. Bifurcation diagram in fr-& parameter plane for K, = 0.0005: (a) the closed curve with two cusps, AGF solutions. The curve ADEF (H is not on this curve) of elsewhere unstable. The locus A B F of homoclinic orbits.
02
-61
tr
03
the simple autocatalytic reaction in the well-stirred case with K> = 0.05 and HA of bifurcation points, within this region there are three stationary Hopf bifurcation points: on section DE the periodic solutions are stable, (b) Blow-up of region close to point B.
[l], which can only be applied to systems of two coupled equations, the return plane in the Poincare map is taken as the line dividing the angle spanned by the two eigenvectors ur and 0,. Periodic solutions only exist on one side of the homoclinic orbit and so this latter method does not permit the efficient bracketing possible with our method. We have not developed the path following routine for tracing the homoclinic curve in twoparameter systems, but this should be possible as should extension of the method to the location of heteroclinic orbits. Our systems have been dissipative rather than conservative. For Hamiltonian systems perturbed by a small time-dependent periodic forcing term the Melnikov technique is applicable. Our method is essentially based on a one-dimensional unstable manifold. For systems such as the Lorenz equations the trajectories shoot
in and
spiral
out
near
the non-zero
stationary
states. To examine this we either need to modify the method here or to reverse time in the equations. Evidence
for homoclinic
orbits
at non-hyper-
bolic saddle points in real chemical systems has been reported by Gray et al. [9, lo] in the hydrogen-oxygen and carbon monoxide-oxygen reactions. Both of these systems display an explosion limit in well-stirred flow reactors separating steady reaction (low reaction rate) from ignition (high stationary state reaction rate). On crossing this curve of bifurcation points (by raising the reaction temperature) the system jumps into large amplitude oscillations. Further increase in ambient temperature reveals a classic Hopf bifurcation extinction of the limit cycle, which reappears as the temperature is reduced again, without hysteresis.
C. Kaa.-Petersen
and S.K. Scott / Homochic
no=02
10
orbits
1 0 -1
1
no=0
4
0.5
OOC~
0
-
-T
01
02
+,
03
OOL-0
I
02
01
03
tr
(b)
-.. f-R no=005
,/I
1.0
00’0 02
.
I
d
_*--
/
---
,/-
/
/’
^
:
1
I I \ \
00
I--
O2
01
t,
n
.03
Cc) Fig. 8. Accurately computed stationary and periodic solution diagrams for the system various p,,: full lines - stable stationary solutions; dashed lines unstable stationary periodic solutions; open circles- amplitude of unstable periodic solutions.
01
02
t,
03
Cd) in fig. 5 with K~ = 0.05 and K,, = 0.0005 for solutions: full circles - amplitude of stable
C. Kaar-Petersen
and S. K. Scott/
Homoclimc
00 .Ol
02
tr
03
orbits
469
01
02
t,
03
(0
(e) Fig.
8. Continued
Further decrease in temperature causes a lengthening in period, with an apparent faster than logarithmic growth. The large amplitude oscillations are extinguished at exactly the same ambient temperature as that at which they appear on the upward traverse, so there is no hysteresis which would normally be associated with homoclinic orbit formation. Acknowledgements We thank Dr. C. Sparrow for providing us with helpful insight into the origin of homoclinic orbits with saddle points. We are grateful to Professors J.K. Hale and R. Aris for suggesting this problem and to SERC for support (grant GR/D 56990).
References [l] R. At-is and I.G. Kevrekidis, private communication (1987). [2] K. Bar-Eli and R.M. Noyes, Relevance of a two-variable
Oregonator to stable and unstable steady states and limit cycles, to thresholds of excitability and to Hopf vs. SNIPER bifurcations, J. Chem. Phys. 86 (1987) 1927-1937. [3] J. Brindley. C. Kaas-Petersen. J.H. Merkin and SK. Scott, in preparation. [4] E. Doedel, AUTO: Software for continuation and bifurcation problems in ordinary differential equations, Caltech (1986). [5] P. Glendinmng and C. Sparrow, Local and global behaviour near homoclinic orbits, J. Stat. Phys. 35 (1984) 645-697. [6] B.F. Gray and S.K. Scott, Multistability and sustained oscillations in isothermal open systems, J. Chem. Sot. Faraday Trans. 1, 81 (1985) 1563-1567. [7] P. Gray and S.K. Scott, Autocatalytic reactions in the isothermal continuous stirred-tank reactor: isolas and other forms of multistability, Chem. Eng. Sci. 38 (1983) 29-43. [8] P. Gray and SK. Scott, Autocatalytic reactions in the isothermal continuous stirred-tank reactor: oscillations and instabilities, Chem. Eng. Sci. 39 (1984) 1087-1097. [9] P. Gray, J.F. Griffiths and SK. Scott, Proc. R. Sot. London A 394 (1984) 243-258. [lo] P. Gray, J.F. Griffiths and SK. Scott, Proc. R. Sot. London A397 (1985) 21-44. [Ill J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields (Springer, New York, Berlin, Heidelberg, 1983).
470
C. Km.-Petersen
und S.K. Scott/
[12] J.K. Hale, Introduction to dynamic bifurcation, in: L. Salvadori, ed. Bifurcation Theory and Applications, Springer Lecture Notes in Mathematics, no. 1057 (1984). [13] B.D. Hassard, N.D. Kazarinoff and Y.-H. Wan, Theory and Applications of Hopf Bifurcation (Cambridge Univ. Press, 1981). [14] C. Kaas-Petersen, PATH - User’s guide, Centre for Nonlinear Studies, University of Leeds (1987). [15] C. Kaas-Petersen and S.K. Scott, Stationary-state and Hopf bifurcation patterns in isothermal reaction-diffusion systems, Chem. Eng. Sci. 43 (1988) 391-392. [16] S.R. Kay, SK. Scott and P.G. Lignola, The application of singularity theory to isothermal autocatalytic reactions, Proc. R. Sot. London A409 (1987) 433-448. [17] S. Schecter, The saddle-node separatrix-loop bifurcation,
Homoclinic orhts
SIAM J. Math. Anal. 18 (19X7) 1142-1156. [18] S.K. Scott, Isolas, mushrooms and oscillations In isothcrmal autocatalytic reaction-dilfusion equations, Chem. Eng. Sci. 42 (1987) 307-316. [19] SK. Scott and W.W. Farr, Dynamic tine structure in the cubic autocatalator, submitted to Chem. Eng. Sci. [20] C. Sparrow, The Lorenz Equations: Bifurcations, Chaos and Strange Attractors (Springer, New York, Heidelberg. Berlin, 1982). [21] C. Sparrow, The Lorem equations, in: A. Holden, ed. Chaos (Manchester Univ. Press, 19X6). [22] H.B. Stewart and J.M.T. Thompson, Towards a clasaitication of generic bifurcations in dissipative dynamical systems, Dynamics and Stability of Systems 1 (19X6) X7-96.