Reduction of dimension for nonlinear dynamical systems
aa r X i v : . [ n li n . C D ] A ug Noname manuscript No. (will be inserted by the editor)
Reduction of dimension for nonlinear dynamical systems
Heather A. Harrington and Robert A. Van Gorder
Received: date / Accepted: date
Abstract
We consider reduction of dimension for non-linear dynamical systems. We demonstrate that in somecases, one can reduce a nonlinear system of equationsinto a single equation for one of the state variables, andthis can be useful for computing the solution when us-ing a variety of analytical approaches. In the case wherethis reduction is possible, we employ differential elimi-nation to obtain the reduced system. While analytical,the approach is algorithmic, and is implemented in sym-bolic software such as
MAPLE or SageMath . In othercases, the reduction cannot be performed strictly interms of differential operators, and one obtains integro-differential operators, which may still be useful. In ei-ther case, one can use the reduced equation to bothapproximate solutions for the state variables and per-form chaos diagnostics more efficiently than could bedone for the original higher-dimensional system, as wellas to construct Lyapunov functions which help in thelarge-time study of the state variables. A number ofchaotic and hyperchaotic dynamical systems are usedas examples in order to motivate the approach.
Keywords reduction of dimension · differentialelimination · nonlinear dynamics · chaotic attractors · computation of chaos H.A. HarringtonMathematical Institute, University of Oxford, Andrew WilesBuilding, Radcliffe Observatory Quarter, Woodstock Road,Oxford OX2 6GG United KingdomE-mail: [email protected]. Van GorderMathematical Institute, University of Oxford, Andrew WilesBuilding, Radcliffe Observatory Quarter, Woodstock Road,Oxford OX2 6GG United KingdomE-mail: [email protected]
Nonlinear dynamical systems are ubiquitous in mathe-matics, engineering, and the sciences, with many real-world phenomenon governed by such nonlinear processes.In particular, non-equilibrium and chaotic dynamics area continuing area of active research for applied mathe-maticians, as approximating such dynamics accuratelyand efficiently can be quite challenging. In the presentpaper, we shall consider reduction of dimension for non-linear dynamical systems. This approach has previouslybeen employed in the literature in order to enable theconstruction of Lyapunov functions [1] and equilibriumdynamics [2], as well as to allow one to more easilyapproximate chaotic attractors analytically [3,4,5,6].One method for reduction of dimension is differentialelimination, in which one algorithmically reduces thenonlinear dynamical system into a single ordinary dif-ferential equation (ODE) for one of the state variables.However, this is possible only when the system reducesto an ODE; if the reduction is instead to an integro-differential equation, the process is not algorithmic andspecific cases must be handled with more individualcare. Our focus shall be on dynamical systems giv-ing chaotic dynamics, but the approach can certainlybe applied for non-chaotic ODE systems. We give anoverview of reduction of dimension, after which we demon-strate in several ways why one might wish to apply thistechnique.With the wide range of numerical methods availablefor solving nonlinear first-order ODE systems of evenhigh order, one may wonder why it might be advanta-geous to convert such systems into a single higher-orderODE. We shall mention several situations in which thedifferential elimination, and more generally reduction
Heather A. Harrington and Robert A. Van Gorder of dimension, may prove useful. We then outline thepaper.Often times, if one is trying to approximate the so-lution to a nonlinear system through some sort of ana-lytical approximation, via series, perturbation, or morecomplicated approaches, one quickly finds that the cou-pled equations require balancing many terms comingfrom the expansion for each of the state variables. In thecase of a single state variable governed by a higher-orderODE, one need only track terms in a single asymptoticexpansion. This approach has been applied when usingTaylor series, approximate Fourier series, and asymp-totic expansions in other types of basis functions tothe solution of a system of nonlinear ODE. In hybridanalytic-numeric methods, such as the homotopy anal-ysis method [7,8,9], such reductions of a system to asingle equation also simplifies the optimization problemwhich is solved to obtain the error-minimizing solution(see, for instance, [2], where the present approach isused in such a capacity). Therefore, the reduction oforder can greatly reduce the complexity of analyticalcalculations under several frameworks.Contraction maps or Lyapunov functions are usefultools for discussing the convergence of solutions to non-linear dynamical systems to large-time steady or quasi-steady dynamics. In situations where contraction mapsor Lyapunov functions are known for a given dynamicalsystem, the state variable governed by a single higher-order ODE necessarily results in a contraction map inthe single state variable. However, as is well known tothose studying stability of nonlinear systems, it is notoften easy to obtain contraction maps for complicatedsystems. As we shall show here, it is possible to usethe reduction of a system to a single higher-order ODEin order to construct a contraction map for the statevariable governed by the aforementioned higher-orderODE. The existence of such a map can then be used todeduce the large-time dynamics of the state variable,as well as for the other state variables in the originalsystem. One example of this is given in [1], and otherexamples are provided in Section 5.Related to both the topic of analytical approxima-tions and Lyapunov functions would be long time dy-namics and equilibrium behavior of nonlinear dynami-cal systems. Indeed, in order to study the equilibriumstructure of a high-order system of ODEs, one mustsolve a coupled system of nonlinear algebraic equationsin order to recover the fixed points for the state vari-ables. First reducing the system to a single ODE allowsone to obtain a single nonlinear algebraic equation forthe fixed point of a single state variable, which can thenbe used to recover the fixed points of the other statevariables. Therefore, when such a reduction to a single ODE is possible, the need to solve a nonlinear alge-braic system for all of the fixed points simultaneouslyis eliminated, resulting in what is often a far less com-putationally demanding problem.Another topic is great recent interest in nonlinearscience has been both the synchronization of chaos [10,11,12] and the control of chaos. In situations where oneis interested in mitigating the possibility of emergentchaos, one can couple a chaotic system to various con-trol terms, or indeed to additional dynamical systems,which may lend a degree to stability. Under such ap-proaches, one often increases the complexity or eventhe dimension of the dynamical system being solved.As such, methods to reduce the dimension of such sys-tems could improve compatibility. Furthermore, sincethe control of chaos is often linked to a control termwhich itself is determined by a Lyapunov function, theconstruction of contraction maps through the reductionapproach outlined here could be of great use.As stated before [13], the competitive modes anal-ysis gives an interesting link between the geometry ofphase space possibly yielding chaotic trajectories (recallthat the competitive modes requirements appear to bea necessary, albeit not sufficient, condition for chaos[14,15,16,17,18,19]). Conversely, the differential elimi-nation may cast light into the geometry of solutions inthe space of derivations. Since this result of the differen-tial elimination is a single higher-order ODE, and sinceany chaos emergent from the nonlinear system shouldbe encoded in the single higher-order ODE, the differ-ential algebraic structure of such equations may castlight into practical geometric tools by which one maystudy systems in which chaos is observed. In particu-lar, through this reduction approach, the calculation ofmode frequencies in the standard competitive modesanalysis becomes much simpler.This paper is outlined as follows. In Section 2, weprovide an algorithmic approach, to differential elim-ination for nonlinear dynamical systems based upondifferential algebra. First laying out the general the-ory, we then give specific
MAPLE code for performingthe differential elimination in a systematic manner. Thealgorithmic approach is useful in the case where thedynamical system can be reduced to a single ODE interms of only one of the state variables. In Section 3 weimplement the approach in order to reduce a varietyof chaotic and hyperchaotic systems, finding that theform of the nonlinearity in the dynamical system willstrongly influence the reducibility properties. However,in cases where the dynamical system is not reducible us-ing differential elimination, one may still obtain morecomplicated reductions, for instance in terms of inte-grals, resulting in more complicated integro-differential eduction of dimension for nonlinear dynamical systems 3 equations for the reduced state variable. The possibleresults are illustrated through concrete examples forthe R¨ossler system (which is completely reducible), theLorenz system (which is partially reducible - that is,reducible in some but not all state variables), and theQi-Chen-Du-Chen-Yuan (which is irreducible under dif-ferential elimination, but which can be reduced to anintegro-differential equation). We give summarizing ob-servations regarding the reducibility of dynamical sys-tems in Section 4.The remainder of the paper is devoted to applica-tions of reduction of dimension for dynamical systems.In Section 5, we demonstrate that reduction of dimen-sion can be useful for obtaining contraction maps andLyapunov functions, which in turn may be used to de-termine asymptotic stability of dynamical systems andalso to control chaos in such systems. In Section 6 wedemonstrate that reduction of dimension can be usedto simplify calculations involved in certain techniquesfor studying the solutions of nonlinear dynamical sys-tems. Indeed, when applicable, we find that the ap-proach greatly reduces the number of nonlinear alge-braic equations required to be solved when constructingtrajectories in state space via undetermined coefficientmethods by a factor of 1 /n , where n is the dimensionof the dynamical system, meaning that the number ofequations needing to be solved will not increase withthe size of the system. Furthermore, when applying thecompetitive modes analysis (which is a type of diagnos-tic criteria for finding chaotic trajectories in nonlineardynamical systems), we find that only one binary com-parison is needed of one first reduces the dimension ofthe dynamical system so that there is a single equationfor one state variable. In contrast, there are normallyof order 2 n − comparisons needed for an n -dimensionaldynamical system. In Section 7 we provide a discussionand possible avenues for future work. Systems of differential equations are ubiquitous andwidely studied. Ritt [20] and Kolchin [21] pioneeredthe field of differential algebra, an algebraic theory forstudying solutions of ordinary and partial differentialequations. We are particularly interested in differentialelimination, an algorithmic subtheory that can simplifysystems of parameterized algebraic differential equa-tions. This permits one to reduce the dimension of adynamical system so that one is left with a single ODEin the state variable. 2.1 Algebra preliminariesHere we briefly review concepts from algebra and dif-ferential algebra. For reference books in differential al-gebra, see [20,21]. If I is a subset of a ring R then ( I ) isthe (algebraic) ideal generated by I . Let I be an ideal of R . Then √ I denotes the radical of I . A derivation overa ring R is a map R R which satisfies (we write ˙ a isthe derivative of a ), for every a, b ∈ R , ˙( a + b ) = ˙ a + ˙ b and ˙( a b ) = ( ˙ a ) b + a (˙ b ) . The field of differential alge-bra is based on the concept of a differential ring (resp.field), which is a ring (resp. field) R endowed with a setof derivations that commutes pairwise. A differentialideal [ I ] of a differential ring R is an ideal of R stableunder the action of derivation.Differential algebra is more similar to commutativealgebra than analysis. In commutative algebra, Buch-berger solved the membership problem (tests whether agiven polynomial is contained in a given ideal) throughthe theory of Gr¨obner bases [30]. From algebraic geom-etry, we know the set of polynomials which vanish overthe solutions of a given polynomial system form an idealand even a radical ideal [31]. In the case of differentialequations, the set of differential polynomials which van-ish over the analytic solutions of a given system of dif-ferential polynomial equations form a differential idealand even a radical differential ideal [20]. Ritt solved thetheoretical problem (of membership for radical differ-ential ideals) and developed algorithmic tools to solvesystems of polynomial ODE and PDEs; however, Ritt’salgorithm requires factorization.Due to the complexity of factorization, Boulier andco-authors avoided it by developing the Rosenfeld-Gr¨obneralgorithm, based on the work of Seidenberg and Rosen-feld, and incorporating Gr¨obner bases [27,25,24]. Sincethen, the algorithm has been improved both theoreti-cally and practically [23,22,26] and it no longer requiresGr¨obner bases. It is available in the DifferentialAlgebra package in
MAPLE [22] and
SageMath as an interfacefor the BLAD and BMI libraries [32,33].Algorithmically, differential elimination involves ma-nipulation of finite subsets of a differential polynomialring R = K { U } where K is the differential field of coef-ficients (i.e. K = Q ), and U is a finite set of dependentvariables. The elements of R are differential polynomi-als, which are polynomials built over the infinite setof all derivatives ΘU , of the dependent variables. Con-sider a system Σ of polynomial differential equations,here, we consider the Lorenz system of three ordinary Heather A. Harrington and Robert A. Van Gorder differential equations:˙ x = a ( x − x ) , ˙ x = x ( b − x ) − x , ˙ x = x x − cx . (1)The Lorenz system can be re-written as: Σ = − ˙ x + a ( x − x ) = 0 , − ˙ x + x ( b − x ) − x = 0 , − ˙ x + x x − cx = 0 . (2)The Rosenfeld-Gr¨obner algorithm takes as an input adifferential system Σ and a ranking . A ranking > is anytotal ordering over the set ΘU of all derivatives of the el-ements of U which satisfies the following axioms: a < ˙ a and a < b ⇒ ˙ a < ˙ b for all a, b ∈ ΘU . The Rosenfeld-Gr¨obner algorithm transforms Σ into finitely many sys-tems called regular differential systems, which reducesdifferential problems to purely algebraic ones that aretriangular. The next step is purely algebraic and trans-forms the regular differential system into finitely manycharacteristic presentations, C , . . . C r . Rosenfeld-Gr¨obneroutputs this finite family C , . . . C r of finite subsets of K { U }\ K , where each C i defines a differential ideal [ C i ].The radical p [ Σ ] of the differential ideal generated by Σ is the intersection presented by characteristic sets: p [ Σ ] = [ C i ] ∩ · · · ∩ [ C r ] . Note differential ideals [ C i ] do not need to be prime,however by Lazard’s lemma, they are necessarily rad-ical. Differential algebra elimination has proven usefulfor parameter estimation, identifiability, and model re-duction of biological and chemical systems [29,28].2.2 Computational MethodWe demonstrate reduction in dimension via differen-tial elimination algorithm RosenfeldGroebner in the
DifferentialAlgebra package implemented in
MAPLE .For sake of using a concrete example, we choose theLorenz system. First, we call the package: with(DifferentialAlgebra):
Next we input the Lorenz system: sys := [-(diff(x1(t),t))+a*(x2(t)-x1(t)),-(diff(x2(t),t))+x1(t)*(b-x3(t))-x2(t),-(diff(x3(t),t))+x1(t)*x2(t)-c*x3(t)]
Next we form our differential ring, embedding the rankof dependent variables in blocks and independent vari-ables in derivations . Since we are considering ordi-nary differential equations, derivations is set to one or-dering, time t . We remark that the DifferentialAlgebra package enables differential elimination of PDEs by in-cluding additional inputs for the derivations (e.g., derivations=[u,x,t] . Note, sys is assumed to havecoefficients in the field Q [ x , x , x ] obtained by adjoin-ing the independent variables to the field of rationalsand symbolic parameters a, b, c are considered arbitraryin the coefficient field. To form the differential ring, weinput: R := DifferentialRing(blocks = [x3, x2, x1],derivations = [t])
Note that x stands to the rightmost place on the listwhich identifies that we are attempting to reduce thedifferential equation to only one variable, i.e. x ( t ). Thisranking eliminates x with respect to x , and then elim-inates x with respect to x . We now call the RosenfeldGr¨obner algorithm for our system and differential poly-nomial ring: G := RosenfeldGroebner(sys, R)simplify(Equations(G[1], solved))
This will return the characteristic presentation (whichshould be understood as an intersection), with the equa-tions given by the ranking, with the final equation asingle ODE for x1(t) , provided that it exists and canbe computed by the algorithm. In some cases the algo-rithm will keep running and therefore should eventuallybe terminated by the user. For such cases, it is unlikelythat a reduction of the specified form exists. However,as we shall consider in the next section, when the reduc-tion is to an integro-differential equation, rather thanan ODE, the approach will not identify the reducedequation.
Here we apply the method of differential eliminationto several nonlinear dynamical systems known to givechaos, in order to see if these equations can be reduced.We first apply the algorithmic approach outlined inSection 2, finding that the approach gives a completereduction (all state variables can be isolated and ex-pressed as the solution to single uncoupled ODEs), apartial reduction (one or more, but not all, state vari-ables can be isolated and expressed as the solution tosingle uncoupled ODEs), or returns no reduction (thealgorithm does not complete in a fixed amount of time),in which case none of the state variables can be ex-pressed as a solution to a single ODE reducible fromthe original system. For simplicity, we shall only con-sider autonomous systems.We consider a number of examples of chaotic sys-tems in Table 1, with the results of the differential elim-ination algorithm given. We also give a summary of thedynamics of the example equations selected. Since the eduction of dimension for nonlinear dynamical systems 5 form of these equations may vary through the literature,we give a list of the specific form of the equations con-sidered, in Appendix A. Note that we have consideredthe differential elimination algorithm for the arbitraryparameter values listed in Appendix A. Table 1 demon-strates that the structure of the dynamical system tendsto play a strong role in whether the system can be re-duced. Indeed, equations with a single nonlinearity tendto be completely or partially reducible, hence at leastone state variable can be solved for via a single nonlin-ear ODE. On the other hand, the equations with manynonlinear terms or higher-order degree of nonlinearity(we consider only equations with polynomial nonlinear-ities) tend more often to be irreducible using the ap-proach. Of the listed equations, note that the R¨osslersystem is one of the few completely reducible systems,lending validity to the belief that it is indeed one of thesimplest possible continuous- time dynamical systemsgiving chaos. Meanwhile, the commonly studied Lorenzsystem is only partially reducible under the approach.More complicated systems tend to be irreducible underthe algorithm, and many of these give more complicateddynamics such as multiple scroll attractors.Note that the algorithm returns a ‘No’ if a reductionis not obtained within a given time interval. For caseswhere the algorithm found a reduction, the computa-tion time was fairly quick. We are therefore comfortablein assuming that a reduction to an ODE does not existin cases where the the algorithm times out. For suchcases, the system may still admit a reduction, but notstrictly in derivatives of one of the state variables. Onesuch example would be a system which is reducible toan integro-differential equation in one of the state vari-ables, but never to simply an ODE.We next consider hyperchaotic systems (chaotic sys-tems giving two or more positive Lyapunov exponents)in Table 2. Again, we find that the more complicatedthe functional form of the nonlinarities, the less likelya system seems to be reducible. Furthermore, hyper-chaotic generalizations of known chaotic systems ap-pear to maintain their reducibility properties, since of-ten a simple additional equation is added to make achaotic system hyperchaotic. The hyperchaotic R¨osslersystem is completely reducible, as was the related chaoticR¨ossler system, again suggesting that the chaotic andhyperchaotic R¨ossler systems are some of the simplestsystems which still exhibit chaos and hyperchaos, re-spectively.The results indicate that completely reducible sys-tems are perhaps the simplest systems giving chaosor hyperchaos. Again, this would support the qualita-tive and topological claims that the R¨ossler systemsare some of the simplest possible equations permitting
System Dynamics Reducible?
Lorenz[34,35] 3D:1-2-2 PartialModified Chua’s circuit[36,37,38] 3D:3-1-1 CompleteChen-Lee[39,40] 3D:2-2-2 PartialRabinovich-Fabrikant[41,42,43] 3D:3-3-3 NoR¨ossler[44,45] 3D:1-1-2 CompleteChen[46,47] 3D:1-2-2 PartialL¨u[48,49] 3D:1-2-2 PartialT-system[50,51,52] 3D:1-2-2 PartialQi-Du-Chen-Chen-Yuan[53] 4D:3-3-3-3 NoQi-Chen-Du-Chen-Yuan[54] 3D:2-2-2 NoGeneralized Lorenz[55,56,57] 3D:2-2-2 NoBlue-sky catastrophe[58,59,60] 3D:3-3-3 NoLorenz-Stenflo[61,62,63] 4D:1-2-2-1 PartialGenesio-Tesi[64,65] 3D:1-1-2 PartialArneodo-Coullet-Tresser[66] 3D:1-1-3 Partial
Table 1
List of chaotic systems and their reduction prop-erties. The numbers in the dynamics column indicate thedimension of the system and then degree of each poly-nomial in the respective reaction functions. For instance,if ˙ x = A ( x, y, z ), ˙ y = B ( x, y, z ), ˙ z = C ( x, y, z ), then3D:deg( A ) − deg( B ) − deg( C ) is reported, where deg( A ) de-notes the degree of A , and so on. When the system has theproperty that it may be reduced to a single ODE in any statevariable, we say that it is completely reducible, and recorda ‘Complete’. If the system may be reduced to an ODE inone or more, but not all, state variables, we say the systemsis partially reducible, and record a ‘Partial’. Finally, when asystem is not reducible to a single ODE in any state vari-able, we record a ‘No’. We note that specific forms of someequations can change from paper to paper, so we record thespecific equations used in Appendix A. System Dynamics Reducible?
R¨ossler [67] 4D:1-1-2-1 CompleteChen [68,69] 4D:1-2-2-1 PartialL¨u [70] 4D:1-2-2-2 NoModified L¨u [71] 4D:2-2-2-1 NoWang-Liu [72] 4D:1-2-2-1 PartialJia [73,74] 4D:1-2-2-2 PartialQWWC system [75,76] 4D:2-2-2-2 No
Table 2
List of hyperchaotic systems and their reductionproperties. The labeling is the same as was given in Table1. We note that specific forms of some equations can changefrom paper to paper, so we record the specific equations usedin Appendix B. chaos [77], as they each involve only a single quadraticnonlinearity. On the other hand, systems with strongerpolynomial nonlinearities, or systems with many non-linear terms, appear to often be irreducible under dif-ferential elimination. Note that for cases where the re-duction might involve integrals, resulting in a type ofintegro-differential equation, the differential eliminationalgorithm would miss such a reduction, even though itexists. This is due to the fact that the differential elim-ination algorithm is working over the ring of deriva-tions, which does not include integrals. Indeed, since
Heather A. Harrington and Robert A. Van Gorder integral operators are fairly cumbersome to introducecompared to their differential operator counterparts (wediscuss this later in Section 7), obtaining an algorith-mic approach including integrals would be challenging.Therefore, the differential elimination algorithm out-lined in Section 2 appears to be a very useful tool forreducing the dimension of dynamical systems, providedthat a reduction to a single ODE exists. For the morecomplicated models, we find the need to proceed on acase-by-case basis with manual manipulations due toany integration needed.We demonstrate reduction of dimension for chaoticsystems into single higher-order ODEs in the next threesubsections. We pick a case where all state variables canbe isolated (the R¨ossler system), a case where one of thestate variables can be isolated in terms of a differentialequation (the Lorenz system), and finally a case wherenone of the state variables can be isolated in terms of adifferential equation (the Qi-Chen-Du-Chen-Yuan sys-tem) so that any reduction would necessarily involve in-tegrals. For all cases considered, we let x, y, z ∈ C n ( R )where n is the dimension of the relevant dynamical sys-tem, and we take a, b, c ∈ R to be parameters.3.1 R¨ossler SystemThe R¨ossler equation [44,45] reads˙ x = − y − z , ˙ y = x + ay , ˙ z = b + x ( x − c ) . (3)We first obtain the ODE for y ( t ). Note from thesecond equation that x = ˙ y − ay , so that ˙ y = ¨ y − a ˙ y andhence from the first equation we have z = − ¨ y + a ˙ y − y .Placing these into the third equation, and performingalgebraic manipulations, we obtain... y − a ¨ y + ˙ y − (¨ y − a ˙ y + y ) ( ˙ y − ay − c ) + b = 0 . (4)Note that this equation is third order, and therefore theinformation of the three-dimensional system (3) can beencoded in this single ODE. By similar manipulations,one may arrive at an equation for x ( t ),( a + c − x ) (cid:26)(cid:18) ddt − ( x − c ) (cid:19) ¨ x − a ˙ x + x + ba + c − x − b (cid:27) = 0 , (5)and an equation for z ( t ), z (cid:18) d dt − a ddt + c (cid:19) ˙ z − bz + z ˙ z − az + cz = 0 . (6) 3.2 Lorenz SystemThe Lorenz system [34,35] is given by˙ x = a ( y − x ) , ˙ y = x ( b − z ) − y , ˙ z = xy − cz . (7)Observing from the first two equations that y = x + 1 a ˙ x (8)and z = b − ¨ x + (1 + a ) ˙ x + xax , (9)the third equation can be used to obtain a single ODEfor the state variable x ( t ), viz., x (cid:18) ddt + c (cid:19) ¨ x + (1 + a ) ˙ x + xx + ax + x ˙ x − abcx = 0 . (10)This agrees with what one obtains from the differentialelimination. On the other hand, we observe that thealgorithmic approach to differential elimination is use-ful for situations in which there is no obvious route toreduce a system into a single equation (through elim-inations and substitutions). A good example of this isfound when trying to obtain a differential equation forthe state variable z ( t ) alone. Using the differential elim-ination, we arrive at a rather complicated equation ofthe form( b − z )(... z ) + P ( z, ˙ z, ¨ z )... z + P ( z, ˙ z, ¨ z ) = 0 , (11)where P and P are complicated polynomials that wedo not list for sake of brevity. Interestingly, this is a fullynonlinear equation, since the highest order derivativeenters into the equation nonlinearly. In contrast, theequation obtained for the state variable x ( t ) is quasi-linear, since it is linear in the highest derivative. Onecould differentiate the equation for z ( t ) in order to iso-late the highest derivative, but by doing so one wouldincrease the differential order of the system, therebydecreasing the regularity of the system. This is partic-ularly important in cases where the solution z ( t ) mayonly be C ( R ).When a system is nonlinear, there may of course beforms of the nonlinearity which do not permit one toobtain an equation for a single state variable in terms ofthat state variable and its derivatives. A good exampleof this is the state variable y ( t ) in the Lorenz system.The algorithmic differential elimination finds no closed eduction of dimension for nonlinear dynamical systems 7 differential equation for y ( t ). As it turns out, the rea-son for this is that any equation governing y ( t ) alonewill necessarily involve integral terms which cannot beeliminated (due to the nonlinearity of the equation). Tosee this, note that if we consider the first equation inthe Lorenz system, which may be written in the form( e at x ) ′ = ae at y , we find x ( t ) = x + ae − at Z t e as y ( s ) ds . (12)Here x is the initial value of the state x ( t ), that is x (0) = x . Yet, from the second equation in the Lorenzsystem, we have z = b − ( ˙ y + y ) /x , which yields z ( t ) = b − ˙ y + yx + ae − at R t e as y ( s ) ds . (13)Placing the representations for x ( t ) and z ( t ) into thethird equation of the Lorenz system, and performing al-gebraic manipulations to simplify the resulting expres-sion, we obtain(¨ y + (1 + c ) ˙ y + cy ) (cid:18) x + ae − at Z t e as y ( s ) ds (cid:19) + ( ˙ y + y ) (cid:18) ay − a e − at Z t e as y ( s ) ds (cid:19) + y (cid:18) x + ae − at Z t e as y ( s ) ds (cid:19) − cb (cid:18) x + ae − at Z t e as y ( s ) ds (cid:19) = 0 . (14)Note that the equation both involves an integral and isnon-autonomous.3.3 Qi-Chen-Du-Chen-Yuan SystemWe now consider the Qi-Chen-Du-Chen-Yuan (QCDCY)system [54], which is given by˙ x = a ( y − x ) + yz , ˙ y = bx − y − xz , ˙ z = xy − cz . (15)The differential elimination algorithm indicates there isno reduction to a single ODE in any of the three statevariables. This system has a quadratic nonlinearity ineach equation, and this added complication is behindthe difficulties in obtaining such a reduction. However,we may still obtain an equation for a single state vari-able, if we are willing to include integral terms. Due tothe complexity in obtaining such an equation, we shallrestrict our attention to finding a single equation forthe state variable z ( t ), noting that similar approaches can be used to find a single equation for either of theother two state variables, x ( t ) or y ( t ).Let us begin by noting that the second equation inthe QCDCY system implies ( e t y ) ′ = e t ( b − z ) x , whileplacing this into the third equation in the QCDCY sys-tem gives ( a + z )( ˙ z + cz ) = xe − at ( e at x ) ′ = x ( ˙ x + ax ).This, in turn, implies that state variables x ( t ) and z ( t )satisfy e at ( x ( t )) = x + 2 Z t e as ( a + z ( s ))( ˙ z ( s ) + cz ( s )) ds . (16)where x (0) = x . The first equation in the QCDCYsystem has not been used, and we place this relationinto that equation to obtain a single equation for thestate variable z ( t ). After several algebraic and differ-ential manipulations, we arrive at the single equation2 e at (cid:18) − a − ˙ za + z (cid:19) ( ˙ z + cz )( a + z ) (cid:0) x + J [ z, ˙ z ] (cid:1) + 2 (cid:0) x + J [ z, ˙ z ] (cid:1) ddt (cid:0) e at ( ˙ z + cz )( a + z ) (cid:1) − e at ( ˙ z + cz ) ( a + z ) − b − z )( a + z ) (cid:0) x + J [ z, ˙ z ] s (cid:1) = 0 , (17)where we have defined the integral operator J [ z, ˙ z ] = 2 Z t e as ( a + z ( s ))( ˙ z ( s ) + cz ( s )) ds . (18)Similar results can be obtained for the other statevariables. The fact that the obtained equations involvean integral operator which cannot simply be differen-tiated away demonstrates why the differential elimina-tion algorithm was not useful for this case. Still, per-forming the manipulations by hand, we have reducedthe fairly complicated QCDCY system into a singleintegro-differential equation, thereby reducing the di-mension of the original system. n -dimensional dynamicalsystems We now give some summarizing remarks based on whatwe have seen in the previous sections. We shall assumethat each system is coupled through at least one statevariable (otherwise the state variables naturally sepa-rate into distinct lower-dimension equations, and theapproach is not needed).
Heather A. Harrington and Robert A. Van Gorder n , there isalways a reduction into a single higher-order ODE. Thisfollows from the process of Gaussian elimination. In thecase that the matrix of coefficients for such a first-ordersystem is full rank, the resulting higher-order ODE willbe of order n . If the matrix of coefficients is singular,then the resulting higher-order ODE will be of orderless than n .4.2 Reducible Nonlinear SystemsFor first order nonlinear systems of dimension n , thereare multiple possibilities, owing to the structure of thenonlinearity.In cases where the system permits the completedifferential elimination (an example being the Rosslerequation), all state variables in a first order nonlin-ear system can be expressed in terms of a higher-orderODE. Note, however, that it is possible for the orderof the single ODE to be different from the dimension n of the first order system. As an example of this point,consider the system˙ x = x − y − z , ˙ y = x , ˙ z = x − x . (19)Clearly, differentiation of the first equation gives ¨ x =˙ x − ˙ y − ˙ z = ˙ x − x + x − x . So, we obtain¨ x − ˙ x − x + x + x = 0 , (20)which is a second order equation for the state variable x ( t ), even though the original system was first order. Asimilar example can be found in [1], where a fourth or-der nonlinear dynamical system was reduced to a singlesecond order nonlinear ODE.It is possible for a system to be reduced to a singleequation, which is not an ODE. This was evident evenfor the Lorenz equation, where an equation for one ofthe state variables involves an integral term in additionto derivative terms. Note that the equation was notclosed under any number of differentiations, due to theform of the intergral terms. As such, the single reducedequation for the state variable could never be expressedstrictly as an ODE of any finite order. Note that thiscan occur for one of the state variables, while a dif-ferent state variable might satisfy a finite order ODE.For such cases, the nonlinearity in the system results intheir being certain favored state variables with whichto perform the reduction to a single ODE. 4.3 Differentially Irreducible Nonlinear SystemsWe have observed that for more complicated nonlin-ear dynamical systems, there is no reduction to a singleODE in one state variable. While it may be the casethat differential elimination does not pick up an ODEthat does exist, it seems as though the failure of dif-ferential elimination is a sign that integrations will beneeded in order to reduce the dimension of such sys-tems. Indeed, when integrations of this kind are calledfor, the manipulations are no longer confined to thespecified differential ring, and the differential elimina-tion cannot be performed. While one can attempt theseintegrations manually, as opposed to algorithmically,obviously it would be desirable to have some kind ofalgorithmic approach. Perhaps one may adjoin inte-gral operators to the differential ring, in order to per-form reductions for more complicated nonlinear sys-tems. This would likely work in cases like the Lorenzsystem, for which there is partial reducibility under dif-ferential elimination. For instance, if one were to definea new variable Y ( t ) = R t e as y ( s ) ds , then one wouldobtain a non-autonomous ODE for Y ( t ) from equation(14). Therefore, this fairly simple integral transforma-tion, in addition to differential operators, can reducethe dimension of the Lorenz system with respect tothe state variable y ( t ). However, in cases like that ofthe QCDCY system, note that the form of the integraloperator given in (18) is rather complicated, depend-ing nonlinearly on both the state variable z ( t ) and itsderivative ˙ z ( t ). For such cases, there is no combinationof elementary integral transforms that can be adjoinedto the differential ring which would permit reductionof dimension to a single ODE. As such, it appears asthough reduction of dimension for certain more com-plicated systems will result in reductions to integro-differential equations, rather than ODEs, for some fun-damental reason related to how complicated the orig-inal dynamical system is. Therefore, the study of pos-sible algorithmic methods for the reduction of dimen-sion for dynamical systems into single scalar integro-differential equations would be an interesting and po-tentially very useful area of future work. Turning out attention now to practical applications forreduction of dimension, recall that contraction mapsand Lyapunov functions are useful tools for studyingthe asymptotic stability of nonlinear dynamical sys-tems. In this section, we use the three examples workedexplicitly in Section 3 in order to demonstrate the util-ity of reduction of dimension for finding Lyapunov func- eduction of dimension for nonlinear dynamical systems 9 tions. Using these results, we can recover stability re-sults for these dynamical systems which were obtainedthrough other approaches, and which agree with exist-ing results in the literature.5.1 R¨ossler SystemThe R¨ossler system has two equilibrium values, ± y ∗ ,for y ( t ), and the constant y ∗ must satisfy the quadraticequation a ( y ∗ ) + cy ∗ + b = 0 . (21)In order to discover a Lyapunov function for theR¨ossler system, it is tempting to assume a bowl-shapedmap of the form αx + βy + γz , or minor variationson this theme involving higher power polynomials ofeven order, but the approach evidently proves fruitless.Therefore, we shall use one of the three equations ob-tained for the isolated state variables of the R¨osslersystem.Consider equation (4) for the R¨ossler system (3).Let us write Y ( t ) = y ( t ) − y ∗ in the neighborhood ofeither equilibrium value y ∗ . This transformation willprove useful, as the Lyapunov function needs to vanishat the equilibrium value selected. (There is therefore theneed to construct such a function in a neighborhood ofeach equilibrium point.) Under this transformation, (4)is put into the form... Y − a ¨ Y + ˙ Y − (cid:16) ¨ Y − a ˙ Y + Y (cid:17) (cid:16) ˙ Y − aY (cid:17) = 0 . (22)Let us define a function m = ¨ Y − a ˙ Y + Y so that (22)is put into the form˙ m − ( ˙ Y − aY ) m = 0 . (23)Observe that (23) can be written as˙ m − e at (cid:0) e − at Y (cid:1) ′ m = 0 . (24)From this, we recover¨ Y − a ˙ Y + Y = m = m exp (cid:18)Z t e aζ ( e − aζ Y ( ζ )) ′ dζ (cid:19) , (25)where m is a constant of integration. As we are inter-ested in recovering information about the asymptoticstability of the R¨ossler system, let us pick the initialcondition Y (0) = ǫ . This corresponds to setting the ini-tial condition such that it is contained within a neigh-borhood of the equilibrium value. Let us also restrict | a | < Y ( t ) = ǫe at/ ( cos √ − a t ! + C sin √ − a t !) + m Z t K ( t, s ) exp (cid:18)Z s e aζ ( e − aζ Y ( ζ )) ′ dζ (cid:19) ds , (26)where C is a constant that will depend on the initialvalue of ˙ Y (0) (the value of which will not impact ouranalysis) and K ( t, s ) is the kernel K ( t, s ) = e a ( t − s ) sin (cid:18) − a t − s ) (cid:19) . (27)Observe that for − < a < ǫ >
0, for large enough time˜ t ( ǫ ) >
0, the solution Y ( t ) will lie in a neighborhood − ǫ < Y ( t ) < ǫ for all t > ˜ t ( ǫ ). Therefore, Y → t → ∞ . Yet, by definition of Y ( t ), this implies y → y ∗ as t → ∞ . Using this, one may shown that x →− ay ∗ and z → − y ∗ as t → ∞ . Hence, we have shownthat a < y ∗ found from (21). Indeed, wedid not even need to calculate the equilibrium value y ∗ for the present analysis, as the analysis holds for anarbitrary equilibrium value satisfying (21).5.2 Lorenz SystemIn order to find a Lyapunov function for the Lorenz sys-tem in a neighborhood of the zero equilibrium ( x, y, z ) =(0 , , V ( x, y, z ) = αx + βy + γz , (28)where α > β >
0, and γ > a , b , and c are positive. Then, thetime derivative of V is given by12 ˙ V = − αax − βy − γcz + ( αa + βb ) xy + ( γ − βb ) xyz . (29)Clearly, we should take γ = βb . Note that − ( √ αax − p βy ) = − αax − βy + 2 p αβaxy . (30) Then,12 ˙ V = − ( √ αax − p βy ) − γcz +( αa + βb − p αβa ) xy , (31)hence ˙ V ≤ βb < √ αβa − αa (sincethis would imply − αax − βy + ( αa + βb ) xy < β = αa . Then, the condition reduces to b < α > α = . This means thatwhenever a >
0, 0 < b <
1, and c >
0, there exists aLyapunov function V ( x, y, z ) = 12 x + a y + ab z , (32)since V (0 , ,
0) = 0, | V | → ∞ as | ( x, y, z ) | → ∞ (ra-dially unbounded), and ˙ V < x, y, z ) = (0 , , < b < V ( x, ˙ x, ¨ x ) = 12 x + a (cid:18) x + ˙ xa (cid:19) + ab (cid:18) b − ¨ x + (1 + a ) ˙ x + axax (cid:19) (33)for the state variable x ( t ). Then, one may verify ˙ˆ V < x = 0. One can obtainsimilar contraction maps in either of the other two statevariables.5.3 Qi-Chen-Du-Chen-Yuan SystemIn order to find a contraction map for the Qi-Chen-Du-Chen-Yuan (QCDCY) system, we begin with the bowlshaped assumption for a Lyapunov function about theequilibrium ( x, y, z ) = (0 , , V ( x, y, z ) = αx + βy + γz , (34)where α > β >
0, and γ > t andusing the three constituent equations of the QCDCYsystem, we have˙ V = ( α − β + γ ) xyz + ( αa + βb ) xy − αax − βy − γcz . (35)Since we need α > β >
0, and γ >
0, we shouldconsider model parameters satisfying a > c >
0. To remove the first term, which is hyperbolic in nature,we should choose β = α + γ . Meanwhile, the removethe second term, which is also hyperbolic, we shouldset β = − ab α for non-zero b . Since all other parametersare positive, we must require b <
0. Then, β = a | b | α ,and placing this into β = α + γ gives γ = a −| b || b | . Aswe need γ >
0, this gives the added restriction a > | b | .The parameter α is arbitrary, so we take α = . Wetherefore obtain V ( x, y, z ) = 12 x + a | b | y + a − | b | | b | z , (36)and this candidate function is indeed a contraction mapsatisfying V (0 , ,
0) = 0, ˙
V < x, y, z ) = (0 , , | V | → ∞ as | ( x, y, z ) | → ∞ , provided that theparameter restrictions a > | b | , b <
0, and c > z ( t ) in the QCDCY sys-tem in Section 3, we shall choose to construct a con-traction map for that state variable here. Doing so, wefind thatˆ V ( z, ˙ z ) = a | b | e at ( ˙ z + cz ) x + 2 R t e as ( a + z ( s ))( ˙ z ( s ) + cz ( s )) ds + e − at (cid:18) x + 2 Z t e as ( a + z ( s ))( ˙ z ( s ) + cz ( s )) ds (cid:19) + a − | b || b | z (37)satisfies ˙ˆ V < z = 0, given that a > | b | , b < c >
0. Hence, ˆ V ( z, ˙ z ) is a contraction map for thestate variable z ( t ) when a > | b | , b <
0, and c >
There are a variety of methods available for trying tofind chaotic trajectories in nonlinear dynamical sys-tems, and the approach highlighted in this paper doesnot add to collection of tools, explicitly. However, thereduction of dimension approach outlined in Section 2can be used to make finding chaos in dynamical systemsmore efficient. To demonstrate this, we shall consider eduction of dimension for nonlinear dynamical systems 11 two rather distinct approaches, namely, the undeter-mined coefficients method for obtaining chaotic trajec-tories and the competitive modes analysis for identifi-cation of chaotic parameter regimes. For each of theseapproaches, we show that an application of reductionof dimension results in a simplification of each test forchaos.6.1 Calculation of trajectories via undeterminedcoefficientsWhen attempting to analytically calculate chaotic tra-jectories, even in an approximate sense, one often re-duces the dimension of the governing equations. Thereason for this lies in the fact that it is easier to con-sider an expansion for one state variable, rather thanmultiple state variables. To best illustrate this point,let us return to the R¨ossler equation (3).One popular method for approximating trajectoriesof chaotic systems analytically is the undetermined co-efficient method [3,4,5,6]. Since Taylor series expan-sions for nonlinear systems often have a finite regionof convergence centered at the origin, yet the chaoticdynamics remain bounded in space, one often consid-ers non-polynomial base functions. One popular choicewould be a function of the form S ( t ; { A j } j = ∞ j = −∞ , α ) = (P ∞ j =0 A j e − αjt for t ≥ , P ∞ j =0 A − j e αjt for t < . (38)In this expression, the A j ∈ R and the parameter α > x ( t ) = S ( t ; { A j } j = ∞ j = −∞ , α ), y ( t ) = S ( t ; { B j } j = ∞ j = −∞ , α ), and z ( t ) = S ( t ; { C j } j = ∞ j = −∞ , α ).Placing these equations into (3), one would obtain aninfinite system of nonlinear algebraic equations for allof the coefficients and the temporal scaling α >
0. Inpractice, one would truncate these expansions, takingthe sum over − J ≤ j ≤ J for some J >
0. As thesolutions may converge slowly - if they converge at all(owing to the nonlinearity), one would need to solve6 J + 1 nonlinear algebraic equations.Assume, instead, that we wish to solve (4) by theapproach described above. We would then insert theexpansion for y ( t ) into (4). Assuming that we can solvethe resulting nonlinear algebraic equations for the con-stants { B j } j = ∞ j = −∞ and α , we can then recover x ( t ) and z ( t ) by recalling x = ˙ y − ay and z = − ¨ y + a ˙ y − y . Fromthese expressions, it is simple to show A n = − ( α | n | + a ) B n and C n = − ( α n + aα | n | + 1) B n for all n ∈ Z . If we were to truncate the expansion for y ( t ), inthe manner described above, we would need to solve2 J + 1 nonlinear algebraic equations for { B j } j = Jj = − J and α , while the coefficients for x ( t ) and z ( t ) are immedi-ately found once we know these parameters. This meansthat by first reducing the dimension of the ODE sys-tem, we would be able to reduce the computationalcomplexity of the problem by a factor of three. Forhigher-dimensional system, the reduction in computa-tional complexity will scale as the dimension of the sys-tem itself. In other words, a solution in term of the un-determined coefficient method will not depend on thesize of the dynamical system provided that the dynam-ical system can be reduced in dimension to a singleequation governing one state variable.6.2 Competitive modes analysis: A check for chaosThe method of competitive modes involves recastinga dynamical system as a coupled system of oscillators[13,14,15,16,17,18]. Consider the general nonlinear au-tonomous system of dimension n given by˙ x i = f i ( x , x , ..., x n ) . (39)Differentiation of (39) once gives a coupled system ofsecond order equations,¨ x i = n X j =1 f j ∂f i ∂x j = − x i g i ( x , x , ..., x i , ..., x n )+ h i ( x , x , ..., x i − , x i +1 , ..., x n ) . (40)When a g i is positive, its respective i th equation be-haves like an oscillator. The following conjecture is posedin [14]: Competitive Modes Requirements : The conditionsfor dynamical systems to be chaotic are given by:(A) there exist at least two modes, labeled g i in thesystem;(B) at least two g ’s are competitive or nearly competi-tive, that is, for some i and j , g i ≈ g j > t ;(C) at least one of the g ’s is a function of evolutionvariables such as t ; and(D) at least one of the h ’s is a function of system vari-ables.The requirements (A)-(D) essentially tell us that acondition for chaos is that two or more equations in (40) behave as oscillators ( g i > g i = g j , i = j , i, j = 1 , , , . . . , n . Accounting for symmetry, this gives2 n − − p > p need not be equal to n , as we have seen in earlier sections). Then, associating y ( t ) to this single state variable, we have d p ydt p = F (cid:18) y, dydt , . . . , d p − ydt p − (cid:19) . (41)Since the competitive modes analysis relies on us ob-taining a system of oscillator equations, let take y = y and write the equation (41) as the system˙ y = y . . . , ˙ y p − = y p , ˙ y p = F ( y , . . . , y p ) . (42)Differentiation of (42) results in the system of secondorder equations given by¨ y = y , ...¨ y p − = y p , ¨ y p − = ˙ y p = F ( y , . . . , y p ) , ¨ y p = p X i =1 ∂F∂y i ˙ y i = p − X i =1 ∂F∂y i y i +1 + ∂F∂y p F ( y , . . . , y p ) . (43)The right hand side of the first p − g = · · · = g p − = 0. Hence, these equations are never oscillators. Meanwhile, we can de-compose the right hand sides of the latter two equa-tions, so that F ( y , . . . , y p ) = − y p − g p − + h p − (44)and p − X i =1 ∂F∂y i y i +1 + ∂F∂y p F ( y , . . . , y p ) = − y p g p + h p . (45)Note that there are now exactly two mode frequen-cies, g p − and g p , and there is always an h i depend-ing on state variables. In order to determine if the sys-tem (42) satisfies the competitiveness conditions (A)-(D) (and therefore if the original system satisfies thesecompetitiveness conditions), it is sufficient to check if g p − = g p > n − − The construction of Lyapunov functions for nonlineardynamical systems is often either simple, or quite chal-lenging, with little room in between. Aside from choos-ing some standard forms (such as the common bowlshape centered about an equilibrium value), there ismore an art to the selection of such function. However,as we have demonstrated for the R¨ossler system, it ispossible to use differential elimination to obtain a con-traction map in a single state variable, which can thenbe used to obtain a stability result for all state vari-ables. A similar approach was also employed in [1] tostudy Michaelis-Menten enzymatic reactions, and afterusing a reduced equation, the proof of global asymp-totic stability for a positive steady state was rathersimple. Differential elimination, and reduction of di-mension more generally, can therefore be useful in help-ing one determine asymptotic properties of solutions tononlinear dynamical systems. On the other hand, it isnow known that there are autonomous chaotic systemswhich lack any equilibrium [78]. The approach outlinedhere could still be used to construct maps which demon-strate the boundedness of trajectories in time for suchsystems, even if such trajectories do not approach anyfixed points.While useful for obtaining contraction maps, reduc-tion of dimension is also a promising tool for finding eduction of dimension for nonlinear dynamical systems 13 and studying chaotic trajectories in nonlinear dynam-ical systems. We were able to show that reduction ofdimension can be used to simplify calculations involvedin using undetermined coefficient methods [3,4,5,6] bya factor of 1 /n , where n is the dimension of the dy-namical system, meaning that the number of equationsneeding to be solved will not increase with the size ofthe system - i.e. the computational complexity of theapproach will not scale with the size of the system butrather will remain fixed. This means that one may ap-proximate chaotic trajectories through such approachesin systems of rather large dimension, without the com-putational problem becoming unwieldy, provided thatthe system can be reduced to a single ODE governingonly one state variable.While chaos is often studied numerically, recentlyanalytical approaches have been employed to constructtrajectories approximating chaotic orbits. When em-ploying these methods, it can be beneficial to considera single equation rather than a system of equations,even if the single equation is more complicated. This istrue for series and perturbation approaches, as the re-duction of order requires one to track fewer functions,which is particularly useful when dealing with messyequations. Furthermore, analytical approaches permit-ting the control of error, such as the optimal homotopyanalysis method, rely on assigning an error control pa-rameter to each state variable. Reduction of dimensioncan allow one to minimize error via a single control pa-rameter [2], rather than over multiple parameters [79],which is computationally less demanding.The reduction of dimension can also be useful for di-agnostic tests for chaos. When applying the competitivemodes analysis, a type of diagnostic criteria for findingchaotic in nonlinear dynamical systems, one performsbinary comparisons between the mode frequencies ofan oscillator corresponding to each equation. However,if one first applies reduction of dimension and reducesthe system of a single ODE for one state variable, weprove that only one comparison would needed. This isparticularly beneficial when studying large dimensionalsystems, since the number of naive comparisons neededscales like 2 n − in dimension n .In the future, it would be interesting to consideran algorithm that considered elimination not only ele-ments ∂ n ( ∂ = ddt ) of a differential ring R [[ ∂ ]], but alsointegral operators ∂ − n (where ∂ − n satisfies ∂ − n ∂ n = ∂ and hence is the inversion of the operator ∂ n ). Indeed,in cases where the differential elimination algorithmfailed to give a reduction of the system to a single ODE,we found by manual substitutions that one can arriveat an integro-differential equation. While more compli-cated, such integro-differential equations can still cast light on the behavior of solutions, and can prove usefulin obtaining Lyapunov functions. More generally thanfor dynamical systems, these inverse operators ∂ − n playa role in the study of operators and integrable hierar-chies arising in nonlinear evolution PDEs [80,81,82,83,84,85,86]. Therefore, the extension of the algorithm tothe ring of formal Laurent series in ∂ would be a fruit-ful area for future work, not only for dynamical systemsbut also for integrable partial differential equations.AcknowledgementsHAH gratefully acknowledges the support of EPSRCFellowship EP/K041096/1. References
1. K. Mallory and R. A. Van Gorder, A transformed time-dependent Michaelis-Menten enzymatic reaction modeland its asymptotic stability, Journal of MathematicalChemistry 52, 222-230 (2014).2. R. Li, R. A. Van Gorder, K. Mallory and K. Vajrav-elu, Solution method for the transformed time-dependentMichaelis-Menten enzymatic reaction model, Journal ofMathematical Chemistry 52, 2494-2506 (2014).3. T. Zhou, G. Chen, and Q. Yang, Constructing a newchaotic system based on the Silnikov criterion, Chaos,Solitons & Fractals 19, 985-993 (2004).4. Y. Jiang and J. Sun, Si’lnikov homoclinic orbits in a newchaotic system, Chaos, Solitons & Fractals 32, 150-159(2007).5. R.A. Van Gorder and S.R. Choudhury, Shil’nikov Anal-ysis of Homoclinic and Heteroclinic Orbits of the T Sys-tem, Journal of Computational and Nonlinear Dynamics6, 021013 (2011).6. T. Zhou and G. Chen, Classification of chaos in 3-Dautonomous quadratic systems-I: basic framework andmethods, International Journal of Bifurcation and Chaos16(09), 2459-2479 (2006).7. S.J. Liao, Beyond Perturbation: Introduction to the Ho-motopy Analysis Method, Chapman & Hall/CRC Press,Boca Raton, 2003.8. S.J. Liao, Homotopy Analysis Method in Nonlinear Dif-ferential Equations. Springer & Higher Education Press,Heidelberg, 2012.9. S. Liao, An optimal homotopy-analysis approach forstrongly nonlinear differential equations, Communica-tions in Nonlinear Science and Numerical Simulation 15,2315-2332 (2010).10. C.K. Ahn, Neural Network H ∞ Chaos Synchronization,Nonlinear Dynamics 60, 295-302 (2010).11. C.K. Ahn, S.T. Jung, S.K. Kang, and S.C. Joo, AdaptiveH ∞ Synchronization for Uncertain Chaotic Systems withExternal Disturbance, Communications in Nonlinear Sci-ence and Numerical Simulation 15, 2168-2177 (2010).12. C.K. Ahn, Output Feedback H ∞ Synchronization for De-layed Chaotic Neural Networks, Nonlinear Dynamics 59,319-327 (2010).13. S. R. Choudhury and R. A. Van Gorder, Competitivemodes as reliable predictors of chaos versus hyperchaosand as geometric mappings accurately delimiting attrac-tors, Nonlinear Dynamics 69, 2255 (2012).4 Heather A. Harrington and Robert A. Van Gorder14. W. Yao, P. Yu, and C. Essex, Estimation of chaotic pa-rameter regimes via generalized competitive mode ap-proach, Communications in Nonlinear Science and Nu-merical Simulation 7, 197 (2002).15. P. Yu, W. Yao, and G. Chen, Analysis on topologicalproperties of the Lorenz and the Chen attractors usingGCM, International Journal of Bifurcation and Chaos 17,2791 (2007).16. Z. Chen, Z.Q. Wu, and P. Yu, The critical phenomenain a hysteretic model due to the interaction between hys-teretic damping and external force, Journal of Sound andVibration 284, 783 (2005).17. W. Yao, P. Yu, C. Essex, and M. Davison, Competitivemodes and their application, International Journal of Bi-furcation and Chaos 16, 497 (2006).18. R.A. Van Gorder and S.R. Choudhury, Classification ofchaotic regimes in the T system by use of competitivemodes, International Journal of Bifurcation and Chaos20, 3785-3793 (2010).19. K. Mallory and R.A. Van Gorder, Competitive modes forthe detection of chaotic parameter regimes in the gen-eral chaotic bilinear system of Lorenz type, InternationalJournal of Bifurcation and Chaos 25, 1530012 (2015).20. J. F. Ritt. Differential algebra. Dover Publications Inc.,New York, (1950).21. E. R. Kolchin. Differential Algebra and AlgebraicGroups. Academic Press, New York (1973).22. F. Boulier, D. Lazard, F. Ollivier and M. Petitot. Com-puting representations for radicals of finitely generateddifferential ideals, AAECC 20(1), 73-121 (2009).23. F. Boulier and F Lemaire. Computing canonical repre-sentatives of regular differential ideals. ISSAC (2000)24. F. Boulier, D. Lazard, F. Ollivier, M. Petitot. Represen-tation for the radical of a finitely generated differentialideal. ISSAC. (1995)25. A. Rosenfeld. Specializations in differential algebra.Transactions of the American Mathematical Society 90,394-407 (1959).26. A. Hashemi and Z. Touraji. An Improvement ofRosenfeld-Gr¨obner Algorithm. In H. Hong and C. Yap(Eds): ICMS Springer-Verlag Berlin Heidelberg. LNCS8592, 466–471 (2014) .27. A. Seidenberg. An elimination theory for differential al-gebra. Univ. California Publ. Math. (New Series). (1956).28. F. Boulier. Differential Elimination and Biological Model-lling. M. Rosenkranz and D. Wang. Workshop D2.2 of theSpecial Semester on Gr¨obner Bases and Related Meth-ods, May 2006, Hagenberg, Austria. de Gruyter, RadonSeries Comp. Appl. Math. 2, 111-139 (2007).29. N. Meshkat, C. Anderson, and J. J. DiStefano III. Al-ternative to Ritt’s pseudo division for finding the input-output equations of multi-output models. MathematicalBiosciences 239, 117–123 (2012).30. B. Buchberger. An Algorithm for Finding a Basis for theResidue Class Ring of a Zero-Dimensional PolynomialIdeal (German). PhD thesis, Math. Inst. Univ of Inns-bruck, Austria (1965).31. O. Zariski and P. Samuel. Commutative Algebra. Grad-uate Texts in Mathematics, vol 28, Springer Verlag, NewYork (1958).32. F. Boulier. The BLAD libraries. ∼ boulier/BLAD (2004).33. W.A. Stein et al. Sage Mathematics Software. The SageDevelopment Team. (2015).34. E.N. Lorenz, Deterministic nonperiodic flow, Journal ofthe Atmospheric Sciences 20(2), 130-141 (1963). 35. M.W. Hirsch, S. Smale, Stephen, and R. Devaney, Dif-ferential Equations, Dynamical Systems, & An Introduc-tion to Chaos (Second ed.). Boston, MA: Academic Press.(2003) ISBN 978-0-12-349703-1.36. L.O. Chua, T. Matsumoto, and M. Komuro, The Dou-ble Scroll, IEEE Transactions on Circuits and Systems(IEEE). CAS-32 (8), 798-818 (1985).37. C.C. Hwang, H.-Y. Chow, and Y.-K. Wang, A new feed-back control of a modified Chua’s circuit system, PhysicaD: Nonlinear Phenomena 92, 95-100 (1996).38. M.T. Yassen, Adaptive control and synchronization ofa modified Chua’s circuit system, Applied Mathematicsand Computation 135(1), 113-128 (2003).39. H.K. Chen and C.I. Lee, Anti-control of chaos in rigidbody motion Chaos, Solitons & Fractals, 21, 957-965(2004).40. C.-M. Chang and H.K. Chen, Chaos and hybrid pro-jective synchronization of commensurate and incommen-surate fractional-order ChenLee systems, Nonlinear Dy-namics 62, 851-858 (2010).41. M.I. Rabinovich and A.L. Fabrikant, On the stochasticself-modulation of waves in nonequilibrium media, Sov.Phys. JETP 77, 617-629 (1979).42. M.F. Danca and G. Chen, Bifurcation and Chaos ina Complex Model of Dissipative Medium, InternationalJournal of Bifurcation and Chaos 14(10), 3409-3447(2004).43. P. Grassberger and I. Procaccia, Measuring thestrangeness of strange attractors, Physica D 9, 189-208(1983).44. O.E. R¨ossler, An Equation for Continuous Chaos, PhysicsLetters A 57(5), 397-398 (1976).45. P. Zgliczynski, Computer assisted proof of chaos in theR¨ossler equations and in the H´enon map, Nonlinearity10, 243-252 (1997).46. G. Chen and T. Ueta, Yet another chaotic attractor, In-ternational Journal of Bifurcation and Chaos 9, 1465-1466 (1999).47. T. Ueta and G. Chen, Bifurcation analysis of Chen’sequation, International Journal of Bifurcation and Chaos10(08), 1917-1931 (2000).48. J. L¨u and G. Chen, A new chaotic attractor coined, Inter-national Journal of Bifurcation and Chaos 12(03), 659-661 (2002).49. J. L¨u, G. Chen, D. Cheng, and S. Celikovsky, Bridge thegap between the Lorenz system and the Chen system,International Journal of Bifurcation and Chaos 12(12),2917-2926 (2002).50. G. Tigan, Analysis of a dynamical system derived fromthe Lorenz system, Sci. Bull. Politehnica Univ Timisoara50, 61-72 (2005).51. G. Tigan and D. Opris, Analysis of a 3D chaotic system,Chaos, Solitons & Fractals 36(5), 1315-1319 (2008).52. R.A. Van Gorder and S.R. Choudhury, Analytical Hopfbifurcation and stability analysis of T system, Commu-nications in Theoretical Physics 55(4), 609-616 (2011).53. G. Qi, S. Du, G. Chen, Z. Chen, and Z. Yuan, On a four-dimensional chaotic system, Chaos, Solitons & Fractals23(5), 1671-1682 (2005).54. G. Qi, G. Chen, S. Du, Z. Chen, and Z. Yuan, Analysis ofa new chaotic system, Physica A: Statistical Mechanicsand its Applications 352(2), 295-308 (2005).55. S. Celikovsky and G. Chen, On a generalized Lorenzcanonical form of chaotic systems, International Journalof Bifurcation and Chaos 12(08), 1789-1812 (2002).56. S. Celikovsky and G. Chen, On the generalized Lorenzcanonical form, Chaos, Solitons & Fractals 26(5), 1271-1276 (2005).eduction of dimension for nonlinear dynamical systems 1557. R. A. Van Gorder, Emergence of chaotic regimes in thegeneralized Lorenz canonical form: a competitive modesanalysis, Nonlinear Dynamics 66, 153 (2011).58. D.V. Turaev and L.P. Shilnikov, Blue sky catastrophes,Dokl. Math. 51, 404 (1995).59. A. Shilnikov and G. Cymbalyuk, Transition betweentonic spiking and bursting in a neuron model via theblue-sky catastrophe, Physical Review Letters 94, 048101(2005).60. R. A. Van Gorder, Triple mode alignment in a canonicalmodel of the blue-sky catastrophe, Nonlinear Dynamics73, 397-403 (2013).61. L. Stenflo, Generalized Lorenz equations for acoustic-gravity waves in the atmosphere, Physica Scripta 53(1),83-84 (1996).62. L. Stenflo, Nonlinear acoustic-gravity waves, Journal ofPlasma Physics 75, 841-847 (2009).63. R.A. Van Gorder, Shilnikov chaos in the 4D Lorenz-Stenflo system modeling the time evolution of nonlinearacoustic-gravity waves in a rotating atmosphere, Nonlin-ear Dynamics 72(4), 837-851 (2013).64. R. Genesio and A. Tesi, Harmonic balance methods forthe analysis of chaotic dynamics in nonlinear systems,Automatica 28(3), 531-548 (1992).65. S.S. Motsa, P. Dlamini, and M. Khumalo, A new mul-tistage spectral relaxation method for solving chaoticinitial value systems, Nonlinear Dynamics 72, 265-283(2013).66. A. Arneodo, P. Coullet, and C. Tresser, Possible newstrange attractors with spiral structure, Communicationsin Mathematical Physics 79, 573-579 (1981).67. O.E. R¨ossler, An Equation for Hyperchaos, Physics Let-ters A 71 155-157 (1979).68. T.G. Gao, Z.Q. Chen, and G. Chen, A hyper-chaos gener-ated from Chen’s system, International Journal of Mod-ern Physics C 17, 471-478 (2006).69. T.G. Gao, Z. Chen, and Z.Q. Chen, Analysis of the hyper-chaos generated from Chen’s system, Chaos Solitons &Fractals 39, 1849-1855 (2009).70. A. Chen, J. Lu, J. L¨u, and S. Yu, Generating hyper-chaotic L¨u attractor via state feedback control, PhysicaA 64, 103-110 (2006).71. G. Wang, X. Zhang, Y. Zheng, and Y. Li, A new modifiedhyperchaotic L¨u system, Physica A 371, 260-272 (2006).72. F.-Q. Wang and C.-X. Liu, Hyperchaos evolved fromthe Liu chaotic system, Chinese Physics 15(5), 963-968(2006).73. Q. Jia, Hyperchaos generated from the Lorenz chaoticsystem and its control, Physics Letters A 366(3), 217-222(2007).74. Q. Jia, Projective synchronization of a new hyperchaoticLorenz system, Physics Letters A 370(1), 40-45 (2007).75. G. Qi, M.A. van Wyk, B.J. van Wyk, and G. Chen, On anew hyperchaotic system, Physics Letters A 372(2), 124-136 (2008).76. G. Qi, M.A. van Wyk, B.J. van Wyk, and G. Chen, Anew hyperchaotic system and its circuit implementation,Chaos, Solitons & Fractals 40(5), 2544-2549 (2009).77. C. Letellier, E. Roulin, and O.E. R¨ossler, Inequivalenttopologies of chaos in simple equations, Chaos, Solitons& Fractals 28, 337-360 (2006).78. S. Jafari, J. C. Sprott, and M.R. Hashemi-Golpayegani,Elementary quadratic chaotic flows with no equilibria,Physics Letters A 377(9), 699-702 (2013).79. R.A. Van Gorder, Analytical method for the constructionof solutions to the F¨oppl - von Karman equations govern-ing deflections of a thin flat plate, International Journalof Non-Linear Mechanics 47, 1-6 (2012). 80. M. Baxter, S. R. Choudhury, and R. A. Van Gorder,Zero curvature representation, bi-Hamiltonian structure,and an integrable hierarchy for the Zakharov-Ito system,Journal of Mathematical Physics 56, 063503 (2015).81. P.G. Drazin and R.S. Johnson, Solitons: an introduc-tion. Vol. 2. Cambridge University Press, Cambridge,UK, 1989.82. W.X. Ma, The algebraic structures of isospectral Lax op-erators and applications to integrable equations, Journalof Physics A: Mathematical and General 25, 5329-5343(1992).83. W.X. Ma, Lax representations and Lax operator alge-bras of isospectral and nonisospectral hierarchies of evo-lution equations, Journal of Mathematical Physics 33,2464-2476 (1992).84. W.X. Ma and M. Chen, Hamiltonian and quasi-Hamiltonian structures associated with semi-direct sumsof Lie algebras, Journal of Physics A: Mathematical andGeneral 39, 10787-10801 (2006).85. D. Robertz, Formal algorithmic elimination for PDEs.Vol. 2121. Springer, Cham, Switzerland 2014.86. S. V. Paramonov, Checking existence of solutions of par-tial differential equations in the fields of Laurent se-ries, Programming and Computer Software 40(2), 58-62(2014). Appendix A: List of Chaotic Systems
When testing our approach, we considered a variety of specificchaotic systems. Results for these are listed in Table 1. As theform and scaling of such equations can vary in the literature,we list the specific form of these equations used in our work.Let x, y, z, w ∈ C n ( R ) where n is the dimension of therelevant dynamical system, and let a, b, c, d ∈ R be parame-ters.Lorenz system [34,35]:˙ x = a ( y − x ) , ˙ y = x ( b − z ) − y , ˙ z = xy − cz . (46)Modified Chua’s circuit [ ? ,37,38]:˙ x = a (cid:18) y − (cid:0) x − x (cid:1)(cid:19) , ˙ y = x − y + z , ˙ z = − by . (47)Chen - Lee system [39,40]:˙ x = ax − yz , ˙ y = by + xz , ˙ z = cz + 13 xy . (48)Rabinovich - Fabrikant equations [41,42,43]:˙ x = y (cid:0) z − x (cid:1) + ax , ˙ y = x (cid:0) z + 1 − x (cid:1) + ay , ˙ z = − z ( b + xy ) . (49)R¨ossler system [44,45]:˙ x = − y − z , ˙ y = x + ay , ˙ z = b + z ( x − c ) . (50)6 Heather A. Harrington and Robert A. Van GorderChen system [46,47]:˙ x = a ( y − x ) , ˙ y = ( b − a ) x − xz + by , ˙ z = xy − cz . (51)L¨u system [48,49]:˙ x = a ( y − x ) , ˙ y = by − xz , ˙ z = xy − cz . (52)T system [50,51,52]:˙ x = a ( y − x ) , ˙ y = ( b − a ) x − axz , ˙ z = xy − cz . (53)4D Qi-Du-Chen-Chen-Yuan system [53]:˙ x = a ( y − x ) + yzw , ˙ y = b ( x + y ) − xzw , ˙ z = − cz + xyw , ˙ w = − dw + xyz . (54)Qi-Chen-Du-Chen-Yuan system [54]:˙ x = a ( y − x ) + yz , ˙ y = bx − y − xz , ˙ z = xy − cz . (55)Generalized Lorenz canonical form [55,56,57]:˙ x = ax − ( x − y ) z , ˙ y = − by − ( x − y ) z , ˙ z = − cz + ( x + y )( x + dy ) . (56)Two-parameter model for the blue-sky catastrophe [58,59,60]:˙ x = (cid:0) a − (cid:0) x + y (cid:1)(cid:1) x + y + 2 y + z , ˙ y = − z − (1 + y ) (cid:0) y + 2 y + z (cid:1) − x + ay , ˙ z = (1 + y ) z + x − b . (57)4D Lorenz - Stenflo system [61,62,63]:˙ x = a ( y − x ) + bw , ˙ y = cx − xz − y , ˙ z = xy − dz , ˙ w = − x − aw . (58)Genesio-Tesi system [64,65]:˙ x = y , ˙ y = z , ˙ z = ax + by + cz + x . (59)Arneodo-Coullet-Tresser system [66]:˙ x = y , ˙ y = z , ˙ z = ax − by − cz − x . (60) Appendix B: List of Hyperchaotic Systems
When testing our approach, we considered a variety of specifichyperchaotic systems. Results for these are listed in Table 2.As the form and scaling of such equations can vary in theliterature, we list the specific form of these equations used inour work.Let x, y, z, w ∈ C ( R ) and let a, b, c, d, e, f ∈ R be param-eters.4D R¨ossler flow [67]:˙ x = − y − z , ˙ y = x + 0 . y + w , ˙ z = 3 + xz , ˙ w = − . z + 0 . w . (61)Hyperchaotic Chen system [68,69]:˙ x = a ( y − x ) , ˙ y = − bx − xz + cy − w , ˙ z = xy − dz , ˙ w = x . (62)Hyperchaotic L¨u system [70]:˙ x = a ( y − x ) + w , ˙ y = by − xz , ˙ z = xy − cz , ˙ w = xz + dw . (63)Modified hyperchaotic L¨u system [71]:˙ x = a ( y − x + yz ) , ˙ y = by − xz + w , ˙ z = xy − cz , ˙ w = − dx . (64)Hyperchaotic Wang-Liu system [72]:˙ x = a ( y − x ) , ˙ y = bx − cxz + w , ˙ z = − dz + ex , ˙ w = − fx . (65)Hyperchaotic Jia system [73,74]:˙ x = a ( y − x ) + w , ˙ y = bx − xz − y , ˙ z = xy − cz , ˙ w = dw − xz . (66)Hyperchaotic Qi - van Wyk - van Wyk - Chen system [75,76]:˙ x = a ( y − x ) + yz , ˙ y = b ( x + y ) − xz , ˙ z = − cz − dw + xy , ˙ w = ez − fw + xy .xy .