Quasi Equilibrium Grid Algorithm: geometric construction for model reduction
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r Quasi Equilibrium Grid Algorithm: geometricconstruction for model reduction
Eliodoro Chiavazzo, Iliya V. Karlin
Aerothermochemistry and Combustion Systems Laboratory (LAV), ETHZCH-8092 Zurich, Switzerland.
Abstract
The Method of Invariant Grid (MIG) is an iterative procedure for model reductionin chemical kinetics which is based on the notion of Slow Invariant Manifold (SIM)[1]-[4]. Important role, in that method, is played by the initial grid which, once re-fined, gives a description of the invariant manifold: the invariant grid. A convenientway to get a first approximation of the SIM is given by the Spectral Quasi Equi-librium Manifold (SQEM) [1]-[2]. In the present paper, a flexible numerical methodto construct the discrete analog of a Quasi Equilibrium Manifold, in any dimen-sion, is presented. That object is named Quasi Equilibrium Grid (QEG), while theprocedure Quasi Equilibrium Grid Algorithm. Extensions of the QEM notion arealso suggested. The QEG is a numerical tool which can be used to find a grid-basedapproximation for the locus of minima of a convex function under some linear con-straints. The method is validated by construction of one and two-dimensional gridsfor model hydrogen oxidation reaction.
Key words:
Chemical kinetics, model reduction, invariant manifold, entropy,nonlinear dynamics, Lagrange multipliers method, variational problem.
Relaxation of complex systems is often characterized by a fast dynamics dur-ing a short initial stage, while the remaining period lasts much longer andit evolves along low-dimensional surfaces in the phase space known as SlowInvariant Manifolds (SIM). In that scenario, a simplified macroscopic descrip-tion of a complex system can be attained by extracting only the slow dy-namics and neglecting the fast one. For this reason, much effort was spent
Email addresses: [email protected] (Eliodoro Chiavazzo), [email protected] (Iliya V. Karlin).
Preprint submitted to Elsevier 30 October 2018 o develop model reduction methods (the Method of Invariant Grid (MIG)[1,2,3,4], the Intrinsic Low Dimensional Manifold method (ILDM) [9,10], theComputational Singular Perturbation method (CSP) [11,12,13], etc.) based onthe notion of SIM. The introduction of a convex Lyapunov function G , when-ever the complex system is supported by such a function, also proves to bevery helpful in model reduction [4,8]. Indeed, it was shown that, through a G function, good approximations of the SIM can be found (e.g. by constructingthe Spectral Quasi Equilibrium Manifold - SQEM - or the Symmetric En-tropic Intrinsic Low Dimensional Manifold - SEILDM - [1,2]) and refined bysome efficient MIG iterations. Moreover, it has been shown that the notionof QEM is also very useful in different fields. For example, it was used in theimplementation of Lattice Boltzmann schemes [6,7]. Construction of a QEMis analytically possible by using the Lagrange multipliers method. However,its implementation becomes too complicated as soon as the number of vari-ables of the problem becomes large: efficient methods, for constructing largedimensional QEM, are still missing. Therefore, in the present paper the no-tion of Quasi Equilibrium Grid (QEG) will be introduced, as a discrete analogof QEM, and a constructive algorithm, applicable in any dimension, will bedeveloped. The procedure suggested proves to be a very flexible tool, so itis possible to get some other SIM approximations, all based on the previousalgorithm, even more accurate than the QEG itself. The paper is organized as follows. In Section 3, some basic notions are out-lined: in particular, the general equations of dissipative reaction kinetics arereviewed, in the notations which are used throughout the paper. At the end ofthat Section, the Method of Invariant Grid (MIG) and the notion of thermo-dynamic projector are briefly discussed (Section 3.2). In Section 4, the QEMdefinition and its geometrical interpretation is given, while in Section 5 the1D Quasi Equilibrium Grid Algorithm is presented. That algorithm is also il-lustrated, by means of an example, in Section 6. The 1D Algorithm extensionto multi-dimensional grids is developed in Section 7. In particular, two pos-sible extension strategies are analyzed: the straightforward extension (Section7.1) and, by following the general idea given in [3], the flag extension (Sec-tion 7.2). Here, it is also shown how the flexibility of the flag extension allowsto get SIM approximation which is better than the Spectral-QEG (Section7.3): the notions of
Guided-QEG and
Symmetric Entropic Guided-QEG areintroduced. An illustrative example, in Section 8, shows how those differentextension techniques work in practice. In order to find out how accurate istheir SIM description, they are also compared on the base of the invariancedefect (Section 8.3). Finally, results are discussed in Section 9.2
Theoretical background
In a closed system with n chemical species A , ..., A n , participating in a com-plex reaction, a generic reversible reaction step can be written as a stoichio-metric equation: α s A + ... + α sn A n ⇋ β s A + ... + β sn A n , (1)where s is the reaction index, s = 1 , ..., r ( r steps in total), and the integers α si and β si are stoichiometric coefficients of the step s . For each reaction step, wecan introduce n -component vectors α s and β s , with components α si and β si ,and the stoichiometric vector γ s = β s - α s . For every A i the extensive variable N i describes the number of particles of that species. If V is the volume, thenthe concentration of A i is c i = N i /V . Dynamics of the species concentrationaccording to the stoichiometric mechanism (1) reads:˙ N = V J ( c ) , J ( c ) = P rs =1 γ s W s ( c ) , (2)where dot denotes the time derivative and W s ( c ) is the reaction rate functionof the step s . In particular, the polynomial form of the reaction rate functionis provided by the mass action law : W s ( c ) = W + s ( c ) − W − s ( c ) = k + s ( T ) n Y i =1 c α i i − k − s ( T ) n Y i =1 c β i i , (3)where k + s ( T ) and k − s ( T ) are the constants of the direct and of the inversereactions rates of the step s respectively. The most popular form of theirdependence is given by the Arrhenius equation: k ± s ( T ) = a ± s T b ± s exp( S ± s /k B ) exp( − H ± s /k B T ) . In the latter equation, a ± s , b ± s are constants and H ± s , S ± s activation enthalpiesand entropies respectively. The rate constants are not independent. Indeed,the principle of detail balance gives a relation between these quantities: W + s ( c eq ) = W − s ( c eq ) , ∀ s = 1 , ..., r, (4)where the positive vector c eq ( T ) is the equilibrium of the system (2). In orderto obtain a closed system of equations, one should supply an equation forthe volume V . For an isolated system the extra-equations are U , V = const (where U is the internal energy), for an isochoric isothermal system we get3 , T = const , and so forth. For example, equation (2) in the latter case simplytakes the form: ˙ c = r X s =1 γ s W s ( c ) = J ( c ) . (5)Finally, also other linear constraints, related to the conservation of atoms,must be considered. In general such conservation laws can have the followingform: Dc = const , (6)where l fixed and linearly independent vectors d i are the rows of the l × n matrix D , and const is a constant vector. In this section, we give an outline of the MIG for chemical kinetics. For detailssee Refs. [1,2,3,4,8].
If we turn our attention to perfectly stirred closed chemically active mixtures,then dissipative properties of such systems can be characterized with a ther-modynamic potential which is the Lyapunov function of equation (2). Thatfunction implements Second Law of thermodynamics: it means that during theconcentrations evolution in time, from the initial condition to the equilibriumstate, the Lyapunov function must decrease monotonically. Therefore if G ( c )is the Lyapunov function, c eq (equilibrium state) is its point of global mini-mum in the phase space. A simple example of a function G is given by the freeenergy of ideal gas in a constant volume and under a constant temperature: G = n X i =1 c i [ln( c i /c eqi ) − . (7)When G is known, also its gradient ∇ G and the matrix of second derivatives H = k ∂ G/∂c i ∂c j k can be evaluated, so that it is possible to introduce thethermodynamic scalar product as follows: h x , y i = ( x , Hy ) , (8)where the notation ( , ) is the usual Euclidean scalar product.4 .2.2 The invariance condition Let us consider Ω as a manifold of a reduced description. The invariancerequirement reads: c (0) ∈ Ω ⇒ c ( t ) ∈ Ω , ∀ t ≥ . (9)Let P be a projector on the tangent bundle of the manifold Ω . The manifold Ω is invariant with respect to the system (2) if and only if the followinginvariance equation (IE) holds:[1 − P ] J ( c ) = 0 , ∀ c ∈ Ω . (10)When the manifold is not invariant, it is not able to satisfy the invariancecondition so that: ∃ c : ∆ = [1 − P ] J ( c ) = 0 , (11)where ∆ is the defect of invariance. One way to find the SIM is to solve theIE iteratively starting from an appropriate initial manifold. Let us now discuss further the projector appearing in the invariance equation.It is an operator which for each point c ∈ Ω projects the vectors J ( c ) onto thetangent subspace of the manifold producing, in this way, the induced vectorfield P J ( c ). In general, condition (10) does not require any special constraintfor the projector P . However, the thermodynamic properties of the kineticequations (2) define the projector unambiguously [1,4,8]. To this end, let usdefine a differential of G , that is linear functional: DG ( x ) = ( ∇ G ( c ) , x ) . (12)A special class of projectors is the thermodynamic one. If a projector belongsto this class then the induced vector field respects the dissipation inequality: DG ( P J ) ≤ , ∀ c ∈ Ω . (13)It has been shown that a projector P respects the (13) if and only if [8]:ker P ⊆ ker DG, ∀ c ∈ Ω , (14)where ker denotes the null-space of an operator. It is clear now that if onewants to solve equation (10), then a projector must be specified. Here weremind the way to construct the thermodynamic projector which will be usedin MIG procedure [1]. This projector depends on the concentration point c and on the tangent space to the manifold Ω .5e are looking for a grid approximation of a q -dimensional SIM. Let G bea discrete subset of a q -dimensional parameter space R q and let F | G be amapping of G into the concentration space. If we select an approximationprocedure to restore the smooth map F from the discrete map F | G (we needa very small part of F , derivatives of F in the grid points only), then thederivatives f i = ∂F/∂y i are available, and for each grid point the tangentspace is: T y = Lin { f i } , i = 1 , ..., n. (15)We assume that one of points y ∈ G maps into the equilibrium, and inother points intersection of the manifold with G levels is transversal (i.e.( DG ) F ( y ) ( x ) = 0 for some x ∈ T y ). Let us consider the subspace T y =( T y ∩ ker DG ). In order to define the thermodynamic projector, it is required,if T y = T y , to introduce the vector e y which satisfies the following conditions: e y ∈ T y , h e y , x i = 0 , ∀ x ∈ T y ,DG ( e y ) = 1 . Let P be the orthogonal projector on T y with respect to the entropic scalarproduct (8), then the thermodynamic projection of a vector x is defined as: T y = T y ⇒ P x = P x + e y DG ( x ) T y = T y ⇒ P x = P x . (16) When MIG method is applied, not a manifold is searched as a solution, but aset of concentration points whose defect of invariance is sufficiently small: let Ω denote that solution (invariant grid). MIG is an iterative procedure: this meansthat, at the beginning, only an initial approximation Ω of Ω is available. Ingeneral, Ω does not respect the invariance condition (10) satisfactorily so the(11) holds: for this reason the position of c ∈ Ω must be changed. We canthink to correct its position and get a new point ( c + δ c ) with a lower defectof invariance ∆ = [1 − P ] J ( c + δ c ). If the initial node is “not far” from theinvariant manifold, a reasonable way to get the node correction δ c is to solvethe linearized invariance equation where the vector field J is expanded to thefirst order and the projector P to the zeroth order:[1 − P ( c )][ J ( c ) + L ( c ) δ c ] = 0 . (17) L is the matrix of first derivatives of J (Jacobian matrix). The Newton methodwith incomplete linearization consists of the equation (17) supplied by the6xtra condition [8]: P δ c = 0 . (18)The additional condition (18) and the atoms balances (6) automatically can betaken into account choosing a basis { b i } in the subspace S = (ker P ∩ ker D ).Let h = dim ( S ), then the correction can be cast in the form δ c = P hi =1 δ i b i ,so that the linearized invariance equation (17) becomes the linear algebraicsystem in terms of δ i : P hi =1 δ i ((1 − P ) Lb i , b k ) = − ((1 − P ) J , b k ) , k = 1 , ..., h. (19) Remark.
Here the usual scalar product ( , ) was used to get the components ofthe left-hand side of (17) in the basis vectors { b i } . Nevertheless, a differentscalar product can be also used without a loss of generality.In the case of the thermodynamic projector, it proves convenient to choosethe basis { b i } orthonormal with respect to the entropic scalar product (8) andwrite the (19) as: P hi =1 δ i h (1 − P ) Lb i , b k i = − h (1 − P ) J , b k i , k = 1 , ..., h. (20)The projector (16) is “almost” h , i − orthogonal ( h im P , ker P i ∼ = 0) close tothe SIM. Because of that special feature, equation (20) can be approximatedand simplified as follows: P hi =1 δ i h Lb i , b k i = − h J , b k i , k = 1 , ..., h. (21)Note that, in general, an approximation carried out by eq. (21) leaves a residualdefect (11) in the grid nodes which cannot be completely annihilated. There-fore, when a higher accuracy in the SIM description is required, equation (19)is recommended. An alternative approach to solve eq. (17) is the relaxation method . Accordingto that method the correction is written as c = c + τ ( c )∆( c ), and the quantity τ ( c ) is obtained from the condition: h ∆ , [1 − P ][ J + τ ( c ) L ∆] i = 0 , and solving with respect to τ : τ ( c ) = − h ∆ , ∆ ih ∆ , L ∆ i . (22)Equation (22) shows that the relaxation method is explicit, but as it adjuststhe node position acting only along the direction of the defect ∆, typically we7xpect it to be less efficient in comparison with the Newton method. On theother hand, this method is particularly easy to implement. Any iterative procedure needs to be supplied by a first approximation. Sinceit plays an important role for both the convergence and efficiency, that ap-proximation must be chosen very carefully. It was shown that a reasonableway, for initializing the MIG, is to construct the Quasi Equilibrium Manifold(QEM) [1,2].
Solution trajectories in the phase-space must obey the set of ODE equations(5). Moreover, all trajectories also satisfy a subset of linear equations (6)which represent the atom conservation. Among all the concentration points,that fulfill the latter constraints, we can choose those points which minimizethe Lyapunov function G of the system we are dealing with. Such pointslie on a manifold that is called Quasi Equilibrium Manifold (QEM). Let ussuppose that some steps of a complex reaction are faster than some others.Since the Lyapunov function G must decrease during the fast dynamics then,when the fast motion is exhausted, the G value is expected to be the minimumon that fast hyperplane. In such a situation, a QEM attempts to achieve amotion decomposition into fast - toward the QEM - and slow - along theQEM. If the invariant manifold exists, the QEM can be taken as a reasonableapproximation of it. In order to be more specific, let a chemical system have n reactive species. The degrees of freedom of that system are ( n − l ) becauseof the atom balances (6). If q < ( n − l ) is the dimension of the QEM, thenthe macroscopic variables for its description are ξ , ..., ξ q so that: ( m , c ) = ξ , ..., ( m q , c ) = ξ q . Here, the n -dimensional vectors m i are related to thehypothetic fast directions. From a mathematical standpoint, the solution of avariational problem: G → min ( m i , c ) = ξ i , ∀ i = 1 , ..., q Dc = const (23)represents the QEM corresponding to the vector set { m i } . We want to stressthe geometry behind (23) because it will be extensively exploited in the fol-lowing. The geometric interpretation of a QEM is illustrated in Fig. 1 for a8 A i C A j C A i (a) (b) Equilibrium Equilibrium t t Fig. 1. Quasi Equilibrium Manifold: the geometrical interpretation. Two differentQE-manifolds (bold lines in (a) and (b)) corresponding to two different linear con-straints subsets in the problem (23). c A i , c A j ). Let us consider the points where G levelcurves (convex curves in Figures 1 (a)-(b)) are cut by the QE-manifolds (boldcurves): in those points the inclination of the tangent to the G -level curvesis constant. Different QEM can be obtained by choosing different vector sets { m i } . A special choice is done when m , ..., m q are the q left eigenvectors ofthe Jacobi matrix L ( c eq ) corresponding to the q smallest absolute eigenval-ues. In that particular case, solution of (23) has its own name: Spectral QuasiEquilibrium Manifold (SQEM) [1,2].
The minimization problem (23) can be, in principle, solved by the method ofLagrange multipliers. However, it is also well known that, when the number ofconstraints and variables increases, then that method becomes prohibitivelycomplicated to implement. Since the number of species and elementary reac-tion steps is usually quite high, the Lagrange multipliers method is not suitablefor most of the practical cases. For this reason, in the sequel a new procedureto overcome that issue is presented. An algorithm (
Quasi Equilibrium GridAlgorithm -QEGA), which can be easily implemented in order to get a dis-crete analog of a QEM in any dimension, is developed. This is achieved byinvestigating further the QEM geometrical construction.9 ig. 2. Quasi Equilibrium Grid: the basic idea.
Let us consider a one-dimensional quasi-equilibrium manifold. Let us assumethat the node c belongs to that manifold. One may now imagine to look fora new node c which still lies on the quasi equilibrium manifold. In general,the node c can be obtained from c by adding a shift ˆ δ c : c = c + ˆ δ c .That idea is applicable whenever a QEM-node c n is known and a new one c n +1 must be found (see Fig. 2): c n +1 = c n + ˆ δ c n . (24)First of all, any node c has to fulfill the atom balances (6). Let { ρ i } be a basisin the null space of matrix D . A convenient way to take automatically intoaccount the conditions (6) is to express any shift ˆ δ c n as a linear combinationof vectors ρ i : ˆ δ c n = X zi =1 µ i ρ i , (25)where z = n − l is the dimension of the basis { ρ i } . By referring to Fig.2, let us now discuss further the tangent space T to the G level surface inany quasi-equilibrium point c n +1 . The space T geometrically represents thelinear constraint of the problem (23). Therefore, any point c of T satisfies thatconstraint, but only c n +1 minimizes G function. The line – l passing from c n +1 and c has the parametric form c = ϕ ˜ t + c n +1 , where ˜ t is a vector of T spanning– l while ϕ is a parameter. In general, the linear constraints of the problem (23)can be also written as: ( m , c ) = ϕ ( m , ˜ t ) + ( m , c n +1 ) ⇒ ( m , ˜ t ) = 0 , ∀ ˜ t ( d i , c ) = ϕ ( d i , ˜ t ) + ( d i , c n +1 ) ⇒ ( d i , ˜ t ) = 0 , ∀ ˜ t (26)where m and d i are the reduced variable vector ( q = 1) and the generic row ofmatrix D , respectively. The vector ˜ t , that respects (26), can always be written10s a linear combination of some vectors t j , where { t j } denotes a basis in thenull space of that matrix E , whose first row is given by m and the rest by therows of D : E = mD . (27)Note that the dimension of { t j } is z −
1. By looking at Fig. 2, the quasi-equilibrium requirement simply becomes the orthogonality condition:( ∇ G ( c n +1 ) , ˜ t ) = 0 , ∀ ˜ t ∈ T (28)which also means: ( ∇ G ( c n +1 ) , t j ) = 0 , ∀ j = 1 , ..., z − . (29)The quasi-equilibrium grid algorithm is based on the equation system (29)and two more assumptions. First of all, we suppose that the known node c n is close to the QEM, although it does not necessarily belong to the QEM.Secondly, let the vector ˆ δ c n be small enough, so that the gradient ∇ G ( c n +1 )can be approximated to the first order: ∇ G ( c n +1 ) ∼ = ∇ G ( c n ) + H ( c n )ˆ δ c n , (30)where H ( c n ) = h ∂ G∂c i ∂c j i denotes again the matrix of second derivatives of thefunction G evaluated at the point c n . By substituting equations (30) and (25)in (29), we obtain: P zi =1 ( t j , H ( c n ) ρ i ) µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − . (31)By using the entropic scalar product (8), equations (31) can be cast into theform: P zi =1 h t j , ρ i i µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − . (32)Both the matrix H and the gradient ∇ G are calculated at the known node c n . Note that the right-hand side of (32) vanishes if the node c n belongs tothe QEM. The node collection, subsequently evaluated through (32), will becalled a Quasi Equilibrium Grid (QEG).
Note, however, that the system (32) is not closed ( z unknowns µ i , but z − P zi =1 h t j , ρ i i µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) = ε (33)where ε is a given number and (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) represents the Euclidean norm of thevector ˆ δ c n . The smaller ε is chosen, the more accurate the expression (30)gets. As it will be shown later on, for small ε the Quasi Equilibrium Gridlies very close to the correspondent Quasi Equilibrium Manifold. The extracondition makes (33) a non-linear algebraic system. A way to solve it will benow discussed. The idea is to find the general solution of the linear system(32), and then to choose the one which also fulfills the non linear condition in(33). Let the basis { ρ i } be orthonormal (in the Euclidean sense). That is notcrucial, but it proves to be convenient in the following analysis; indeed thenon-linear system (33) now is cast as follows: P zi =1 h t j , ρ i i µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − P zi =1 µ i = ε . (34)The general solution of (32) can always be written as: µ ... µ z = w ν ... ν z + p ... p z , (35)where w is a free parameter, while ν = [ ν , ..., ν z ] T and p = [ p , ..., p z ] T arethe solution of the homogeneous problem and a special solution of (32), re-spectively. Without any restriction, we assume ( ν , ν ) = 1. Once ν and p areknown, the non linear equation of (34) can be written, in terms of w , as: w + 2( ν , p ) w + ( p , p ) − ε = 0 . (36)If the solvability condition is satisfied,( ν , p ) − ( p , p ) + ε > , (37)then the two real valued solutions of (36) ( w I , w II ), upon substitution into(35), give two possible sets [ µ , ..., µ z ]. Therefore, by using the (24) and (25),two new nodes c In +1 , c IIn +1 (both close to the quasi equilibrium curve) can beevaluated from the previous one c n (see Fig. 3). A criterion, able to choosebetween those two solutions, depends on the phase-space zone where the gridneeds to be constructed. This idea will be clarified in the following example.12 ig. 3. Two solutions for the 1D QEG algorithm. Usually, the equilibrium point is supposed to be a good starting node for theQEG procedure: c = c eq . Remark.
The QEG-equations (32) can be generalized as follows: P zi =1 h t j , ρ i i µ i = − η ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − , (38)where η is a parameter 0 ≤ η ≤
1. When η = 1, (32) is recovered. On the otherhand, if the QEG-nodes are close to the QEM, then the non-homogeneousterms can be neglected (they vanish on the QEM). Therefore, a reasonableapproximation of the system (32) is given when η = 0. In the latter case,the solvability condition (37) is fulfilled. If η = 1 and (37) does not hold,that parameter can be chosen in such a way that the solvability condition issatisfied. In the example below, solvability condition (37) is always satisfiedand we use (33). In this section, an example will be considered in order to illustrate howthe algorithm, described in the previous section, works for finding a one-dimensional SQE-grid. That grid will be compared with the relative spec-tral quasi-equilibrium manifold, too. Let us consider the following four-step13 c B c C Equilibrium
Fig. 4. Reaction (39): some trajectories projected into the phase-subspace ( c C , c B ). three-component reaction (kindly suggested by A.N. Gorban): .A ↔ B, k +1 = 1 , .B ↔ C, k +2 = 1 , .C ↔ A, k +3 = 1 , .A + B ↔ C, k +4 = 50 . (39)The atom balance takes the form: Dc = (cid:20) (cid:21) c A c B c C = 1 , (40)and the equilibrium point is chosen as: c eqA = 0 . c eqB = 0 . c eqC = 0 .
4. As Fig.4 shows, the system is effectively two-dimensional, so the quasi-equilibriummanifold is expected to provide an one-dimensional reduced description ( q =1). Indeed, any solution trajectory, after a rapid initial dynamics, is attractedto a 1D curve and along it reaches the equilibrium point. If the system is closedand the reaction (39) takes place under constant volume and temperature, wecan assume that the system is supported by the Lyapunov function G (7): G = c A " ln c A c eqA ! − + c B " ln c B c eqB ! − + c C " ln c C c eqC ! − . (41)14nce a 3-dimensional vector m has been chosen, the QEM equation can befound by solving the variational problem (23): G → min( m , c ) = ξ Dc = 1 . (42)In the following, the Spectral Quasi Equilibrium Manifold (SQEM) [2,1] willbe constructed. The Jacobian matrix L in the equilibrium point and its slowestleft eigenvector x sl are: L ( c eq ) = − − . . − − . . − . , x sl = (cid:20) . , − . , . (cid:21) . (43)Solution of the problem (42), with the choice m = x sl , delivers the 1D SQEMfor the case shown in Fig 4. To this end, let us rewrite the (42) in a moreexplicit form: c A = 0 . . ξ − . φ ( ξ ) c B = 0 . − . ξ − . φ ( ξ ) c C = φ ( ξ ) ∂G ( φ,ξ ) ∂φ = 0 , ∂ G ( φ,ξ ) ∂φ > , (44)where c = [ c A , c B , c C ] is the solution of the problem (42), while φ denotesthe relation between c C and the reduced variable ξ on the SQEM. By usingthe G function (41), the problem (44) is equivalent to the implicit equation . . ξ − . φ . ! − . . − . ξ − . φ . ! − . φ . ! − . (45)The solution of (45), by means of relations (44), gives the SQEM shown inFig. 5(a). One may now apply the QEG-algorithm described above, in orderto make a comparison with the analytic solution just found. An orthonormalbasis { ρ i } in the null space of the matrix D = (cid:20) (cid:21) has dimension z = 2and can be chosen as follows: ρ = [ − . , . , − . , ρ = [ − . , − . , . . (46)15 c B c C
0 0.1 0.3 0.5 0.7 0.9 1 00.10.20.30.40.5 c B c C SQEMSQE−grid (left branch)SQE−grid (right branch)SQEM Equilibrium Equilibrium (a) (b)
Fig. 5. (a) The bold curve is the SQE-manifold which was analytically evaluatedby solving the (45). In that case the SQEM represents a very good approximationof the invariant manifold. (b) The SQE-manifold is compared with the SQE-gridwhere ε = 10 − . Since the matrix E has the form: E = . − . . , (47)a vector t spanning ker( E ) is: t = [ − . , − . , . . The system (34), in this example, simply reads: h t , ρ i µ + h t , ρ i µ = − ( t , ∇ G ) µ + µ = ε . (48)By solving (48) in a QEG-node c n , the shift vector ˆ δ c n = µ ρ + µ ρ allowsto evaluate the new QEG-node c n +1 = c n + ˆ δ c n . The QEG procedure, startingfrom the equilibrium point c eq = c , was performed twice, keeping uniformlyparameter ε = 10 − . The first time, by choosing the solution in such a waythat c B n +1 < c B n , the left branch of the SQE-grid was obtained; then, byimposing c B n +1 > c B n , also the right branch was calculated. The algorithmwas terminated as soon as at least one component of the new node c n +1 becomes negative. The result, shown in Fig. 5(b), proves that the SQE-grid isin excellent agreement with the analytical curve (SQEM).16 c B c C SQEM ε =10 −3 ε =5 ⋅ −3 ε =10 −2 ε =5 ⋅ −2 Fig. 6. SQEG Left branch of the case in Fig. 5 (b). Different approximations com-pared with the analytical solution (SQEM). Each grid is calculated by using adifferent parameter ε . There is no need to stress the importance of the grid spacing parameter ε forthe QEG accuracy. In the case of Section 6, the SQEG was computed severaltimes with different values of ε . The QEG algorithm is based on the linearapproximation (30). Therefore, the smaller is (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) = ε the more accurate isthe QEM description by means of the QEG. Nevertheless, the smaller is ε thelarger is the number of times that the system (33) must be solved to have agrid of a fixed size. For this reason, we need to keep ε as large as possible.We estimated (at least the order of magnitude) the upper limit of spacing( ε u ) which gives a QEG “not far” from the relative QEM. From our numericalexperiments, a reasonable value for that was ε u ∼ = 10 − . As Fig. 6 shows, theQEG is not far from the QEM even for a quite coarse grid ( ε > ε u ). The QEG algorithm, which has been developed for constructing 1-dimensionalgrids, can be modified in order to get multi-dimensional grids, wheneverneeded. From all reasonable extension strategies, two of them here will beanalyzed: a straightforward extension and a flag extension (flag extension, forinvariant grids, was introduced in Ref. [3]). In the first case, the algorithmof paragraph 5 and the equation system (34) are tuned for a q -dimensionalgrid calculation. Here, the implicit assumption is that the grid dimension q isfixed and uniform everywhere in the phase space (like for the QEM construc-tion). However, a second flexible approach, suitable for SQEG construction,was developed, too. In that case, the grid dimension can be varied at will.17 .1 The straightforward extension According to the straightforward extension, if a node c n close to the q -dimensionalQEM is known, then a new node c n +1 can be added to the QE-grid by shifting c n : c n +1 = c n + ˆ δ c n , ˆ δ c n = X zi =1 µ i ρ i , (49)where { ρ i } is still a basis in the null space of matrix D . The linear constraintsof the problem (23) define the tangent space T to the G level surfaces in thenew node c n +1 . Let c be a generic point of T , the line – l passing from c n +1 and c has the parametric form: c = ϕ ˜ t + c n +1 , where ˜ t is the vector of T whichspans – l and ϕ is the parameter. The generalized form of the relations (26) is: ( m , c ) = ϕ (cid:16) m , ˜ t (cid:17) + (cid:16) m , c in +1 (cid:17) ⇒ (cid:16) m , ˜ t (cid:17) = 0 , ∀ ˜ t ∈ T ...( m q , c ) = ϕ (cid:16) m q , ˜ t (cid:17) + (cid:16) m q , c in +1 (cid:17) ⇒ (cid:16) m q , ˜ t (cid:17) = 0 , ∀ ˜ t ∈ T ( d i , c ) = ϕ (cid:16) d i , ˜ t (cid:17) + (cid:16) d i , c in +1 (cid:17) ⇒ (cid:16) d i , ˜ t (cid:17) = 0 , ∀ ˜ t ∈ T, (50)which means that vector ˜ t belongs to the null space of the matrix E (ker E ): E = m ... m q D . (51)Now, the dimension of basis { t j } in ker( E ) is ( z − q ). Since the quasi equi-librium condition requires that, among all the points c of T , c n +1 has theminimal value of G , the following orthogonality conditions hold:( ∇ G ( c n +1 ) , t j ) = 0 , ∀ j = 1 , ..., z − q. (52)For small vector ˆ δ c n , the approximation (30) can be used, so that the (52)become: X zi =1 h t j , ρ i i µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − q. (53)As the system (53) shows, the larger is the QEM dimension ( q ) the smalleris the set of “mere” quasi-equilibrium equations available, while the numberof unknowns remains constant ( z ). The closure of the rectangular system (53)requires q more equations and has only to do with the geometric structurewhich we want to provide the grid with (e.g. grid spacing, shift vector ori-entation in the phase-space, etc). In general, the geometric structure of the18 ig. 7. 2D Quasi Equilibrium Manifold. Location of solutions of the system (54) inthe phase space. grid under construction can be chosen at will: therefore there is no uniquegeometric closure for that system. However, one possible condition could beimposed, like in (33), by fixing the Euclidean norm of shift vector: (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) = ε .Nevertheless, ( q −
1) geometric constraints are still missing. In order to illus-trate how the geometric closure issue can be overcome, the case q = 2 willbe considered in the following. For that special case, a possible closure, whichcan be easily generalized, will be presented. If a two-dimensional QEG has tobe constructed, then only one extra equation is needed to close the system: P zi =1 h t j , ρ i i µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) = ε. (54)Fig. 7 shows that all the possible solutions of (54) are located, as a “crown”,near the QEM. A way to choose only two of them can be achieved by intro-ducing a new fixed vector ˜ m and imposing a given angle ϑ between ˜ m andˆ δ c n : X zi =1 ( ˜ m , ρ i ) µ i = (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) · k ˜ m k cos ϑ. (55)The choice ϑ = π/ X zi =1 ( ˜ m , ρ i ) µ i = 0 . (56)(56) allows to write a closed system: P zi =1 h t j , ρ i i µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − P zi =1 ( ˜ m , ρ i ) µ i = 0 (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) = ε, (57)where the extra information, through ε and ˜ m , concerns the grid spacingand the phase-space zone of interest where the grid must be constructed. In19eneral, the geometric closure of (53) can be achieved when ( q −
1) independentvectors { ˜ m i } and the parameter ε are fixed. Here, we present an approachwhich allows to get a rectangular structured grid. The general form of (57) is: P zi =1 h t j , ρ i i µ i = − ( t j , ∇ G ( c n )) , ∀ j = 1 , ..., z − q P zi =1 ( ˜ m j , ρ i ) µ i = 0 , ∀ j = 1 , ..., q − (cid:13)(cid:13)(cid:13) ˆ δ c n (cid:13)(cid:13)(cid:13) = ε. (58)The q -dimensional grid construction is split in q subsequent steps. Startingfrom the equilibrium c eq , system (58) is solved by choosing ( q − m j vectorsamong the q available and imposing: ˜ m j = m j ∀ j = 1 , ..., q −
1. In thisway, a first set of QEG nodes is attained as soon as ε is known. Now, startingfrom each of those points, system (58), by using a different combination of m j vectors, gives some more nodes. The procedure ends ( q -th step) when all thepossible different combinations of ( q −
1) vectors { m j } are over. In Section 8,that idea will be explained by means of an illustrative example. A multi-dimensional QE-grid construction becomes non-trivial especially when q becomes large. As reported in section 7.1, the straightforward extension re-quires to introduce some additional vectors ˜ m j . The flag extension can beapplied when a SQEG is searched. That procedure is strongly based on thealgorithm presented in paragraph 5 and it naturally leads to a rectangularstructured grid. The idea, which is behind, is simple and makes this methodvery flexible and really suitable for constructing high-dimensional rectangu-lar grids. Let us suppose that q is the grid dimension and the q SQE-vectors { m , ..., m s } are fixed. Here, the assumption is that m is the slowest eigen-vector (corresponding to the smallest eigenvalue by absolute value), m thesecond slowest and so forth. The grid construction is achieved in s subse-quent steps. In each step one more dimension is added to the grid. At thebeginning, by using m = m , the algorithm in section 5 provides the 1Dquasi-equilibrium grid. Now, starting from any node c ∗ of that grid, a new 1DQEG is constructed where m = m . In this case, the second QEG representsa trajectory on the 2D-manifold attracted to the slowest 1D-manifold in thenode c ∗ , once the fast dynamics is exhausted (see Fig. 8). G function dependson the equilibrium point c eq : G = G ( c , c eq ). Since c ∗ can be considered asa “local equilibrium” for the fast motion, the second 1D grid is obtained byminimizing G = G ( c , c ∗ ). Once the previous step is completed, the grid canbe extended in the third dimension by adding, in each node c ′ of the new 2Dgrid, a 1D QEG where m = m and G = G ( c , c ′ ). In this way, the procedureis performed up to a q -dimensional grid. By extending partially a previous20 ig. 8. A 2D flag. Once the 1D quasi equilibrium grid is found, from each node c ∗ ,new 1D quasi equilibrium grids are added. The second slowest 1D grid representsthat trajectory collected by the first 1D quasi equilibrium grid in the node c ∗ . grid, it becomes really easy to have some grids whose dimension is differentin different phase space zones. It is worth to stress that the straightforwardand the flag extension deliver two different objects: the first one just givesthe quasi-equilibrium grid “brute force”, while the second one is its conve-nient “approximation” which has some useful features as it will be illustratedin the following. First of all, the flag grid does not demand any extra vec-tor for the geometric closure and the grid dimension can be easily varied indifferent phase-space zones. Secondly, if a grid refinement procedure (MIG)is used in order to get an invariant grid out of the quasi-equilibrium one [2],then the flag extension reveals to be a very useful tool. Indeed, let us assumethat a multi-dimensional invariant grid is required in order to reduce a givenmodel. A possible strategy might be given by a “hybrid procedure” where theQEG algorithm and the MIG method are alternatively used according to thesequence: •
1D quasi-equilibrium grid construction (slowest grid); • MIG refinements until the 1D invariant grid is obtained; • flag extension from 1D invariant grid to 2D quasi-equilibrium grid; • MIG refinements until the 2D invariant grid is obtained; • flag extension from 2D invariant grid to 3D quasi-equilibrium grid; • MIG refinements...
The latter suggestion sheds light on one more option which, if implementedduring the flag extension, allows to go beyond the SQEG approximation of theinvariant manifold. Let us assume that the hybrid procedure of Section 7.2 is21tilized and a k -dimensional invariant grid (let c ∗ be its generic node) has to beextended to a ( k + 1)-dimensional grid. That grid will approximate the ( k + 1)-dimensional invariant grid, better than the SQEG does, if in each invariantnode c ∗ the vector m is chosen as the ( k + 1)-th slowest left eigenvector (byabsolute value) of Jacobi matrix L ( c ∗ ). According to [1], here a considerablesimplification can be achieved by replacing the full Jacobian L ( c ∗ ) with: L sym ( c ∗ ) = 12 (cid:16) L ( c ∗ ) + H − L T ( c ∗ ) H (cid:17) , (59)where L T is the ordinary transposition, and H is evaluated at the point c ∗ ,too. Matrix L sym is symmetric with respect to the entropic scalar product(8). For that reason the spectral decomposition will be much more viable (seealso Ref. [2]). Those two new approximations will be named: Guided QuasiEquilibrium Grid (GQEG) when the full Jacobian L ( c ∗ ) is used, while Sym-metric Entropic Guided Quasi Equilibrium Grid (SEGQEG) if L sym replacesthe full matrix. In order to give an idea about the effort needed, for examplein a SEGQEG construction, let us consider a 2-dimensional grid. In that case,the spectral decomposition of a symmetric operator is performed only overthe nodes of an one-dimensional grid. Moreover, also a possible criterion, forgetting a multi-dimensional grid, naturally applies: if at the node c ∗ of the k -dimensional invariant grid, the ratio | λ k +1 | / | λ k | (between eigenvalues of L or L sym , respectively) is not larger than a fixed threshold, the ( k +1)-dimensionalgrid will not be extended at that point. In this way, the grid dimension q isgenerally not uniform in the phase space. Several techniques suggested aboveare only some reasonable ones. The flexibility of the method proposed allowsto set up different procedures, still based on the Quasi Equilibrium Grid ap-proach: the QEG system (53) supplied by a geometrical closure. Let us consider a model for hydrogen oxidation reaction where six species H (hydrogen), O (oxygen), H O (water), H , O , OH (radicals) are involved insix steps in a closed system under constant volume and temperature (see Ref.224], p. 291): .H ↔ H, k +1 = 2 , .O ↔ O, k +2 = 1 , .H O ↔ H + OH, k +3 = 1 , .H + O ↔ H + OH, k +4 = 10 , .O + H ↔ O + OH, k +5 = 10 , .H + O ↔ H O, k +6 = 10 . (60)The conservation laws are: c H + 2 c H O + c H + c OH = b H = 22 c O + c H O + c O + c OH = b O = 1 . (61)When the equilibrium point is fixed, for example c eqH = 0 . , c eqO = 0 . , c eqH O = 0 . , c eqH = 0 . , c eqO = 0 . , c eqOH = 0 . , (62)then the rest of the rate constants k − i are calculated using the detailed balanceprinciple (4). The system under consideration is fictitious in the sense thatthe subset of equations corresponds to the simplified picture of this chemicalprocess and the rate constants reflect only orders of magnitude for relevantreal-word systems. We can assume that the Lyapunov function G has the form: G = X i =1 c i " ln c i c eqi ! − . (63)Here, we are interested in the 2D SQEG construction. Two left eigenvectorsof Jacobian matrix L ( c eq ) are: x s l = [ − . , − . , . , . , . , − . x s l = [0 . , − . , . , − . , − . , . , (64)where x s l and x s l are the slowest and the second slowest one, respectively.23 c H c O c O H Fig. 9. The 2D SQEG constructed by using the straightforward extension with ε = 0 . · − : projection into the phase-subspace ( c H , c O , c OH ). In order to get a 2D SQEG for that example, the straightforward extensionwas used as first strategy. Matrices D and E take now the form: D = , E = − . − .
568 0 .
225 0 . − . − . . − . . − . − .
713 0 . . (65)As suggested in the end of section 7.1, the procedure has been started from theequilibrium point and it was split in two subsequent steps. At the beginning,the system (57) was solved by imposing ε = 0 . · − and ˜ m = x s l : in thisway, the grid nodes, denoted by circles, in Fig. 9 were obtained. In the secondstep, (57) was solved by starting from any circle: this time, the geometricconstraints were ε = 0 . · − and ˜ m = x s l . During that step, in each circle,the horizontal dots of Fig. 9 were found, too.24 .03 0.035 0.04 0.045 0.050.0120.0140.0160.0180.020.0220.0240.0260.02800.0050.010.0150.020.0250.03 c H c O c O H Equilibrium1D SQEG (down branch)1D invariant grid (down branch) 1D SQEG (upper branch)1D invariant grid (upper branch)
Fig. 10. 1D Spectral Quasi Equilibrium Grid with ε = 3 · − : comparison withthe 1D invariant grid obtained by MIG refinements. After that, also the flag extension procedure was applied. Now, the 1D spectralquasi equilibrium grid is needed. Matrices D and E are in this case: D = , E = − . − .
568 0 .
225 0 . . − . . (66)Starting from the equilibrium point c eq , the system (34) was solved by fixing ε = 3 · − (see Fig. 10). The flag extension was used to get a 2D grid outof the 1D one. Now, the new matrix E reads: E = . − . . − . − .
713 0 . , while the Lyapunov function G has the form: G = X i =1 c i " ln c i c ∗ i ! − , (67)where c ∗ = [ c ∗ , ..., c ∗ ] is any 1D grid node which is extended in the seconddimension (see Fig. 8). Figures 11(a)-(b) show two different 2D SQE-grids:the first one is obtained by extending the 1D SQE-grid, while in the secondcase the 1D invariant grid is used. In other words, the latter result was attainedby the “hybrid procedure” QEGA + M IG suggested in the end of section 7.2.For both cases, in the second dimension, the grid spacing was ε = 1 . · − .25 c H c O c O H c H c O c O H (a)(b) Fig. 11. The flag extension. (a) 2D SQEG (dots) extended from the 1D SQEG(circles). (b) 2D SQEG (dots) extended from the 1D invariant grid (circles). Thegrid spacing, in the second dimension, was ε = 1 . · − . Finally, the GQEG and SEGQEG approximations are computed for the hy-drogen oxidation reaction (60) (1.case). Here, the grid spacing is uniformlykept ε = 0 . · − . Each grid has 10 ×
15 nodes and it is compared withboth the SQEG (straightforward extension) of similar size (check Table 1) andthe invariant one. The invariant grid was obtained by refining the approxima-tions through the MIG procedure. All those grids lie quite close to each other.However, a “more pathological” case (2.case) is also analyzed (here the SQEG,far from the equilibrium, presents a remarkable deviation from the invariantgrid): now the rate constant set is taken as k +1 = 20, k +2 = 1, k +3 = 1, k +4 = 10 , k +5 = 10 , k +6 = 10 , while the equilibrium point coordinates still are givenby (62). For that case, the SQEG, GQEG and SEGQEG were constructed bychoosing the grid spacing and size as for the previous case. Note that all thegrids were partially extended only below the equilibrium point in the phase-space zone where they present the largest deviation from the invariant grid.This time those three approximations have a low invariance defect only nearthe equilibrium. In order to estimate how far each grid is from the invariantone, the following procedure is implemented. A 10 ×
15 matrix, collecting inany grid node an invariance defect measure, is constructed. As suggested by26 .050.10.150.20.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.112345678910 x 10 −3 c H c O c O H Fig. 12. 1.case: k +1 = 2, k +2 = 1, k +3 = 1, k +4 = 10 , k +5 = 10 , k +6 = 10 , ε = 0 . · . Two grids formed by 10 ×
15 nodes. A 2D GQEG (dots) and a2D invariant grid (circles) are reported. The grids are partially extended below theequilibrium point (square). c H c O c O H Fig. 13. 2.case: k +1 = 20, k +2 = 1, k +3 = 1, k +4 = 10 , k +5 = 10 , k +6 = 10 , ε = 0 . · . Two grids formed by 10 ×
15 nodes. A 2D GQEG (dots) and a2D invariant grid (circles) are reported. The grids are partially extended below theequilibrium point (square). [2], that local measure may be q ( ∆ , ∆ ) / ( J , J ), where ∆ and J are the in-variance defect (11) and the vector field of (5), respectively. ∆ is evaluated byusing the thermodynamic projector (16). By averaging over all the invariancedefect measures, the mean invariance defect is provided: results for both casesare condensed in Table 1. Note that the adopted invariance defect measure isdimensionless as it compares the invariance defect with the vector field. Cal-27 .case 2.caseSQEG 0.318 0.645GQEG 0.238 0.460SEGQEG 0.303 0.491Table 1Mean invariance defect (dimensionless): three approximations of the invariant gridunder comparison for the hydrogen oxidation reaction. In 1.case, the parameter setis: k +1 = 2, k +2 = 1, k +3 = 1, k +4 = 10 , k +5 = 10 , k +6 = 10 , ε = 0 . · . In 2.case,the parameter set is: k +1 = 20, k +2 = 1, k +3 = 1, k +4 = 10 , k +5 = 10 , k +6 = 10 , ε = 0 . · . culations prove that the GQEG is better than the SQEG (straightforwardlyextended); nevertheless the SEGQEG construction, since it requires a muchlower computational effort and still has an error similar to the GQEG, isrecommended when the SQEG is considered not satisfactory (e.g. big meandefect). In this paper, the problem of Quasi Equilibrium Manifold approximation bymeans of a grid description is addressed. To this end, the notion of Quasi Equi-librium Grid (QEG) is introduced and a proper algorithm to construct it, inany dimension, is suggested (QEGA). It has been shown, through illustrativeexamples, that the QEGA gives a very good QEM approximation without fac-ing the analytical difficulties of Lagrange multipliers method implementationin large dimension. Since the QEGA is a completely numerical procedure, itreveals particularly suitable for providing the MIG procedure with the firstSIM approximation. As it has been illustrated, some proper hybrid procedures
QEGA + M IG , where both methods are alternatively used, allow to obtainaccurate SIM approximations. It was proved that two special
QEGA + M IG procedures deliver enhanced approximations of SIM: the
Guided Quasi Equi-librium Grid and the
Symmetric Entropic Guided Quasi Equilibrium Grid .Here, we want to stress the two major advantages of the method proposed.First of all, it is a completely numerical algorithm which only deals with nodessets. Moreover, it is a local construction: namely, the computation of a newnode c n +1 , which has to be added to the grid, only depends on the previousneighbor c n . Those two points make the QEG construction suitable for numer-ical applications and parallel realizations. Finally, it is not excluded that theQEGA is applicable not only for model reduction, but in some very differentfields, too. Indeed, it was mentioned that the QEM notion already is exploitedfor some applications in Lattice Boltzmann schemes simulations. More gen-erally, the QEGA is a numerical tool which can be used to find a grid-based28pproximation for the locus of minima of a convex function under some lin-ear constraints. In this paper we focused only on the geometry of the modelreduction, that is, construction of slow invariant manifolds approximations.The implementation of grid-based integrators for dynamic equations will bepresented in a separate publication.
10 Acknowledgments
Prof. A. N. Gorban is gratefully acknowledged for the fruitful discussions aboutthe concept of Quasi Equilibrium Manifold and for suggesting the reactionmechanism (39). We thank Prof. K. B. Bouluochos, Dr. C. E. Frouzakis forseveral discussions and suggestions. This work was partially supported bySNF, Project 200021-107885/1 (E.C.) and by BFE, Project 100862 (I.V.K.).
References [1] A. N. Gorban, I. V. Karlin, Method of invariant manifold for chemicalkinetics, Chem. Eng. Sci. 58 4751-4768 (2003).[2] E. Chiavazzo, A. N. Gorban, I. V. Karlin, Comparison of invariant manifoldsfor model reduction in chemical kinetics, Comm. in Comput. Physics Vol. 2,No. 5, pp. 964-992 (2007).[3] A. N. Gorban, I. V. Karlin, A. Y. Zinovyev, Invariant grids for reactionkinetics, Physica A 333 106-154 (2004).[4] A. N. Gorban, I. V. Karlin, Invariant Manifolds for Physical and ChemicalKinetics, Lect. Notes Phys. 660 (Springer Berlin Heidelberg 2005), doi:10.1007/b98103.[5] A. N. Gorban, I. V. Karlin, V. B. Zmievskii, T.F.Nonnenmacher,Relaxational trajectories: global approximations, Physica A 231 648-672(1996).[6] I. V. Karlin, A. Ferrante, H. C. ¨Ottinger, Perfect entropy functions of thelattice Boltzmann method, Europhys. Lett. 47, 182-188 (1999).[7] S. Arcidiacono, J. Mantzaras, S. Ansumali, I. V. Karlin, C. Frouzakis, K. B.Boulouchos, Phys. Rev. E 74, 056707 (2006).[8] A. N. Gorban, I.V.Karlin, Thermodynamic parametrization. Physica A 190393-404 (1992).[9] U. Maas, S. B. Pope, Simplifying Chemical Kinetics: Intrinsic Low-Dimensional Manifolds in Composition Space, Comb. and Flame 88 239-264(1992).
10] V. Bykov, I. Goldfarb, V. Gol’dshtein, U. Maas, On a modified version ofILDM approach: asymptotic analysis based on integral manifolds, IMA Jour.of Applied Math. 1-24 (2005).[11] S. H. Lam, D. A. Goussis, Conventional asymptotic and computationalsingular perturbation for symplified kinetics modelling , in: M.O. Smooke(Ed.), Reduced Kinetic Mechanisms and Asymptotic Approximations forMethane-Air Flames, Springer Lecture Notes, (Springer Berlin 1991), pp.227-242.[12] S. H. Lam, D. A. Goussis, The CSP Method for Simplifying Kinetics, Intern.Jour. of Chemical Kinetics 26 461-486 (1994).[13] D. A. Goussis, M. Valorani, An efficient iterative algorithm for theapproximation of the fast and slow dynamics of stiff systems, Jour. ofComputational Physics 214 316-346 (2006).[14] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Runborg,C. Theodoropoulos, Equation-free: Coarse-grained Multiscale ComputationEnabling Microscopic Simulators to Perform System-level Analysis, Comm.Math. Sci. 1 715-762 (2003).[15] A. M. Lyapunov, The general problem of the stability of motion, Taylor &Francis, London, (1992).[16] N. Kazantzis, Singular PDEs and the problem of finding invariant manifoldsfor nonlinear dynamical systems. Phys. Lett. A 272, 257-263 (2000).[17] J. C. Robinson, A concise proof of the “geometric” construction of inertialmanifolds, Phy. Lett. A 200, 415-417 (1995).[18] L. B. Ryashko, E. E. Shnol, On exponentially attracting invariant manifoldsof ODEs, Nonlinearity 16, 147-160 (2003).[19] R. S. MacKay, Slow Manifolds, in: Energy Localisation and Transfer, ed. byT. Dauxois, A. Litvak-Hinenzon, R. S. MacKay, A. Spanoudaki (World Sci.,2004), 149-192.[20] A. N. Gorban, I. V. Karlin, Uniqueness of thermodynamic projector andkinetic basis of molecular individualism, Physica A 336, 3-4, 391–432 (2004).[21] P. D. Christofides, P. Daoutidis, Finite-dimensional control of parabolic PDEsystems using approximate inertial manifolds, J. Math. An. & Appl. 216,398-420 (1997).[22] A. Degenhard, J. Rodriguez-Laguna, Projection Operators for NonlinearEvolutionary Dynamics, SIAM Multiscale Modeling and Simulation 4, 641-663 (2005).[23] H. G. Kaper, T. J. Kaper, Asymptotic analysis of two reduction methods forsystems of chemical reactions, Physica D 165 66-93 (2002).
24] N. Kazantzis, Singular PDEs and the problem of finding invariant manifoldsfor nonlinear dynamical systems, Physics Letters, A272(4), 257-263 (2000).[25] M. R. Roussel, S. J. Fraser, On the geometry of transient relaxation, J. Chem.Phys. 94, 7106-711 (1991).[26] M. R. Roussel, S. J. Fraser, Geometry of the steady-state approximation:Perturbation and accelerated convergence methods, J. Chem. Phys. 93, 1072-1081 (1990).[27] S. J. Fraser, The steady state and equilibrium approximations: A geometricalpicture, J. Chem. Phys. 88 4732-4738 (1988).24] N. Kazantzis, Singular PDEs and the problem of finding invariant manifoldsfor nonlinear dynamical systems, Physics Letters, A272(4), 257-263 (2000).[25] M. R. Roussel, S. J. Fraser, On the geometry of transient relaxation, J. Chem.Phys. 94, 7106-711 (1991).[26] M. R. Roussel, S. J. Fraser, Geometry of the steady-state approximation:Perturbation and accelerated convergence methods, J. Chem. Phys. 93, 1072-1081 (1990).[27] S. J. Fraser, The steady state and equilibrium approximations: A geometricalpicture, J. Chem. Phys. 88 4732-4738 (1988).