Slow Manifold Structure in Explosive Kinetics. 2 ... - ACS Publications

J. Phys. Chem. A 2006, 110, 13463-13474

13463

Slow Manifold Structure in Explosive Kinetics. 2. Extension to Higher Dimensional Systems M. Giona,*,† A. Adrover,† F. Creta,‡ and M. Valorani‡ Dipartimento di Meccanica e Aeronautica and Dipartimento di Ingegneria Chimica, Facolta` di Ingegneria, UniVersita` di Roma “La Sapienza” Via Eudossiana 18, 00184 Roma, Italy ReceiVed: June 9, 2006; In Final Form: September 11, 2006

This Article extends the geometric analysis of slow invariant manifolds in explosive kinetics developed by Creta et al. to three-dimensional and higher systems. Invariant manifolds can be characterized by different families of Lyapunov-type numbers, based either on the relative growth of normal to tangential perturbations or on the deformation of m-dimensional volume elements (if the manifold is m-dimensional) and of the complementary (n - m)-elements in the normal orthogonal complement. The latter approach, based on elementary concepts of exterior algebra, is particularly simple because the evolution of the relevant volume elements can be related to suitable local stretching rates, and local analysis can be performed directly from the knowledge of the Jacobian matrix of the vector field. Several examples of bifurcations of the points-atinfinity, which modify the manifold structure, are discussed for 3-D models of exothermic reactions.

1. Introduction The geometric description of the structure of invariant manifolds in chemical reacting systems, divorced from perturbative analyses and expansions, provides many useful suggestions for the understanding of global dynamics.1,2 This is because complex reacting schemes of physicochemical interest are almost never expressed in a canonical singularly perturbed form3,4 in which slow and fast variables are a priori identified. Conversely, geometric methods display a sufficiently high degree of generality to make their results directly applicable to reaction schemes of practical interest. Geometric methods applied to invariant manifold reconstruction have been proposed by Roussel and Fraser,5,6 Davis and Skodje,1 and Adrover et al.7 by focusing on different specific features characterizing the invariant structures. This Article develops further the geometric approach proposed by Creta et al.8 for 2-D combustion models and extends its range of applicability to generic dynamical systems of the form

dz ) F(z) dt

(1)

where z ∈ Rn and n > 2. The direct extension of the approach developed by Creta et al.,8 based on the scaling of the ratio of normal to tangent perturbations, shows some technical and practical shortcomings for higher dimensional systems, for the simple reason that in Rn, with n > 2, the tangent spaces and the normal spaces to an invariant manifold are no longer 1-D. This issue is addressed in section 2. The extended method proposed in section 3 makes use of the tools of exterior algebra10-12 by considering the evolution of m-dimensional measure elements in the tangent spaces to the manifold (where m is the dimension of the manifold) and (n - m)-dimensional measure elements in the complementary * Corresponding author. E-mail address: [email protected] (M.G.). † Dipartimento di Ingegneria Chimica. ‡ Dipartimento di Meccanica e Aeronautica.

normal spaces. This approach, which is proposed and analyzed in section 3, provides a simple and efficient characterization of the dynamic features of invariant manifolds based on the properties of suitable Lyapunov-type numbers. Moreover, this approach is particularly simple and suitable to practical application to generic kinetic models, because the relevant quantities can be obtained exclusively from the local Jacobian matrix, i.e., from the local stretching rates for measure-element evolution. The occurrence of slow manifolds of higher dimensions makes the bifurcational analysis associated with the behavior of the equilibria at infinity8 more rich and articulated than in the 2-D case. Bifurcational properties of the points-at-infinity controlling the structure of the slow invariant manifold are illustrated in section 4 by considering the case of a 3-D system associated with the dynamics of two exothermic reactions in series. For this system, a complete bifurcational analysis is presented and the results are explained by means of the stretching properties experienced by normal and tangential measure elements along the invariant manifolds. For a presentation of the relevant mathematical tools associated with vector dynamics and with the definitions of slow invariant manifolds, the reader is referred to ref 8. 2. From n ) 2 to Higher Dimensional Systems Moving from 2- to 3-dimensional systems (or higher), the extension of the definitions and the characterization of slow invariant manifolds presented in ref 8 requires a significant (additional) amount of conceptual and formal complexity. This is due to the fact that although in 2-D systems the invariant manifolds W of dynamic interest are 1-D structures, and consequently, both the tangent subspaces Cz and the normal subspaces Nz at any point z ∈ W are 1-D, this is intrinsically not true for n > 2. Consider an m-dimensional invariant manifold W (m < n) for the system eq 1 with n > 2, and let z ∈ W be a generic point of the manifold. In this case, the tangent subspace Cz to the manifold at z is m-dimensional and invariant under vector dynamics, and the normal subspace Nz is (n - m)-dimensional and not invariant. Following Fenichel9 and Creta et al.,8 it is

10.1021/jp063608o CCC: $33.50 © 2006 American Chemical Society Published on Web 11/29/2006

13464 J. Phys. Chem. A, Vol. 110, No. 50, 2006

Giona et al.

TABLE 1: r/ω-Lyapunov Numbers along the Invariant Manifolds of a Linear Autonomous 3-D System Λω

ΛR

W (1) 1

λ2/λ1 > 1

λ3/λ1 > 1

(1) 2 (1) 3 (1)

λ1/λ2 < 1

λ3/λ2 > 1

λ1/λ3 < 1

λ2/λ3 < 1

λ2/λ1 > 1 λ3/λ1 > 1

λ2/λ3 < 1 λ3/λ2 > 1

W (2) 1,3

λ2/λ1 > 1

λ2/λ3 < 1

λ1/λ2 < 1

λ1/λ3 < 1

λ3/λ1 > 1

λ1/λ3 < 1

manifold W

W W W (2) 1,2 W W

(2) 2,3 (2)

still possible to introduce some Lyapunov-type numbers based on the relative properties of normal and tangent vectors to W. Let

nt(z) ) Πφt(z)[φ/t (z) n0] ct(z) ) φ/t (z)c0

n0 ∈ N z c0 ∈ Cz

(2)

where φt is the phase flow associated with eq 1, φ/t (z) ) ∂φt(z)/∂z, and Πz indicates the normal projector at z, i.e, the operator mapping any vector into its component lying in the normal subspace Nz. The R/ω-Lyapunov-type number can be defined in Rn as

Λω ) lim tf∞

sion, is the behavior for t f -∞, i.e., the occurrence of the highest possible value of the R-Lyapunov number ΛR. On the basis of the data in Table 1, it is possible to provide the following definitions of global and generalized slow manifolds of dimension m for a generic n-dimensional system. Given the dynamical system eq 1, a global slow manifold of dimension m is an invariant, exponentially attracting, stable m-dimensional manifold for which ΛR attains its maximum value greater than 1, and Λω > 1. A global m-dimensional fast manifold is an invariant, exponentially attracting m-dimensional manifold for which ΛR < 1 and Λω attains its smallest value less than 1. The concept of generalized slow manifolds can be defined by removing the condition that Λω should be greater than 1 from the definition of a global slow manifold. According to these definitions, the global 1-D slow manifold (2) for the above linear system is given by W (1) 1 , and W 1,2 is the global 2-D slow manifold. To give a numerical example of the definition of R/ωLyapunov-type numbers to nonlinear models, let us consider the evolution of two exothermic reactions in series A f B f C in a batch-jacketed reactor,13,14 henceforth referred to as the 3-D Semenov model. This model will be used throughout this Article as a prototype for higher dimensional combustion systems as it regards the geometry of the invariant manifolds. Each reaction step is elementary and of first-order with respect to its reactant. The balance equations for reactant concentrations cA and cB read

supn0∈Nz log||nt(z)||

dcA o ) -kA e-EA/RTcA dτ

supc0∈Cz log||ct(z)|| ΛR ) lim

tf-∞

supn0∈Nz log||nt(z)|| supc0∈Cz log||ct(z)||

dcB o o ) kA e-EA/RTcA - kB e-EB/RTcB dτ

z ∈ W (3)

where the supremum is taken over all the initial vectors lying in the normal subspace (in the numerator) and in the tangential subspace to W at z (in the denominator). To analyze how this definition of the Lyapunov-type numbers applies in practical calculations, consider a 3-D linear constant coefficient dynamical system, F(z) ) Az, where the coefficient matrix possesses three distinct negative eigenvalues {-λ1, -λ2, -λ3}, with λ1 < λ2 < λ3, and let e1, e2, and e3 be the corresponding eigenvectors. For this system, the invariant manifolds associated with directions of the eigenspaces are the three 1-D eigenmanifolds W (1) h ) {z|z ) ξeh, ξ ∈ (-∞, ∞)}, h ) 1, 2, 3 and the three 2-D eigenmanifolds W (2) h,k ) {z|z ) ξeh + ηek, ξ, η ∈ (-∞, ∞)}, h ) 1, 2, 3 < k. By definition, W (1) 1 and W (2) 1,2 are respectively the global slow 1- and 2-D manifolds of the system. Table 1 reviews the values of the R/ω-Lyapunov-type numbers on these eigenmanifolds and on generic 1-D (W (1)) and 2-D (W (2)) invariant manifolds. We observe that for t f ∞ the controlling normal contraction rate is λ1 everywhere except (2) on W (1) 1 and on W h,k with h ) 1. Conversely, for t f -∞, the controlling normal elongation rate, because of stability reversal, (2) is λ3 everywhere except on W (1) 3 and on W h,k with k ) 3. Clearly, for a generic 1-D invariant manifold W (1) (an orbit), its behavior for t f ∞ will mimic that of W (1) 1 , and for t f -∞, that of W (1) 3 . A similar reasoning applies for a generic 2-D invariant manifold W (2). As in the 2-D case,8 a slow invariant manifold is characterized by R/ω-Lyapunov-type numbers both greater than 1. Moreover, the discriminating feature of a slow invariant manifold, which distinguishes it from other invariant manifolds of equal dimen-

FcV

dT o o ) (-∆HA)kA e-EA/RTcA + (-∆HB)kB e-EB/RTcB dτ Ua(T - Tc) (4)

where T is the temperature and the specific reaction enthalpies ∆HA and ∆HB are negative. With the dimensionless quantities o γA ) EA/RTc, γB ) EB/RTc, t ) τkA e-γA, x ) cA/cref, y ) cB/ cref, and z ) (T - Tc)γA/Tc where cref is a reference concentration and Tc the coolant temperature, eq 4 becomes

dx ) -fA(z)x dt dy ) fA(z)x - ζ fB(z)y dt dz fA(z)x + hfB(z)y - δz ) dt

(5)

where

fA(z) ) exp o

( ) γAz γA + z o

fB(z) ) exp

( ) γB z γA + z

(6)

and ζ ) kB eγA-γB/kA > 0, h ) (-∆HB)ζ/(-∆HA) > 0, P ) o (-∆HA)crefγA/FcVTc, Q ) UaeγA/FcVkA , ) 1/P, and δ ) Q/P. Figure 1A depicts the phase plot of the 3-D Semenov model for a fixed set of parameter values, which gives rise to a global 1-D slow manifold. The global 1-D slow manifold is depicted in panel A with a thicker line and is obtained by means of material line advection (MLA).8 Figure 1B depicts the ratio log||nt||/log||ct|| starting from a generic normal vector n0 ∈ Nz.

Slow Manifold Structure in Explosive Kinetics. 2

J. Phys. Chem. A, Vol. 110, No. 50, 2006 13465

Figure 1. (A) 1-D global slow manifold (thicker line) for the 3-D Semenov model with γA ) 5, γB ) 8, ζ ) 10, h ) 2, δ ) 1, ) 10-3. The figure depicts some orbits attracted by the manifold. (B) log||nt||/log||ct|| vs t along the invariant manifold depicted in panel A.

Equation 10 is the basic equation for understanding the difference between two-dimensional and higher models. In fact, for n ) 2, vn belongs to a 1-D subspace, the unit tangent vector of which nˆ is uniquely determined by the vector field F at the point. Correspondingly, from eq 10, one obtains Figure 2. Schematic evolution of normal and tangential vectors to an (1) invariant 1-D manifold W .

In this case, > 1. A similar analysis performed on the backward evolution of tangent and normal vectors, not shown for brevity, yields ΛR > 1, confirming the global nature of the 1-D slow manifold, depicted in panel A. The application of the definitions of R/ω-Lyapunov-type numbers expressed by eq 3 exhibits some practical shortcomings in higher dimensional (n > 2) systems that are worth addressing. Consider the dynamics of tangent and normal vectors along a 1-D invariant manifold W (1) (i.e., an orbit) of a generic n-dimensional system (Figure 2). The tangent sub-bundle to W (1) is invariant under vector dynamics (i.e., under φ/t (z)), and the normal sub-bundle is not, this being the reason for the application of the local normal projector Πφt(z) in the definition of nt. Let us analyze in detail the dynamics of a generic vector v, starting from an initial vector v0 ∈ Tz having components both in the tangent and in the normal subspaces. For any time t > 0, let v(t) ) φ/t (z)v0. The vector v(t) can be expressed as

||vn(t)|| ) ||vn(0)|| exp(

∫0tων(z(t′)) dt′)

Λω

v ) aF + vn

(7)

where F ) F(φt(z)) is the vector field, a a scalar depending on time, and vn(t) ) Πφt(z)[v(t)] is the normal component of v. By definition, both v and F satisfy the equation for vector dynamics describing the evolution of a generic vector in the tangent bundle under the action of the vector field

dv ) F*v dt

dF ) F*F dt

(8)

where F*(z) ) ∂F(z)/∂z. After differentiating eq 7 with respect to time and substituting into it the expressions for the time derivatives of v and F (eq 8), it follows that

dvn da ) F*vn - F dt dt

(9)

Equation 9 confirms what was stated above, namely that the normal sub-bundle is not invariant, as in eq 9 the extra term F da/dt appears aligned in the tangential direction. To get rid of this term, take the scalar product with respect to vn, to obtain (after some steps)

d||vn|| ) (F*vˆ n, vˆ n)||vn|| dt

where vˆ n ) vn/||vn|| is the unit vector associated with vn.

(10)

ων ) (F*nˆ , nˆ ) (11)

where z(t) indicates the trajectory along the manifold and ων is the normal stretching rate. It follows that the R/ω-Lyapunov numbers can be expressed as a function of the tangential and normal stretching rates because

log

( )

||ct|| ) ||c0||

∫0tωτ(z(t′)) dt′ log

( )

||nt|| ) ||n0||

∫0tων(z(t′)) dt′

(12)

where ωτ ) (F*cˆ , cˆ ) is the tangential stretching rate, with cˆ ) F/||F||. Conversely, if n > 2, the evolution of the normal component of a generic vector cannot be expressed a priori with respect to a given normal stretching rate (unless the vector evolution is not explicitly accounted for through an expression equivalent to eq 11), because the normal subspaces are no longer 1-D, and the unit normal direction vˆ n(t) along the orbit depends on the orientational dynamics within the normal sub-bundle. This makes the analysis of the evolution of normal perturbations more cumbersome, unless some form of further specification in the definition of the characteristic normal stretching rates is not added. A way for bypassing this shortcoming is discussed in the next section. 3. m-Forms, Stretching Rates, and Invariant Manifold Properties An alternate way to provide a geometrical characterization of the invariant manifolds is to consider the evolution of m-dimensional measure elements. This section develops this approach, which leads to a new definition of the Lyapunovtype numbers, and analyzes the differences and the advantages of this approach with respect to the analysis developed in section 2 based on the evolution of the norms of normal and tangential perturbations. Consider an invariant, exponentially attracting m-dimensional manifold W (m) for eq 1 with m < n. The tangent subspace Cz to W (m) at any point z ∈ W is m-dimensional and invariant, and the (n - m)-dimensional normal subspace Nz is not. The idea is to characterize the dynamical properties of W (m) by means of the stretching rates of m-dimensional measure elements constructed upon vectors lying on Cz (where m is the dimension of the manifold itself) and of (n - m)-dimensional measure


Giona et al.

elements deriving from vectors lying on the normal subspace. This can be achieved by introducing exterior algebraic concepts (which are succinctly reviewed in Appendix A). (2) (m) Let c(1) 0 , c0 , ..., c0 be a family of m linearly independent vectors spanning the subspace Cz0 at z0 ∈ W (m), and define c(h)(t) ) φ/t (z0) c0(h), h ) 1, ..., m. Because the tangent subbundle is invariant, each c(h)(t) lies in the subspace Cφt(z0) for any time t g 0. The m-dimensional measure element spanned by c(1), c(2), ..., c(m) is given by the m-form c(1) ∧ c(2) ∧ ... ∧ c(m). The time evolution of this m-form can be obtained by differentiating it with respect to time and by enforcing the evolution equation for each individual vector dc(h)/dt ) F*c(h),

d[c(1) ∧ c(2) ∧ ... ∧ c(m)] dc(1) ) ∧ c(2) ∧ ... ∧ c(m) + dt dt dc(2) dc(m) ∧ ... ∧ c(m) +... + c(1) ∧ c(2) ∧ ... ∧ c(1) ∧ dt dt ) F*,∧m[c(1) ∧ c(2) ∧ ... ∧ c(m)]

(13)

In eq 13, we have introduced the operator F*,∧m acting in the exterior m-space associated with the tangent subspaces Cz to indicate the action of the differential operator on the m-forms. It should be observed that the action of the operator F*,∧m depends solely on the m-form c(1) ∧ c(2) ∧ ... ∧ c(m) and not on c(1), c(2), ..., c(m) individually.15 By taking the scalar product in the m-exterior space (see Appendix A) with respect to c(1) ∧ c(2) ∧ ... ∧ c(m), we find that it follows after some algebraic manipulations that

d||c(1) ∧ ... ∧ c(m)||m

A similar approach applies to the (n - m)-measure element (2) generated by (n - m) linearly independent vectors v(1) 0 , v0 , ..., (n-m) v0 initially lying in the normal (n - m)-dimensional subspace Nz0 to W (m). The analysis is slightly more elaborate in this case, because a generic normal vector generates in its evolution a vector possessing nonvanishing components both in the normal and in the tangential subspaces of the image point. (h) (h) Let v(h)(t) ) φ/t (z0)v(h) 0 , vν (t) ) Πφt(z0)[v (t)], h ) 1, ..., n m, where Πφt(z0) is the normal projector at the point φt(z0). The final result of this calculation is that

(n-m) d||v(1) ||n-m ν ∧ ... ∧ vν

)

dt n-m

[

(n-m) ων,h(z(t))]||v(1) ||n-m ∑ ν ∧ ... ∧ vν h)1

where ων,h, h ) 1, ..., n - m are the n - m normal stretching rates defined starting from a generic basis of n - m orthonormal unit vectors nˆ 1, ..., nˆ n-m spanning Nz

ων,h(z) ) (F*(z) nˆ h(z), nˆ h(z))

h ) 1, ..., n - m

n-m m

)[

(F* ˆt(h)(z(t)), ˆt(h)(z(t)))]||c(1) ∧ ... ∧ c(m)||m ∑ h)1 m

)[

ωτ,h(z(t))]||c(1) ∧ ... ∧ c(m)||m ∑ h)1

(14)

where ˆt(h), h ) 1, ... m is a system of m orthonormal unit vectors spanning Cz(t), where z(t) ) φt(z0) and

ωτ,h(z) ) (F*tˆ(h)(z), ˆt(h)(z))

h ) 1, ..., m

dµτ,m(t) dt

∑

m

ω(m) τ (z) )

(16)

h)1

∑

ω(n-m) (z) ) ν

m

µτ,m(t) ) µτ,m(0) exp[

∑∫0 ωτ,h(z(t′)) dt′] h)1 t

(17)

∑ ων,h(z)

h)1

z ∈ W (m) (21) which provide an intrinsic dynamical characterization of the local behavior at points z ∈ W (m). 3.1. Exterior Lyapunov-Type Numbers and Invariant Manifolds. On the basis of the evolution of the measures µτ,m(t) and µν,(n-m)(t) associated with an invariant m-dimensional manifold W (m), it is possible to introduce the exterior R/ωR ω and ΛE,m defined as Lyapunov-type numbers ΛE,m

m log µν,(n-m)(t) (n - m) log µτ,m(t) R ) lim ΛE,m tf∞

and therefore

(20)

n-m

ωτ,h(z)

h)1

tf∞

ωτ,h(z(t))]µτ,m(t)

t

To sum up, the local stretching properties along an mdimensional invariant manifold can be expressed by means of the m- and (n - m)-dimensional stretching rates

ω ΛE,m ) lim

m

)[

∑ ∫0 ων,h(z(t′)) dt′] h)1

(15)

are the tangential stretching rates that can be defined starting from this orthonormal tangential system. The value of each individual ωτ,h depends on the chosen basis, and their sum (appearing in eq 14) is independent of the orthonormal basis chosen and is a local intrinsic property of the action of the dynamical system along the invariant manifold. If we indicate with µτ,m(t) ) ||c(1)(t) ∧ ... ∧ c(m)(t)||m the measure of the m-dimensional tangential measure element, from eq 14 it follows that

(19)

As for the tangential stretching rates, the value of each individual ων,h depends on the chosen normal basis for Nz, and their sum is independent of the basis itself. From eq 19, it follows that the (n - m)-measure of the normal (n-m) measure element µν,(n-m)(t) ) ||v(1) (t)||n-m, ν (t) ∧ ... ∧ vν (1) (n-m) spanned by vν , ..., vν , satisfies the equation

µν,(n-m)(t) ) µν,(n-m)(0) exp[

dt

(18)

m log µν,(n-m)(t) (n - m) log µτ,m(t)

(22)

The dimension m and the co-dimension (n - m) of the manifold W (m) enter explicitly in this definition in order to “homogenize” the scaling of the measure elements of different



dimensions. Indeed, the argument of the limits appearing in eq 22 can be expressed as

m log µν,(n-m)(t) (n - m) log µτ,m(t)

)

log[µν,(n-m)(t)]1/(n-m)

(23)

log[µτ,m(t)]1/m

and therefore, the exterior Lyapunov-type numbers can be viewed as the ratio of the logarithms of two characteristic vector lengths in the normal and tangential subspaces, [µν,(n-m)(t)]1/(n-m) and [µτ,m(t)]1/m, defined starting from the evolution of complementary measure elements. It follows from the analysis developed above that the logarithms appearing in eq 22 can be expressed by means of (n-m) the integrals of the stretching rates ω(m) (eqs 17 and τ and ων (m) 20) along system trajectories lying in W , i.e.,

log

( ) µτ,m(t)

µτ,m(0)

)

∫

log

(

µν,(n-m)(t)

µν,(n-m)(0)

)

)

(z(t′)) dt′ ∫0tω(n-m) ν

(24)

This is the main advantage in adopting the definition of an exterior Lyapunov-type number. Moreover, it follows from the definition of eq 21 that (n-m) (z) ) Trace[F*(z)] ω(m) τ (z) + ων

(25)

n where Trace[F*(z)] ) ∑h)1 F*h,h(z) is the trace of the Jacobian matrix. Therefore, in practical applications it is not necessary to estimate both ω(m) and ω(n-m) , but solely one of these τ ν stretching rates, because the remaining one follows from eq 25. This result is particularly useful in the analysis of 1-D invariant manifolds because

ˆ (z), Fˆ (z)) ω(1) τ (z) ) (F*(z) F (z) ) Trace[F*(z)] - ω(1) ω(n-1) ν τ (z) (26) where Fˆ ) F/||F||. Let us apply the definition of the exterior Lyapunov numbers to an n-dimensional linear system dz/dt ) Az, the coefficient matrix of which admits n distinct negative eigenvalues -λ1, -λ2, ..., -λn with λ1 < λ2 < ... < λn associated with the eigenvectors e1, e2, ..., en. For fixed m > 1, there exist n!/m!(n - m)! different m-dimensional eigenmanifolds W i(m) 1,...,in spanned by ei1, ei2, ..., ein with i1 < i2 < ... < in and passing through the origin. On each of these manifolds, the R- and ω-exterior Lyapunov numbers coincide and are given by ω R (W i(m) ) ) ΛE,m (W i(m) )) ΛE,m 1,...,in 1,..., in

[

respectively), it follows that ω R (1) ΛE,1 (W (1) 1 ) ) ΛE,1(W 1 ) )

ω R (2) (W (2) ΛE,2 1,2) ) ΛE, 2(W 1,2) )

t

ω (m)(z(t′)) dt′ 0 τ

Figure 3. (A) ω(2) ν vs the curvilinear abscissa s along the global slow invariant manifold for the 3-D Semenov model depicted in Figure 1. (B) |r(1) E | vs s along the same global slow manifold.

]

Trace(A) - (λi1 + ... + λin) m (27) λi1 + ... + λin (n - m) (2) (n-1) Specifically, for W (1) 1 , W 1,2, ...,W 1,...,n-1 (corresponding to the slow m ) 1, 2, ..., (n - 1)-dimensional manifolds,

λ2 + ... + λn (n - 1)λ1

2(λ3 + ... + λn) (n - 2)(λ1 + λ2)

ω (n-1) R (n-1) ΛE,n-1 (W 1,...,n-1 ) ) ΛE,n-1 (W 1,...,n-1 ))

>1

>1

(28)

(n - 1)λn >1 λ1 + ... + λn-1

Conversely, for any other m-dimensional invariant manifold W (m), the exterior Lyapunov numbers attain generically the expression ω ) ΛE,m

m(λm+1 + ... + λn) >1 (n - m)(λ1 + ... + λm) m(λ1 + ... + λn-m) R ΛE,m ) < 1 (29) (n - m)(λn-m+1 + ... + λn)

Therefore, the m-dimensional global slow manifolds are characterized by the occurrence of the maximum value of the exterior R/ω-Lyapunov numbers greater than 1. This gives rise to the following definition of global and generalized slow manifolds based on the measure-element scaling. Given the dynamical system eq 1, a global slow manifold of dimension m is an invariant, exponentially attracting, stable R ω m-dimensional manifold for which ΛE,m and ΛE,m are greater than 1. A global m-dimensional fast manifold is an invariant, exponentially attracting m-dimensional manifold for which R ω ΛE,m < 1 and ΛE,m < 1. The concept of a generalized slow ω manifold can be defined by removing the condition that ΛE,m should be greater than 1 from the definition of a global slow manifold. As an example, Figure 3A shows the behavior of the normal stretching rate ω(2) ν (s) along the curvilinear abscissa s of the global slow manifold depicted in Figure 1. As expected, ω(2) ν (s) is uniformly negative, meaning that 2-D normal measure elements shrink exponentially along the manifold. Figure 3B shows the stretching ratio r(1) E ,

r(1) E (z) )

ω(2) ν (z) 2ω(1) τ (z)

(30)

along the manifold. The absolute value of r(1) E is reported, because the tangential stretching rates attain positive values close to the explosion. It can be observed that there are portions of the manifold along which the stretching ratio attains values less


Giona et al.

than 1, indicating that local inversion in the behavior of a normal/tangential stretching rate occurs. This phenomenon is thoroughly addressed in section 4, which analyzes the bifurcations occurring in the 3-D Semenov system and how these bifurcations modify the structure of the slow invariant manifolds. 4. Bifurcations and Slow Manifold Structure In this section, we analyze the structures and properties of slow invariant manifolds in the 3-D Semenov model by combining the compactification technique with the scaling theory of exterior measure elements developed in section 3. The properties of invariant manifolds can be further addressed by considering the compactification of the phase space, i.e., by introducing the following coordinate transformation

uh )

zh

x

h ) 1, ..., n,

(31)

n

1+

∑zk2

k)1

n mapping Rn onto the n-dimensional unit sphere Sn ) {u|∑h)1 2 uh e 1}. Correspondingly, the introduction of the Poincare´ projected system (Pp-system) associated with eq 1

duh dt

n

) (1 -

∑

n

uk2)1/2 [Fh - uh

k)1

∑ uk F k ]

(32)

eigenvalues {µ(i)} and eigenvectors {w(i)}

k)1

and the analysis of its behavior close to the boundary ∂S 1n ) n uh2 ) 1} permits us to investigate the behavior of the {u|∑h)1 original system eq 1 at infinity. The introduction of the Pp-system eq 32 makes it possible to analyze the structure and properties of global invariant manifolds of eq 1 in terms of the properties of the equilibrium points-at-infinity, i.e., the equilibrium points u∞eq ∈ ∂S 1n of the associated Pp-system. For example, for a dynamical system possessing a unique globally attracting equilibrium point zeq ) 0, a global/generalized 1-D invariant manifold is a heteroclinic orbit of the Pp-system connecting ueq ) 0 (corresponding to zeq in the transformed coordinates) to one of the equilibrium points-at-infinity u∞eq, such that the exterior R/ω-Lyapunov numbers ΛωE,1 and ΛRE,1 (controlled by the behavior close to ueq and u∞eq, respectively) possess specific properties. To comment on this issue, let us consider again an autonomous linear system dz/dt ) Az analyzed in section 2. The associated Pp-system attains the form

duh dt

n

)

∑

k)1

n

Ahkuk - uh

∑ Akmukum

Figure 4. Analysis of the Pp-system associated with a 3-D linear system. (A)-(F): points-at-infinity. (A) and (B): u∞,( eq,1 ) ((-1/x5, 2/x5, 0) (stable nodes on ∂S 1n). (C) and (D): u∞,( eq,2 ) ((1/x2, 0, 1/ x2) (saddles on ∂S 1n). (E) and (F): u∞,( eq,3 ) ((0, - 2/x5, 1/x5) (unstable nodes on ∂S 1n). The thick line connecting A and B is the 1-D slow manifold. The curve γ is the set of heteroclinic orbits connecting ∞( 1 1 u∞,( eq,1 with ueq,2 on ∂S n. The shaded region intersecting ∂S n at γ is the 2-D slow manifold. Some orbits of the Pp-system are also drawn in order to highlight the role of 1- and 2-D slow manifolds.

(33)

k,m)1

Elementary algebraic manipulations yield 2n + 1 equilibrium points for the Pp-system: (i) the stable equilibrium point ueq ) 0 (corresponding to the unique stable equilibrium point zeq ) 0 of the original system) characterized by the same eigenvalues { - λ1, - λ2, ..., - λn}, λ1 < λ2 < ...< λn and eigenvectors {e1, e2, ..., en} of the original linear system eq 1 and (ii) 2n equilibrium points-at-infinity u∞,( eq ) (eh, h ) 1, ..., n. The equilibrium points-at-infinity correspond to the invariant directions associated with the 1-D eigenmanifolds of the system. ∞ Let ueq,i be a generic equilibrium point-at-infinity. It can be ∞ is characterized by the following set of shown that ueq,i

µ(i) k ) -λk + λi µ(i) i ) 2λi > 0

w(i) k ) ek - (ek, ei)ei k ) 1, 2, ..., i - 1, i + 1, ..., n w(i) i ) ei

i ) 1, ..., n

(34)

This implies: (1) Each equilibrium point-at-infinity is unstable on the n-dimensional unit sphere Sn. (2) Each equilibrium pointat-infinity can be either a stable one or an unstable one or a saddle point if one considers the dynamics of the Pp-system restricted to the boundary ∂S 1n of the unit sphere Sn. More ∞,( precisely, u∞,( eq,1 are stable nodes, ueq,k , k ) 2, ..., n - 1 are ∞,( saddle points, and ueq,n are unstable nodes. (3) The global invariant 1-D slow manifold of the original linear system W (1) 1 ) span{e1} is, in the transformed coordinates {uh}, the heteroclinic orbit of the corresponding Pp-system connecting ueq with u∞eq,1, i.e., with the stable node on ∂S 1n. (4) The global invariant 2-D slow manifold of the original linear system W (2) 1,2 ) span{e1, e2} is, in the transformed coordinates {uh}, the 2-D invariant manifold of the corresponding Pp-system intersecting ∂S 1n at the heteroclinic connection on ∂S 1n between the stable ∞ and the first saddle node u∞eq,2 (characterized by only node ueq,1 one positive eigenvalue for the dynamics of the Pp-system restricted to the boundary ∂S 1n). For example, Figure 4 shows the 1-D and 2-D slow manifolds of the Pp-system associated with a 3-D linear system with eigenvalues {-λ1, -λ2, -λ3} ) {-1, -10, -100} and eigenvectors e1 ) (-1/x5, 2/x5, 0), e2 ) (1/x2, 0, 1/x2), e3 ) (0, - 2/x5, 1/x5). The points-at∞,( ) (e1 (points A and B, stable nodes on ∂S 1n), infinity are ueq,1 ∞,( ∞,( ueq,2 ) (e2 (points C and D, saddle points on ∂S 1n), and ueq,3 1 ) (e3 (points E and F, unstable nodes on ∂S n). The 1-D slow invariant manifold (thick line connecting A and B) is the union of the two heteroclinic orbits connecting ueq ) 0 with u∞,+ eq,1 and ueq ) 0 with u∞,. The 2-D slow invariant manifold is the eq,1



plane (shaded region) passing through ueq and intersecting ∂S 1n at the circle γ representing the set of heteroclinic orbits ∞,( ∞,( with ueq,2 on ∂S 1n. connecting ueq,1 The introduction of the Pp-system, the computation of equilibrium points-at-infinity, and the analysis of their stability play an important role also for nonlinear systems as it regards the identification and characterization of global/generalized slow manifolds. For a dynamical system possessing a unique globally attracting equilibrium point zeq ) 0, the global/generalized 1-D slow invariant manifold is, in the transformed coordinates, the heteroclinic orbit of the Pp-system connecting ueq ) 0 to the equilibrium point-at-infinity u∞eq,1 that is a stable node for the dynamics of the Pp-system restricted to ∂S 1n. This implies that, in the original coordinates, the global/ generalized 1-D slow invariant manifold is a curve passing through zeq ) 0 and tangent, at infinity, to the direction associated with the equilibrium point u∞eq,1. Analogously, the global/generalized 2-D slow invariant manifold is, in the transformed coordinates, a 2-D surface intersecting ∂S 1n along the curve γ representing the heteroclinic connection on ∂S 1n ∞ between the stable node ueq,1 and the first saddle point u∞eq,2. This implies that, in the original coordinates, the global/ generalized 2-D slow invariant manifold is a surface tangent, at infinity, to the plane spanned by the directions associated ∞ . with the equilibrium points-at-infinity u∞eq,1 and ueq,2 Similarly, the m-dimensional slow manifold (with m < n) turns out to be a hypersurface tangent, at infinity, to the hyperplane spanned by the directions associated with the stable ∞ ∞ node u∞eq,1 and the first m - 1 saddle points ueq,2 , ..., ueq,m . The nature “global” or “generalized” of invariant slow manifolds can be defined by computing the exterior R/ωω (controlled by the behavior close to Lyapunov numbers ΛE,m R zeq and ΛE,m (controlled by the behavior at infinity), i.e., along the asymptotic directions associated with the equilibrium pointsat-infinity). Moreover, given the fundamental role of the equilibrium points-at-infinity, it is possible to investigate the influence of model parameters on the global behavior of the system and on its invariant manifold structures by analyzing possible local “bifurcations” of the points-at-infinity occurring when model parameters vary. This issue will be addressed in the next subsection where the 3-D Semenov model is analyzed in detail enforcing the observations presented above. 4.1. Three-Dimensional Semenov Model. The 3-D Semenov model eq 5 possesses a unique stable equilibrium point zeq ) 0 characterized by the following eigenvalues/eigenvectors

λ1 ) -1

(

e1 ) C1

(δ - )(ζ - 1) (δ - ) , ,1 ζ-1+h ζ-1+h

(

λ2 ) -ζ

e2 ) C2 0,

λ3 ) -δ/

)

δ-ζ ,1 h

)

e3 ) (0, 0, 1)

(35)

where C1 and C2 are normalization constants. The analysis of the Pp-system is particularly simple for the 3-D Semenov model as for z f ∞, fA(z) f eγA ) f ∞A, fB(z) f eγB ) f ∞B , and therefore, the behavior at infinity is described by the linear system dz/dt ) A∞z, where

(

-f ∞A 0 0 -ζf ∞B 0 A∞ ) f ∞A f ∞A-1 hf ∞B -1 -δ-1

)

(36)

TABLE 2: Position and Stability of Equilibrium Points-at-Infinity of the 3-D Semenov Model for 1 < ζ < δ/E < ∞a δ/

stable ∞ A

u∞eq,3 u∞eq,1 u∞eq,1

ζ < δ/ < f f ∞A < δ/ < ζf ∞B ζf ∞B < δ/

saddle u∞eq,1 u∞eq,3 u∞eq,2

∈F ∈F ∈F

∉F ∈F ∈F

unstable u∞eq,2 ∉F ∞ ueq,2 ∉F ∞ ueq,3 ∈F

a The symbol F ) {(u , u , u ) | u > 0, u > 0} indicates the 1 2 3 1 2 physical region of positive reactant concentrations.

The matrix A∞ admits the following eigenvalues and eigenvectors

λ∞1 ) -f ∞A e∞1 ) C∞1

(

)

(δ - f ∞A)(f ∞A - ζf ∞B ) -(δ - f ∞A) ,1 , f ∞A(f ∞A - ζf ∞B - hf ∞B ) f ∞A - ζf ∞B - hf ∞B

(

e∞2 ) C∞2 0,

λ∞2 ) -ζf ∞B λ∞3 ) -δ/

)

δ - ζf ∞B ,1 hf ∞B

e∞3 ) (0, 0, 1)

(37)

The equilibrium points-at-infinity of the Pp-system are ∞,( associated with the eigendirections of A∞, i.e., ueq,h ) (e∞h , h ) 1, ..., 3. By considering the dynamics of the Pp-system restricted to the boundary ∂S 13 of the unit sphere, we find that the stability of the equilibrium points-at-infinity, according to eq 34, is controlled by the following eigenvalues ∞,( ) (e∞1 ueq,1

∞ ∞ ∞ ∞ µ(1) 2 ) λ2 - λ1 ) (f A - ζf B ) ∞ ∞ ∞ µ(1) 3 ) λ3 - λ1 ) (f A - δ/)

∞,( ueq,2 ) (e∞2

∞ ∞ ∞ ∞ µ(2) 1 ) λ1 - λ2 ) -(f A - ζf B ) ∞ ∞ ∞ µ(2) 3 ) λ3 - λ2 ) (ζf B - δ/)

∞ u∞,( eq,3 ) (e3

∞ ∞ ∞ µ(3) 1 ) λ1 - λ3 ) -(f A - δ/) ∞ ∞ ∞ µ(3) 2 ) λ2 - λ3 ) -(ζf B - δ/) (38)

In order to illustrate how the relative position and stability of points-at-infinity actually controls the global behavior and the structure of invariant manifolds of the 3-D Semenov model, let us consider the following set of parameter values 1 < ζ and f ∞A < ζf ∞B and analyze system properties by letting the parameter δ/ vary in the range ζ < δ/ < ∞. 4.2. Transcritical Bifurcations at Infinity. Given that 1 < ζ < δ/, at the equilibrium point zeq ) 0, span(e1) represents the 1-D slow eigenspace, and span (e1, e2) is the 2-D slow eigenspace. At infinity, i.e., on the boundary ∂S 13, the position and stability of the equilibrium points for f ∞A < ζf ∞B are summarized in Table 2 where the symbol F indicates the physical region F ) {(u1, u2, u3)|u1 > 0, u2 > 0}, i.e., the region of positive concentration values. ∞ For ζ < δ/ < f ∞A, the equilibrium points-at-infinity ueq,3 ∞ and ueq,1 are a stable node and a saddle point, respectively. The u3-axis (corresponding to the z axis in the original coordinates) is an invariant 1-D manifold for the Pp-system. For ζ < δ/ < f ∞A, it represents the heteroclinic connection between the stable ∞ and the equilibrium point ueq ) 0. Therefore, the node ueq,3 z-axis, in the original coordinates, is the generalized slow manifold of the Semenov model for ζ < δ/ < f ∞A because it is


Giona et al.

Figure 5. Movements of points-at-infinity u∞eq,1 and u∞eq,2 on the boundary ∂S 13 for increasing values of the parameter δ/. Arrows indicate increasing values of δ/. Points u∞eq,1 (saddle point for δ/ < f ∞A, dotted line) and u∞eq,3 (stable point for δ/ < f ∞A) coincide for δ/ ) f ∞A and then exchange their stabilities. Points u∞eq,2 (unstable node for δ/ < ζf ∞B , dot-line curve) and u∞eq,3 (saddle point for f ∞A < δ/ < ζf ∞B ) coincide for δ/ ) ζf ∞B and then exchange their stabilities. The dot-line curve, dotted line, and continuous lines indicate unstable, saddle, and stable points, respectively.

an invariant 1-D manifold characterized by the following exterior R/ω-Lyapunov numbers ΛωE,1 and ΛRE,1 ω ΛE,1

-1 - ζ ) 1

(39)

By increasing the value of δ/, we find that the point-atinfinity u∞eq,1 moves toward u∞eq,3. For δ/ ) f ∞A, the points-atinfinity u∞eq,1 and u∞eq,3 coincide and a transcritical bifurcation occurs on the boundary ∂S 13 corresponding to an exchange of stability between the stable node and the saddle point. Actually, ∞ for δ/ > f ∞A, the points ueq,3 and u∞eq,1 have exchanged their ∞ stabilities and the point ueq,1 has moved into the physical region, becoming the stable node. Figure 5 shows how the point-at-infinity u∞eq,1 moves on ∂S 13 for increasing values of the parameter δ/ and how its stability changes, from saddle point (dotted curve) to stable node (continuous line) after the bifurcation point (δ/ ) f ∞A, u∞eq,1 ) u∞eq,3). Therefore, for δ/ > f ∞A, the u3-axis ceases to be the generalized slow manifold because there exists a global 1-D slow manifold representing the heteroclinic connection between u∞eq,1 and the equilibrium point ueq ) 0. In the original coordinates, this 1-D global slow manifold is a curve passing through zeq ) 0 (tangent at zeq to the slow eigendirection e1) and tangent, at infinity, to the direction e∞1 associated with the stable node u∞eq,1 ∈ F. Figure 6 shows the 1-D global slow manifold for f ∞A < δ/ < ζf ∞B together with some system orbits highlighting also the structure of the 2-D global slow manifold tangent, at infinity, to the hyperplane span(e∞1 , e∞3 ) (spanned by the directions associated with the stable node u∞eq,1 and the saddle point u∞eq,3) and tangent, at the equilibrium point zeq ) 0, to the hyperplane span(e1, e2) (spanned by the slower eigendirections e1 and e2). The dot-line portion of the 1-D slow manifold indicates the region where the global slow manifold is inverting, as can be

Figure 6. 1- and 2-D global slow manifolds for 1 < ζ < δ/ and f ∞A < δ/ < ζf ∞B . The 2-D slow manifold is tangent, at infinity, to the hyperplane span(e∞1 , e∞3 ) (spanned by the directions associated with the stable node u∞eq,1 and the saddle point u∞eq,3) and tangent, at the equilibrium point zeq ) 0 to the hyperplane span(e1, e2) (spanned by the slower eigendirections e1 and e2). The dot-line portion of the 1-D global slow manifold indicates the region where the manifold is inverting.

also noticed by observing the spatial behavior of a system orbit (broken line) that, starting close to the inverting region, evolves almost parallel to the manifold in the inverting region. The presence of an inverting region along the global 1-D slow manifold is confirmed by the analysis of the quantity ω(2) ν /(2 (1) ωτ ) along the manifold itself (in Figure 7 , it is shown as a function of the z-coordinates, monotonically increasing along the manifold from zeq to infinity). It reveals the presence of two consecutive regions where the manifold is inverting: region (1) (2) (1) a where ω(2) ν /(2ωτ ) < 1 and region b where ων /(2ωτ ) < 0. By further increasing the value of δ/, another transcritical bifurcation occurs for δ/ ) ζf ∞B , involving the two points-at∞ that exchange their stability so that, for infinity u∞eq,2 and ueq,3 ∞ δ/ > ζf B , the point u∞eq,2 ∈ F has become the saddle belonging to the physical region and the unstable node. Figure 5 shows ∞ how the point-at-infinity ueq,2 moves on ∂S 13 for increasing values of the parameter δ/ and how its stability changes, from


(1) Figure 7. |ω(2) ν /(2ωτ )| vs z along the 1-D global slow manifold for 1 < ζ < δ/ and f ∞A < δ/ < ζf ∞B . There are two consecutive regions (1) where the manifold is inverting: region a, where ω(2) ν /(2ωτ ) < 1, and (2) (1) region b, where ων /(2ωτ ) < 0 (dotted line).

Figure 8. 1- and 2-D global slow manifolds for 1 < ζ < δ/ and δ/ > ζf ∞B . The 2-D slow manifold is tangent, at infinity, to the hyperplane span(e∞1 , e∞2 ) and tangent, at the equilibrium point zeq ) 0 to the hyperplane span(e1, e2) (spanned by the slower eigendirections e1 and e2). The dot-line portion of the 1-D global slow manifold indicates the region where the manifold is inverting.

unstable point (dot-line curve) to saddle point (dotted line) after ∞ ∞ ) ueq,3 ). the bifurcation point (δ/ ) ζf ∞B , ueq,2 This second transcritical bifurcation influences solely the structure of 2-D slow manifolds. For f ∞A < δ/ < ζf ∞B , the 2-D slow manifold is a surface tangent, at infinity, to the hyperplane span(e∞1 , e∞3 ). After this second transcritical bifurcation, for δ/ > ζf ∞B , the 2-D slow manifold is a surface tangent, at infinity, to the hyperplane span(e∞1 , e∞2 ) (spanned by the directions associated with the stable node u∞eq,1 and the new saddle point u∞eq,2 ∈ F) and tangent, at the equilibrium point zeq ) 0, to the hyperplane span(e1, e2) (spanned by the slower eigendirections e1 and e2, see Figure 8). This bifurcation analysis of the equilibrium points-at-infinity is supported and confirmed by the analysis of the exterior R R-Lyapunov numbers ΛE,m , for m ) 1 and 2, computed along the 1- and 2-D manifolds defined by the asymptotic directions associated with the equilibrium points-at-infinity. Figure 9A shows the behavior of ΛRE,1 computed along the

J. Phys. Chem. A, Vol. 110, No. 50, 2006 13471 three different asymptotic directions e∞1 , e∞2 , and e∞3 , as a function of the bifurcation parameter δ/. It can be observed ∞ that, for δ/ < f ∞A, i.e., when ueq,3 is a stable node, then ΛRE,1( ∞ e3 ) > 1 and it attains the largest value ΛRE,1(e∞3 ) > ΛRE,1(e∞1 ) > ΛRE,2(e∞2 ). For δ/ > f ∞A, i.e., after the first transcritical bifurcation when u∞eq,1 becomes a stable node, then ΛRE,1(e∞1 ) > 1 and it attains the largest value ΛRE,1(e∞1 ) > ΛRE,1(e∞3 ) > ΛRE,2(e∞2 ). In a similar way, we computed the exterior R-Lyapunov numbers ΛRE,2 (see Figure 9B) along the three different planes span(e∞1 , e∞2 ), span(e∞3 , e∞1 ), and span(e∞3 , e∞2 ) spanned by the three different asymptotic directions, as a function of the parameter δ/. It can be observed that, for δ/ < ζf ∞B , ΛRE,2(e∞3 , e∞1 ) > 1 attains the largest values, and for δ/ > ζf ∞B , i.e., after the second transcritical bifurcation, ΛRE,2(e∞1 , e∞2 ) > 1 attains the largest values. To sum up, we have shown that the first transcritical bifurcation (the exchange of stability between the stable node and the saddle point-at-infinity) influences the structure of the 1-D slow manifold, and the second transcritical bifurcation (the exchange of stability between the saddle point and the unstable node) influences the structure of the 2-D slow manifold. This analysis has been performed by assuming 1 < ζ < δ/. The complete bifurcation diagram for the 3-D Semenov model is presented and discussed in the next section. 4.3. Bifurcation Diagram. Figure 10 shows the locus of bifurcation points in the parameter space ζ-δ/ for γA < γB. The general features of the diagram are not dependent upon the value of the parameters and h. Continuous thick lines indicate the occurrence of the first transcritical bifurcation, i.e., the exchange of stability between the stable node and the saddle point. Broken lines indicate the occurrence of the second transcritical bifurcation, i.e., the exchange of stability between the saddle point and the unstable node, that influences the structure of the 2-D slow manifold. In the parameter space ζ-δ/, it is possible to identify 17 different regions corresponding to different features of the 1and 2-D slow manifolds, as reviewed in Tables 3 and 4. Gray regions (A, C, E-G, J-L, O-Q) are characterized by the existence of a global 1-D slow manifold. For example, Figures 6 and 8 show the 1-D global (and inverting) slow manifold for (ζ, δ/) ∈ O and P, respectively. In the three gray regions A, G, and L, the z-axis plays the role of the global noninverting 1-D slow manifold, and in the white regions B, D, H, I, M, and N, the z-axis represents the generalized 1-D slow manifold. All the regions, except E, D, I, and N, are characterized by the existence of a 2-D global slow manifold (see, for example, Figure 6 for region O and Figure 7 for region P). When the 1- or 2-D slow manifold is a generalized manifold, the local behavior of the system close to the equilibrium point

R Figure 9. ΛE,m vs δ/, for m ) 1 and 2, computed along the 1- and 2-D manifolds identified by the asymptotic directions e∞1 , e∞2 , and e∞3 associated with the equilibrium points-at-infinity.


Giona et al.

TABLE 3: 1- and 2-D Slow Manifolds for Different Values of the Two Parameters ζ and δ/Ea parameters A A B B C C D D E E F F G G H H I I

δ/ < ζ < 1 δ/ < ζf ∞B < f ∞A ζ < δ/ < 1 δ/ < ζf ∞B < f ∞A ζ < δ/ < 1 ζf ∞B < δ/ < f ∞A ζ < 1 < δ/ δ/ < ζf ∞B < f ∞A ζ < 1 < δ/ ζf ∞B < δ/ < f ∞A ζ < 1 < δ/ ζfB < f ∞A < δ/ δ/ < ζ < 1 δ/ < f ∞A < ζf ∞B ζ < δ/ < 1 δ/ < f ∞A < ζf ∞B ζ < 1 < δ/ δ/ < f ∞A < ζf ∞B

local eigendirections

1- and 2-D slow manifolds W (1) ) span{e3} ) span{e∞3 } ) (0, 0, z) W (2) ) span{e3, e2} ) span{e∞3 , e∞2 } ) (0, y, z) W (1) ) span{e∞3 } ) (0, 0, z) W (2) ) span{e2, e3} ) span{e∞3 , e∞2 } ) (0, y, z) W (1) ) [e2 r e∞2 ] W (2) ) span{e2, e3} ) span{e∞2 , e∞3 } ) (0, y, z) W (1) ) span{e∞3 } ) (0, 0, z) W (2) ) span{e∞3 , e∞2 } ) (0, y, z) W (1) ) [e2 r e2∞] W (2) ) span{e∞2 , e∞3 } ) (0, y, z) W (1) ) [e2 r e∞2 ] W (2) ) [{e2, e1} r {e∞2 , e∞1 }] W (1) ) span{e3} ) span{e∞3 } ) (0, 0, z) W (2) ) [{e3, e2} r {e∞3 , e∞1 }] W (1) ) span{e∞3 } ) (0, 0, z) W (2) ) [{e2, e3} r {e∞3 , e∞1 }] W (1) ) span{e∞3 } ) (0, 0, z) W (2) ) [{e2, e3} r {e∞3 , e∞1 }]

e3, e2, e1 e∞3 , e∞2 , e∞1 e2, e3, e1 e∞3 , e∞2 , e∞1 e2, e3, e1 e∞2 , e∞3 , e∞1 e2, e1, e3 e∞3 , e∞2 , e∞1 e2, e1, e3 e∞2 , e∞3 , e∞1 e2, e1, e3 e∞2 , e∞1 , e∞3 e3, e2, e1 e∞3 , e∞1 , e∞2 e2, e3, e1 e∞3 , e∞1 , e∞2 e2, e1, e3 e∞3 , e∞1 , e∞2

global global generalized global global global generalized generalized global generalized global global global global generalized global generalized generalized

Refer to Figure 10 for the identification of the different regions, from A to I. The local eigendirections {ei} (at zeq) and {e∞h } (at infinity) are ordered from the slowest to the fastest one. The symbol [ei r e∞h ] indicates that the 1-D slow manifold is tangent at zeq to the eigendirection ei and ∞ tangent, at infinity, to the direction e∞h associated with the point-at-infinity ueq,h . The symbol [{ei, ej} r {e∞h ,e∞k }] indicates that the 2-D slow ∞ ∞ ∞ ∞ manifold is tangent, at infinity, to the plane span (eh , ek ) (spanned by the directions associated with the points-at-infinity ueq,h and ueq,k ) and tangent, at the equilibrium point zeq ) 0 to the hyperplane span(ei, ej) (spanned by the eigendirections ei and ej). a

TABLE 4: 1- and 2-D Slow Manifolds for Different Values of the Two Parameters ζ and δ/Ea parameters ζ < 1 < δ/ f ∞A < δ/ < ζf ∞B ζ < 1 < δ/ fA < ζf ∞B < δ/ δ/ < 1 < ζ δ/ < f ∞A < ζ f ∞B 1 < δ/ < ζ δ/ < f ∞A < ζ f ∞B 1 < ζ < δ/ δ/ < f ∞A < ζf ∞B 1 < ζ < δ/ f ∞A < δ/ < ζf ∞B 1 < ζ < δ/ f ∞A < ζ f ∞B < δ/ 1 < δ/ < ζ f ∞A < δ/ < ζ f ∞B

J J K K L L M M N N O O P P Q Q a

local eigendirections

1- and 2-D slow manifolds e∞1 ]

W (1) ) [e2 r W (2) ) [{e2,e1} r {e∞1 , e∞3 }] W (1) ) [e2 r e∞1 ] W (2) ) [{e2,e1} r {e∞1 , e∞2 }] W (1) ) span{e3} ) span{e∞3 } ) (0, 0, z) W (2) ) [{e3,e1} r {e∞3 ,e∞1 }] W (1) ) span{e∞3 } ) (0, 0, z) W (2) ) [{e1, e3} r {e∞1 , e∞3 }] W (1) ) span{e∞3 } ) (0, 0, z) W (2) ) [{e1, e3} r {e∞3 , e∞1 }] W (1) ) [e1 r e∞1 ] W (2) ) [{e1, e2} r {e∞1 , e∞3 }] W (1) ) [e1 r e1∞] W (2) ) [{e1, e2} r {e∞1 , e∞2 }] W (1) ) [e1 r e∞1 ] W (2) ) [{e1, e3} r {e∞1 , e∞3 }]

e2, e1, e3 e∞1 , e∞3 , e∞2 e2, e1, e3 e∞1 , e∞2 , e∞3 e3, e1, e2 e∞3 , e∞1 , e∞2 e1, e3, e2 e∞3 , e∞1 , e∞2 e1, e2, e3 e∞3 , e∞1 , e∞2 e1, e2, e3 e∞1 , e∞3 , e∞2 e1, e2, e3 e∞1 , e∞2 , e∞3 e1, e3, e2 e∞1 , e∞3 , e∞2

global global global global global global generalized global generalized generalized global global global global global global

Refer to Figure 10 for the identification of the different regions, from J to Q.

zeq is controlled by the occurrence of a transient (or HartmanGrobman) 1- or 2-D slow manifold.8 This phenomenon can be appreciated by observing the phase-space portrait of the 3-D Semenov model corresponding to parameter values ζ - δ/ falling in region D (see Figure 11). ∞ In region D, the points-at-infinity u∞eq,3 and ueq,2 are the stable node and the saddle point, respectively, so that the invariant z-axis and the invariant plane (0, y, z) are the generalized 1- and 2-D slow manifolds. Actually, the z-axis is inverting for z ∈ [0, zmax], where zmax is determined by enforcing the condition

ω(2) ν

|

0,0zmax 2ω(1) τ

)

fA(zmax) + ζfB(zmax) )1 2δ/

(40)

Phase-space orbits starting far from the equilibrium point rapidly relax onto the 2-D generalized manifold (0, y, z), and

orbits starting close to the equilibrium point settle down onto the transient manifold represented by a finite portion of the plane span(e2, e1), e2 and e1 being the two slower eigendirections of the equilibrium point zeq. Similarly, a finite portion of the linear manifold span(e2) is the 1-D Hartman-Grobman manifold. 5. Concluding Remarks This Article has completed the analysis developed in ref 8 on the geometric characterization of slow invariant manifolds extending it to high dimensional systems. The use of tools and methods deriving from exterior algebra provides a simple way to characterize stretching dynamics and the relative contraction of normal measure elements compared to tangential ones along generic m-dimensional invariant manifolds. This Article has thoroughly examined the bifurcational properties of a 3-D prototypical combustion system, showing how the local bifurcations of the equilibrium points-at-infinity


Figure 10. Locus of bifurcation points in the parameter space ζ-δ/ for γA < γB. The general features of the diagram are independent of the value of the parameters and h. Continuous thick lines: exchange of stability between the stable and the saddle points (it affects 1-D slow manifolds). Broken lines: exchange of stability between the saddle and the unstable points (it affects 2-D slow manifolds). The 17 different regions, from A to Q, are characterized by different features of the 1and 2-D slow manifolds, as reviewed in Tables 3 and 4. Gray (white) regions are characterized by the existence of a global (generalized) 1-D slow manifold. Regions D, E, I, and N are characterized by the existence of a generalized 2-D slow manifold.

J. Phys. Chem. A, Vol. 110, No. 50, 2006 13473 of the elementary concepts of the cross-vector product in R3 and determinant algebra.10-12 Let E be the tangent space at a point z of the n-dimensional phase space Rn; E is an n-dimensional vector space. It is possible to introduce a family of new vector spaces, E∧,1, E∧,2, ..., E∧,m, the elements of which are the measure elements of dimension 1, 2, ..., n at z. By definition, E∧,1 ) E, and E∧,p is a vector space composed of elements (referred to as skew-symmetric p-forms or simply p-dimensional measure elements) defined starting from p vectors v(1), ..., v(p) of E, composed through an exterior (or wedge) product ∧:

ap ∈ E∧,p

ap ) v(1) ∧ v(2) ∧ ... ∧ v(p)

v(R) ∈ E (A.1)

The wedge product of vectors belonging to E satisfies the following conditions, which define it uniquely. (1) It is multilinear; i.e., it is linear in all of its p entries. This means that if c1 and c2 are two real constants

(c1v(1) + c2w(1)) ∧ v(2) ∧ ... ∧ v(p) ) c1 (v(1) ∧ v(2) ∧ ... ∧ v(p)) + c2(w(1) ∧ v(2) ∧ ... ∧ v(p)) (A.2) Analogous relationships hold for the linearity referred to the second, third, and pth entry. (2) It is skew-symmetric, i.e.,

v(1) ∧ v(2) ∧ ... ∧ v(p) ) 0 if v(i) ) v(j) for i * j

Figure 11. 1- and 2-D generalized and transient slow manifolds for ζ < 1 < δ/ and δ/ < ζf ∞B < f ∞A (region D). 1-D generalized slow manifold: z-axis. 1-D transient slow manifold: finite portion of the linear manifold span(e2). 2-D generalized slow manifold: plane (0, y, z). 2-D transient slow manifold: finite portion of the plane span(e2, e1).

modify the occurrence and the dynamic properties of invariant slow 1- and 2-D manifolds. A complete bifurcational analysis in the parameter space has been obtained and the results explained by means of the stretching properties experienced by normal and tangential measure elements along the invariant manifolds. In the 3-D Semenov system, the bifurcations of the points-at-infinity, which modify the structure of the slow 2- and 3-D invariant manifolds, are much more complex than those in the classical 2-D Semenov system studied by Creta et al.8 The bifurcational analysis developed in this Article is, to the best of our knowledge, the first complete bifurcational characterization of a 3-D chemical system as it regards the existence and the properties of the global/generalized slow manifolds. The analysis developed for the 3-D Semenov model can be extended in principle to generic reactive schemes of practical interest. Indeed, albeit the apparent formal complexity, the use of exterior Lyapunov-type numbers is sufficiently simple from the computational point of view to be applied without any major problem to higher dimensional kinetic models. This will be developed in future works with particular emphasis on combustion and on biochemical reaction networks. Appendix A Exterior algebra is the algebra of the exterior forms (or measure elements) and is the natural algebraic generalization

(A.3)

The latter condition implies that v(1) ∧ ... ∧ v(p) changes sign whenever two vectors v(i), v(j) permute. Equations 42 and 43 are the natural conditions arising from the definition of p-dimensional measure elements starting from p vectors. For example, if the vectors v(1), ..., v(p) are linearly dependent, it follows from eqs 42 and 43 that the resulting p-form is identically equal to zero. n be an orthonormal basis for E. A basis for E∧,p is Let {ei}i)1 given by the family of wedge products

ei1,...,ip ) ei1 ∧ ei2 ∧ ... ∧ eip with i1 < i2 < ...< ip

(A.4)

and where ik ) 1, ..., n with k ) 1, ..., p. The dimension of E∧,p equals the number of combinations of n elements of class p, i.e., dim(E∧,p) ) n!/(p!(n - p)!). For n ) 3, for example, E∧,1 ) E, dim(E∧,1) ) 3 and the canonical basis for E∧,1 is trivially {e1, e2, e3}. Similarly, dim(E∧,2) ) 3 and the canonical basis for E∧,2 is {e1 ∧ e2, e1 ∧ e3, e2 ∧ e3}, whereas dim(E∧,3) ) 1 with the only element {e1 ∧ e2 ∧ e3} in its basis. A generic measure element belonging to E∧,p can be expressed as a linear combination of elements of this basis

ap )

∑ ai ,...,i ei ,...,i i ... > i2 > i1. From the definition, eqs 42 and 43, it follows that, if ap ) v(1) ∧ ... ∧ v(p), the coefficient ai1,...,ip of ap with respect to the natural basis {ei1,...,ip} of E∧,p is given by the determinant

|

|

13474 J. Phys. Chem. A, Vol. 110, No. 50, 2006 (1) Vi(1) Vi(1) ... ... Vip 1 2

ai1,...,ip )

(2) Vi(2) Vi(2) ... ... Vip 2 2 ... ... ... ... ... (p) Vi(p) Vi(p) ... ... Vip 1 2

Giona et al.

(A.6)

p The orthonormal nature of the basis {ei}i)1 for E induces the orthonormal nature of {ei1,...,ip} as a basis for E∧,p. This means that an inner (scalar) product (‚, ‚)p can be defined for E∧,p with respect to this representation as

(ap, bp)p )

∑ ai ,...,i bi ,...,i i

Slow Manifold Structure in Explosive Kinetics. 2 ... - ACS Publications

Recommend Documents