- Email: [email protected]
Automatic phase space analysis of dynamical systems
Computing Systems in EngineeringVol. 1, Nos 2-4, pp. 607~614, 1990 Printed in Great Britain.
0956-0521/90$3.00+ 0.00 © 1990PergamonPress plc
AUTOMATIC PHASE SPACE ANALYSIS OF D Y N A M I C A L SYSTEMS E. SACKS Department of Computer Science, Princeton University, Princeton, NJ 08544, U.S.A.
(Recewed 27 March 1990) Abstract--Research on automating the analysis of nonlinear ordinary differential equations by combining artificial[intelligencetechniqueswith existingnumericalsoftware is discussed. A program called POINCARE that analyzes one-parameter autonomous planar systems at the level of experts through a combination of theoretical dynamics, numerical simulation, and geometrical reasoning is described. The input is the system, a bounding box for the system state, a bounding interval for the parameter, and error tolerances. POINCARE partitions the parameter interval into open subintervals of equivalent behavior bounded by bifurcation points, classifies the bifurcation points, and constructs representative phase diagrams for the subintervals. It detects local and global generic one-parameter bifurcations. It constructs the phase diagrams by identifying fixed points, saddle manifolds, and limit cycles and partitioning the remaining trajectories into open regions of uniform asymptotic behavior. It is explained how to extend POINCARE to two-parameter families of periodically driven planar equations. The basic algorithms carry over after some modifications and extensions. A survey of open issues in the automatic analysis of other dynamical systems is given.
l. INTRODUCTION This paper describes research on automating the analysis of nonlinear ordinary differential equations (ODEs) by combining artificial intelligence techniques with existing numerical software. Differential equations are a key modeling tool in virtually every field of science and engineering. They provide a rich language for modeling evolving physical systems, such as circuits, mechanisms, chemical reactions, plasma flow, weather patterns, and space vehicles. The solutions to the equations describe the time-varying behaviors of the physical systems. There are two main tasks in solving differential equations: (1) obtaining individual solutions for specific inputs and parameter values; and (2), combining the individual solutions into a global picture. Computers have revolutionized the first task, but have contributed little to the second. This is hardly surprising. The first task falls under the rubric of numerical analysis where computers made their earliest and greatest impact, whereas the second involves mathematical, geometric, and abstraction skills that computers are just beginning to attain. Dynamicists determine the global behavior of differential equations through a combination of theoretical and numerical analysis. Theory constrains the possible solutions, whereas numerics provides empirical answers to questions that theory cannot resolve. Theoretical knowledge guides the numerical analysis and helps interpret the output. Without this guidance, the fastest computer and the most powerful numerical software become slow and unreliable, since blindly searching the entire space of input conditions and parameter values requires a huge (in principal
infinite) amount of computation. Even were such a search possible, the voluminous output would be overwhelming. I aim to automate this human guidance, transforming computers from idiot savants to expert dynamicists. 2. PHASE SPACE
This section contains a brief introduction to the ideas and terminology used in the analysis of differential equations. Fuller treatments appear in Guckenheimer and Holmes ~ and in Hirsch and Smale. 2 I discuss autonomous planar systems .~ = f ( x , y )
y~=g(x, y),
(1)
where x and y are defined on open intervals of the real line. I return to the general case in Sec. 6. The phase space of the planar system (1) is the Cartesian product of the domains of x and y. Points in phase space represent states of the system. Curves on which the equations are satisfied, called trajectories, represent solutions. The topological properties of trajectories characterize the qualitative behavior of solutions. A point trajectory, called a fixed point, indicates an equilibrium solution, whereas a closed curve, called a limit cycle, indicates a periodic solution. A fixed point or limit cycle is an attractor when all nearby points approach it asymptotically and is a repellor when all nearby points move away (or, equivalently, approach it asymptotically in reverse time). The basin of an attractor or a repellor is the set of points that approach it or move away from it.
607
608
E. SACKS
>
(c)
(b)
Fig. 1. (a) Real attractor; (b) spiral repellor; and (c) saddle. The stable and unstable manifolds of the saddle appear as solid and broken thick lines.
h(x)
Fig. 2. Ball on a track and its phase diagram. Fixed points appear as solid circles labeled with the initials of their types. Saddle manifolds appear as thick lines. The Jacobian matrix
~f
ef
G
~y
(2)
~g ~g
determines the behavior of solutions near a fixed point. The fixed point is an attractor when the real parts of both eigenvalues are negative and a repellor when both are positive. Trajectories spiral toward the attractor or away from the repellor when the eigenvalues are complex, but move straight in or out when they are real. A saddle occurs when one eigenvalue is positive and the other is negative. Two trajectories,
called the stable manifold, approach the saddle asymptotically; two others, called the unstable manifold, move away; and all other nearby trajectories approach then move away. Figure 1 illustrates these behaviors. Dynamicists describe the trajectories of a system with a phase diagram, a partition of its phase space into attractors, basins, and unbounded trajectories, The phase diagram of a system concisely describes its qualitative behavior for all possible initial states: it shows where each trajectory ends up while abstracting away the details of transient behavior. F o r example, the system X=v = -- gv -- h ' ( x ) v
~x
!
(a)
(b)
Fig. 3. Alternate phase diagrams for the ball on the track. (a) Ix = 0; and (b) Ix < 0.
(3)
Automatic phase space analysis
Y
/ X
Fig. 4. Phase diagram of an aeroelastic galloping model. Thick lines indicate limit cycles. Thin lines indicate scanning trajectories used to detect limit cycles. models the motion of a ball on a track. The variables x and v represent location and velocity, h(x) is the height of the track, and la is the damping coefficient. Figure 2 shows the ball on the track and its phase diagram for h(x) = x4/4 - x2/2 and I~ = 1. The local minima and maxima o f h produce the attracting fixed points and the saddle, respectively. The stable manifold of the saddle delimits the basins of the attractors. Realistic models contain imprecisely specified parameters. A full description of a parameterized system, called a bifurcation diagram, consists of a partition of its parameter space into open regions of equivalent phase diagram,,; bounded by bifurcation curves where the phase diagram changes qualitatively. In typical bifurcations, fixed points or limit cycles appear, disappear, or change stability type. For example, varying ~t in the track equations yields two intervals of equivalent phase diagrams bounded by the bifurcation ~t = 0 where the attractors change to repellors (Fig. 3). 3. LIMITATIONS OF EXISTING SOFTWARE
Existing software packages assist dynamicists in the construction of phase and bifurcation diagrams.
609
Dynamicists find fixed points by Newton Raphson iteration or its cousins and calculate their eigenvalues and eigenvectors with linear algebra packages. They construct trajectories and sensitivity matrices with numerical integrators such as ODESSA. 3 They find fixed point bifurcations with continuation packages such as AUTO. 4 Several packages combine these capabilities with graphical output, for example INSITE, 5 KAOS, 6 M A C M A T H , 7 and P H A S E R ) The user must determine the inputs to the packages: initial guesses for fixed points, initial states for trajectories, and initial fixed points for continuation. The user must also interpret the numerical output and combine the results into a global picture. In the track example, existing packages can find the fixed points, determine their stability types, generate the saddle manifolds, and detect the bifurcation at I-t= 0. They cannot deduce that the resulting phase diagram and bifurcation diagram are complete. Nor can they construct complete diagrams for systems with limit cycles and other, more complicated attractors. 4. PLANAR SYSTEMS
As a first step toward automating the expertise of dynamicists, I tackled the simplest interesting class of equations: one-parameter families of autonomous planar equations. These systems are important in many applications and serve as a stepping stone to harder problems. I developed a program called P O I N C A R E that analyzes one-parameter families at the level of expert dynamicists. The input is the system, a bounding box of system states, an interval of parameter values, and error tolerances. P O I N C A R E partitions the parameter interval into subintervals of equivalent systems delimited by bifurcation points, picks a representative parameter value in each subinterval, and constructs the portion of its phase diagram in the bounding box. POINCARE constructs phase diagrams by identifying fixed points, saddle manifolds, and limit cycles in the bounding box and partitioning the remaining
bifurcation point •
sa
Fig. 5. Hopf bifurcation. bifurcation point
• sa
Fig. 6. L i m i t cycle bifurcation.
610
E. SACKS bifurcation point •
S • •
S
sr
Fig, 7. Saddle connection bifurcation. trajectories into basins and escape regions. These are the only possible asymptotic behaviors in a planar system by the Poincar6-Bendixson theorem. 2 It calculates the fixed points by setting all derivatives to zero and solving the resulting equations symbolically or numerically. It constructs saddle manifolds and limit cycles by numerical integration. The stable and unstable manifolds of a saddle are tangent to the eigenvectors of its negative and positive eigenvalues respectively. P O I N C A R E approximates the two trajectories in the stable manifold by integrating backward from two nearby points on the negative eigenvector, one on either side of the saddle. It approximates the unstable manifold by integrating forward from points on the positive eigenvector. Each limit cycle encloses a fixed point other than a saddle by an index theorem. P O I N C A R E scans out from the fixed points until it finds all enclosing limit cycles within the bounding box. F o r example, it finds three concentric limit cycles in the phase diagram shown in Fig. 4. P O I N C A R E generates this exact diagram as well as the other phase and bifurcation diagrams in the paper. Theoretical results provide P O I N C A R E with a framework for bifurcation analysis. One-parameter families of planar systems have a small set of possible bifurcations (except for rare nongeneric cases). Fixed points and limit cycles have three generic bifurcations. A fixed point undergoes a saddle node bifurcation at a parameter value where its Jacobian has a zero eigenvalue. As the parameter passes through that value, a saddle and a real attractor or repellor merge into a single degenerate fixed point then vanish. A fixed point undergoes a Hopf bifurcation when a conjugate pair of complex eigenvalues crosses the imaginary axis. Either a spiral attractor becomes a repellor and gives birth to an attracting limit cycle (Fig. 5), or a spiral repeUor becomes an attractor and gives birth to a repelling limit cycle. A limit cycle bifurcates at a parameter value where it attracts nearby trajectories from one side and repels them from the other. As the parameter passes through that value, a stable and unstable limit cycle collide and
co-destruct analogously to a saddle node bifurcation (Fig. 6). Fixed point and limit cycle bifurcations are called local bifurcations because they can be detected by tracking one phase space point through parameter space. Bifurcations that involve entire regions of phase space are called global. The most significant global bifurcations are ones where attractors and repellors appear and vanish. There are two generic possibilities in a planar system: a saddle connection, in which a saddle collides with and destroys a limit cycle (Fig. 7), and a saddle node infinite period (SNIPER), in which a saddle node bifurcation gives birth to a limit cycle (Fig. 8). POINCARE recognizes both types. It currently ignores secondary phenomena, such as changes in the orientation of the stable and unstable manifolds of saddles. POINCARE derives fixed point bifurcations directly from the input equations with a standard algebraic continuation package. 4 It detects the remaining bifurcations by means/ends analysis: it seeks bifurcations that explain the differences between the phase diagrams at every consecutive pair of "interesting" parameter values. It forms an initial list of interesting parameter values containing the endpoints of the parameter range and the midpoints between consecutive pairs of fixed point bifurcations. It compares the phase diagrams of each consecutive pair of parameter values on the list. If the diagrams are equivalent, it assumes that no further bifurcations occur at intermediate parameter values. If a single bifurcation explains the difference between the diagrams, it searches for such a bifurcation point. It looks up fixed point bifurcations and SNIPERs, which occur at fixed point bifurcations, in its precompiled list and finds other bifurcations by a binary search of the parameter interval. Otherwise, it assumes that multiple bifurcations occur between the parameter values and recursively searches the two subintervals formed by their midpoint. Comparing phase diagrams directly in means/ends analysis would be slow and complicated, since typical diagrams contain hundreds of points. Instead, the
bifurcation point
Fig. 8. SNIPER bifurcation.
Automatic phase space analysis m3: Phase diagram at Da=0.099596255: fixed-points: O. saddle at (0.40198523
611
1.9094298)
1. s p i r a l - s o u r c e at (0.8487833 4.031721) 2. s p i r a l - s i n k at (0.2238762 1.063412) limit cycles:
none
gl: Forward saddle-connection b i f u r c a t i o n at Da=0.10236387. nl: Phase diagram at Da:0.10257676: f i x e d - p o i n t s : O. saddle at (0.3358984 1.5955175) 1. s p i r a l - s o u r c e at (0.8578782 4.0749216) 2. s p i r a l - s i n k at (0.26993406 1.2821866) limit cycles:
stable enclosing
O, I, 2
Fig. 9. Symbolic descriptions produced by POINCARE.
Y
scriptors involves pattern-matching their common features and detecting discrepancies. 5. E X A M P L E
The continuously stirred tank reactor (CSTR) model from chemical engineering provides an example of complicated bifurcation patterns in a planar system. The equations are 2 = Da(1
j, = BDa(1
-~X
Fig. 10. Phase diagram for
D a = m 3.
comparison module employs compact symbolic descriptions of phase diagrams. It records the location and stability type of each fixed point, the stability type and enclosed fixed points of each limit cycle, and the forward and backward limit points of each saddle manifold. It finds the enclosed fixed points by standard computational geometry methods and obtains the other information directly from the phase diagram. Comparing two phase diagram de-
x)e' - x
--
-
x)e ) -
([3 +
l)y,
(4)
where x is the reactant concentration, y the reactor temperature, D a the Damk6hler number, B the heat of reaction, and [3 the heat transfer coefficient. Many research papers discuss the dynamics and bifurcations of this important system. POINCARE derives the complete bifurcation diagram for B = 19, 13= 3, parameter range 0.07 ~
/
Y
/ /
sr
,x Fig. 11. Phase diagram for
D a = n,
¢ after a saddle connection bifurcation.
612
E. SACKS
T
/
X
F(q)
Fig. 12. Trajectories of a T-periodic planar ODE and the corresponding return map F. bifurcation. Figures 10 and 11 show the corresponding phase diagrams. 6. LARGER SYSTEMS
P O I N C A R E demonstrates that a computer can match expert dynamicists in the analysis of autonomous planar ODEs, a problem that they solved (in principle) 40 years ago. 1° The next step is to automate the analysis of periodically driven planar equations, a problem at the frontier of dynamics research? ~''2 The Poincar6-Bendixson theorem breaks down, leading to a rich collection of additional asymptotic behaviors, including chaotic ones. The goal is to characterize trajectories and bifurcations, especially chaos and transition paths between order and chaos. Bifurcation analysis now involves a two-parameter space, since the behavior depends upon the frequency and amplitude of the drive function. The additional parameter introduces new generic bifurcations and complicates parameter-space search. We can reduce the analysis of a planar system with period T to that of a map on the state plane, called the return map. For each point p, there exists a unique trajectory whose value at time 0 is p. The return map carries p to the value of that trajectory at time T (Fig. 12). The return map F captures the essential dynamics of the periodic ODEs in a simpler setting. The set of iterates {F/(p)}, called the trajectory of p, describes the long-term behavior of the system given initial state p at time 0. Periodic trajectories of the ODEs with period T correspond to fixed points of the map where F ( p ) = p ; subharmonic periodic trajectories of order k correspond to periodic points where Fk(p) = p. Driven systems can also have quasiperiodic trajectories that contain two incommensurate periodic components and chaotic trajectories that wander erratically through phase space. These appear as closed invariant curves and fractal trajectories of the return map. The theory of planar maps closely resembles that of autonomous planar ODEs. The top-level analysis tasks in both cases are constructing phase diagrams and bifurcation diagrams. Although the maps have a
far richer set of behaviors than do the ODEs, the basic POINCARE analysis algorithm still applies. Finding and classifying fixed points works much as before. A fixed point is an attractor when the absolute values of both eigenvalues of the Jacobian are less than l, a repellor when both are greater than l, and a saddle when one is greater than and the other less than 1. Constructing saddle manifolds is harder because each manifold consists of a family of trajectories (Fig. 13), but still doable. 5 Applying these algorithms to F k yields the k periodic points and their manifolds. Another complication is that saddle manifolds can cross, yielding infinitely many periodic trajectories. In practice, we can safely ignore most of these trajectories because they attract tiny regions in phase space. We can find invariant curves by scanning outward from fixed and periodic points, much as POINCARE finds limit cycles. Once again, the invariant curve consists of a family of trajectories whose dynamics we must examine, whereas the limit cycle is a single trajectory with simple dynamics. Blind search is the only known way to find chaotic trajectories. POINCARE's basic bifurcation analysis strategy carries over to planar maps. Continuation of fixed and periodic points followed by means/ends analysis detects the generic bifurcations of one-parameter planar maps. These include all the ODE bifurcations
Fig. 13. Stable and unstable manifolds of a fixed point of a planar map.
Automatic phase space analysis as well as new bifurcations, such as a flip where an attracting period k point becomes repelling and gives birth to an attracting 2k point. The algorithm cannot find infinite sequences of bifurcations, which are common in maps. For example, an infinite sequence of flips collecting at a finite point produces chaos. We must modify the algorithm to postulate infinite sequences upon encountering a long, but finite initial segment. Adding a second parameter to the map changes the bifurcation analysis problem more dramatically. It increases the size of the search space quadratically and introduces new generic phenomena. The one-parameter bifurcations now occur along curves in parameter space, instead of at isolated points. New, more-complicated fixed point bifurcations occur at points in parameter space where the Jacobian has multiple eigenvalues with modulus one. New global phenomena also arise. We can still construct the curves of fixed point bifurcations by continuation, but finding the other bifurcation curves requires a planar search. I am currently extending the interval search algorithm of P O I N C A R E to the planar case. 7. RELATED WORK
A number of researchers have worked on automatic analysis of ordinary differential equations and maps. I developed a program that constructs piecewise linear approximations of autonomous planar differential equations and combines their local phase diagrams into a joint phase diagram. ~3 It chooses approximations that are accurate enough to reproduce the essential properties of their nonlinear prototypes, yet simple enough to be analyzed completely and efficiently. It derives the local phase diagrams from linear system theory and combines them by examining the boundaries between linear regions for possible transitions. This works well in systems with a small number of nonlinearities, but fails in essentially nonlinear systems, such as 2 = x 2. The program can sometimes deduce the existence of limit cycles, but can neither construct them explicitly nor detect their number or stability type. It also cannot construct attractor basins. It performs no bifurcation analysis. Yip's K A M program TM constructs phase diagrams for one-parameter area-preserving planar maps by numeric experimentation. A map is area-preserving if the area of any measurable set in phase space equals the area of its image. This models conservative phenomena, such as frictionless bouncing balls. KAM treats the equations as a black box for generating trajectories. It classifies trajectories according to their geometry and picks initial points for simulation according to the geometry of existing trajectories. It need not search for attractors and repellors, since these cannot exist in area-preserving maps. It finds local bifurcations by matching changes in the phase diagram against theoretical predictions. It does not find global bifurcations.
613
Abelson's bifurcation interpreter t5 finds periodic orbits and local bifurcations of one-parameter planar maps and one-parameter periodic planar ODEs. It calculates an initial periodic point, traces it in parameter space to find bifurcations, and interprets and summarizes the results according to bifurcation theory. It does not construct phase diagrams or detect global bifurcations. The AUTO software package perlbrms the same tasks, except for extrapolating infinite sequences of bifurcations, such as period-doubling. Unlike AUTO, the bifurcation interpreter searches a bounding box in phase space for all periodic orbits up to a user-specified maximum period. It also tests its conclusions for global consistency, for example merging an adjacent pair of unexplained bifurcations into a single generic bifurcation. Lee ~6classifies the response of a driven circuit over a range of driving frequencies and amplitudes. The circuit can contain any number of components of all types, including nonlinear ones, subject to the measurement constraints of her equipment. Her program works directly with the circuit, rather than with equations. It samples the frequencies and amplitudes at regular intervals, drives the circuit at each sample point, classifies the response as periodic, quasiperiodic, or chaotic, and groups the sample points into regions where the response is qualitatively identical. It classifies a response according to patterns in its time-signal and power-spectrum. For example, a single peak in the spectrum indicates a periodic response at the driving frequency. It chooses an arbitrary initial charge for the capacitors in the circuit, making no attempt to explore the phase space. Although unnecessary for planar ODEs, powerspectrum analysis is a powerful tool for classifying the behavior of planar maps and of larger systems of ODEs. 8. C O N C L U S I O N S
In this paper, I describe a program that analyzes one-parameter families of planar differential equations at the level of human experts. It constructs phase diagrams, performs bifurcation analysis, and summarizes its output into concise symbolic and graphical representations. I also sketch an extension to two-parameter periodic planar equations. What next? The primary dimensions for expansion are system type, phase space topology, phase space dimension, and parameter space dimension. Other system types include partial differential equations, integral equations, and delay equations. The topology of phase space affects the set of likely bifurcations and asymptotic behaviors. For example, planar ODEs on the torus have bounded trajectories that wrap around forever without approaching a fixed point or a limit cycle. Toridal and cylindrical phase spaces arise naturally in systems with periodic state variables, such as angles. We have already seen the effect of adding driving to planar ODEs. Increasing
614
E. SACKS
the dimension of phase space introduces additional phenomena, such as saddles with two-dimensional stable or unstable manifolds. Each new parameter gives rise to new bifurcation types and complicates parameter space search. Continuation techniques do not yet scale up to bifurcation surfaces. Essentially no work has been done on automating any of these extensions, but I expect that to change very soon.
7. 8. 9. 10.
Acknowledgements--The author is very grateful to Hal Abelson, Joel Friedman, John Guckenheimer, Yannis Kevrekidis, Gerry Sussman, Mike Wellman, and Ken Yip for useful advice on dynamics and suggestions regarding this work.
11.
12. REFERENCES
1. J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1986. 2. M. W. Hirsch and S. Smale, Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press College Division, Orlando, FL, 1974. 3. J. R. Leis and M. A. Kramer, "ODESSA--an ordinary differential equation solver with explicit simultaneous sensitivity analysis," ACM Transactions on Mathematical Software 14(1) (1988). 4. E. Doedel, "AUTO 86 user manual: software for continuation and bifurcation problems in ordinary differential equations," Technical Report, Princeton University, NJ, 1986. 5. T. S. Parker and L. O. Chua, Practical Numerical Algorithms for Chaotic Systems, Springer, New York, 1989. 6. J. Guckenheimer and S. Kim, "Kaos dynamical systems
13. 14.
15.
16.
toolkit with interactive graphics interface" (to be submitted). J. H. Hubbard and B. H. West, "Using the disk MacMath 6.1," Department of Mathematics, Cornell University, Ithaca, NY, 1986. H~ Kgcak, Differential and Difference Equations through Computer Experiments, Springer, New York, 1986. A. Uppal, H. Ray and A. B. Poore, "On the dynamic behavior of continuous stirred tank reactors," Chemical Engineering Science 29, 967-985 (1974). S. Lefschetz, Differential Equations: Geometric Theory, John Wiley, New York, 1962. D. G. Aronson, M. A. Chory, G. R. Hall and R. P. McGehee, "Bifurcations from an invariant circle for two-parameter families of maps of the plane: a computer-assisted study," Communications in Mathematical Physics 83, 303-354, 1982. B. B. Peckham, "The closing of resonance horns for periodically forced oscillators," Ph.D. thesis, University of Minnesota, MN, 1988. E. P. Sacks, "Automatic qualitative analysis of dynamic systems using piecewise linear approximations," Artificial Intelligence 41, 313-364 (1990). K. M. Yip, "Generating global behaviors using deep knowledge of local dynamics," Proceedings of the National Conference on Artificial Intelligence, American Association of Artificial Intelligence, 1988, pp. 280-285. H. Abelson, "The bifurcation interpreter: a step towards the automatic analysis of dynamical systems," AIM 1174, Massachusetts Institute of Technology, Artificial Intelligence Laboratory, 545 Technology Square, Cambridge, MA 02139, 1989 (to appear in Computers and Mathematics). M. K. Lee, "Summarizing qualitative behavior from measurements of nonlinear circuits," AI-TR 1125, Massachusetts Institute of Technology, Artificial Intelligence Laboratory, 545 Technology Square, Cambridge, MA 02139, 1989.
0956-0521/90$3.00+ 0.00 © 1990PergamonPress plc
AUTOMATIC PHASE SPACE ANALYSIS OF D Y N A M I C A L SYSTEMS E. SACKS Department of Computer Science, Princeton University, Princeton, NJ 08544, U.S.A.
(Recewed 27 March 1990) Abstract--Research on automating the analysis of nonlinear ordinary differential equations by combining artificial[intelligencetechniqueswith existingnumericalsoftware is discussed. A program called POINCARE that analyzes one-parameter autonomous planar systems at the level of experts through a combination of theoretical dynamics, numerical simulation, and geometrical reasoning is described. The input is the system, a bounding box for the system state, a bounding interval for the parameter, and error tolerances. POINCARE partitions the parameter interval into open subintervals of equivalent behavior bounded by bifurcation points, classifies the bifurcation points, and constructs representative phase diagrams for the subintervals. It detects local and global generic one-parameter bifurcations. It constructs the phase diagrams by identifying fixed points, saddle manifolds, and limit cycles and partitioning the remaining trajectories into open regions of uniform asymptotic behavior. It is explained how to extend POINCARE to two-parameter families of periodically driven planar equations. The basic algorithms carry over after some modifications and extensions. A survey of open issues in the automatic analysis of other dynamical systems is given.
l. INTRODUCTION This paper describes research on automating the analysis of nonlinear ordinary differential equations (ODEs) by combining artificial intelligence techniques with existing numerical software. Differential equations are a key modeling tool in virtually every field of science and engineering. They provide a rich language for modeling evolving physical systems, such as circuits, mechanisms, chemical reactions, plasma flow, weather patterns, and space vehicles. The solutions to the equations describe the time-varying behaviors of the physical systems. There are two main tasks in solving differential equations: (1) obtaining individual solutions for specific inputs and parameter values; and (2), combining the individual solutions into a global picture. Computers have revolutionized the first task, but have contributed little to the second. This is hardly surprising. The first task falls under the rubric of numerical analysis where computers made their earliest and greatest impact, whereas the second involves mathematical, geometric, and abstraction skills that computers are just beginning to attain. Dynamicists determine the global behavior of differential equations through a combination of theoretical and numerical analysis. Theory constrains the possible solutions, whereas numerics provides empirical answers to questions that theory cannot resolve. Theoretical knowledge guides the numerical analysis and helps interpret the output. Without this guidance, the fastest computer and the most powerful numerical software become slow and unreliable, since blindly searching the entire space of input conditions and parameter values requires a huge (in principal
infinite) amount of computation. Even were such a search possible, the voluminous output would be overwhelming. I aim to automate this human guidance, transforming computers from idiot savants to expert dynamicists. 2. PHASE SPACE
This section contains a brief introduction to the ideas and terminology used in the analysis of differential equations. Fuller treatments appear in Guckenheimer and Holmes ~ and in Hirsch and Smale. 2 I discuss autonomous planar systems .~ = f ( x , y )
y~=g(x, y),
(1)
where x and y are defined on open intervals of the real line. I return to the general case in Sec. 6. The phase space of the planar system (1) is the Cartesian product of the domains of x and y. Points in phase space represent states of the system. Curves on which the equations are satisfied, called trajectories, represent solutions. The topological properties of trajectories characterize the qualitative behavior of solutions. A point trajectory, called a fixed point, indicates an equilibrium solution, whereas a closed curve, called a limit cycle, indicates a periodic solution. A fixed point or limit cycle is an attractor when all nearby points approach it asymptotically and is a repellor when all nearby points move away (or, equivalently, approach it asymptotically in reverse time). The basin of an attractor or a repellor is the set of points that approach it or move away from it.
607
608
E. SACKS
>
(c)
(b)
Fig. 1. (a) Real attractor; (b) spiral repellor; and (c) saddle. The stable and unstable manifolds of the saddle appear as solid and broken thick lines.
h(x)
Fig. 2. Ball on a track and its phase diagram. Fixed points appear as solid circles labeled with the initials of their types. Saddle manifolds appear as thick lines. The Jacobian matrix
~f
ef
G
~y
(2)
~g ~g
determines the behavior of solutions near a fixed point. The fixed point is an attractor when the real parts of both eigenvalues are negative and a repellor when both are positive. Trajectories spiral toward the attractor or away from the repellor when the eigenvalues are complex, but move straight in or out when they are real. A saddle occurs when one eigenvalue is positive and the other is negative. Two trajectories,
called the stable manifold, approach the saddle asymptotically; two others, called the unstable manifold, move away; and all other nearby trajectories approach then move away. Figure 1 illustrates these behaviors. Dynamicists describe the trajectories of a system with a phase diagram, a partition of its phase space into attractors, basins, and unbounded trajectories, The phase diagram of a system concisely describes its qualitative behavior for all possible initial states: it shows where each trajectory ends up while abstracting away the details of transient behavior. F o r example, the system X=v = -- gv -- h ' ( x ) v
~x
!
(a)
(b)
Fig. 3. Alternate phase diagrams for the ball on the track. (a) Ix = 0; and (b) Ix < 0.
(3)
Automatic phase space analysis
Y
/ X
Fig. 4. Phase diagram of an aeroelastic galloping model. Thick lines indicate limit cycles. Thin lines indicate scanning trajectories used to detect limit cycles. models the motion of a ball on a track. The variables x and v represent location and velocity, h(x) is the height of the track, and la is the damping coefficient. Figure 2 shows the ball on the track and its phase diagram for h(x) = x4/4 - x2/2 and I~ = 1. The local minima and maxima o f h produce the attracting fixed points and the saddle, respectively. The stable manifold of the saddle delimits the basins of the attractors. Realistic models contain imprecisely specified parameters. A full description of a parameterized system, called a bifurcation diagram, consists of a partition of its parameter space into open regions of equivalent phase diagram,,; bounded by bifurcation curves where the phase diagram changes qualitatively. In typical bifurcations, fixed points or limit cycles appear, disappear, or change stability type. For example, varying ~t in the track equations yields two intervals of equivalent phase diagrams bounded by the bifurcation ~t = 0 where the attractors change to repellors (Fig. 3). 3. LIMITATIONS OF EXISTING SOFTWARE
Existing software packages assist dynamicists in the construction of phase and bifurcation diagrams.
609
Dynamicists find fixed points by Newton Raphson iteration or its cousins and calculate their eigenvalues and eigenvectors with linear algebra packages. They construct trajectories and sensitivity matrices with numerical integrators such as ODESSA. 3 They find fixed point bifurcations with continuation packages such as AUTO. 4 Several packages combine these capabilities with graphical output, for example INSITE, 5 KAOS, 6 M A C M A T H , 7 and P H A S E R ) The user must determine the inputs to the packages: initial guesses for fixed points, initial states for trajectories, and initial fixed points for continuation. The user must also interpret the numerical output and combine the results into a global picture. In the track example, existing packages can find the fixed points, determine their stability types, generate the saddle manifolds, and detect the bifurcation at I-t= 0. They cannot deduce that the resulting phase diagram and bifurcation diagram are complete. Nor can they construct complete diagrams for systems with limit cycles and other, more complicated attractors. 4. PLANAR SYSTEMS
As a first step toward automating the expertise of dynamicists, I tackled the simplest interesting class of equations: one-parameter families of autonomous planar equations. These systems are important in many applications and serve as a stepping stone to harder problems. I developed a program called P O I N C A R E that analyzes one-parameter families at the level of expert dynamicists. The input is the system, a bounding box of system states, an interval of parameter values, and error tolerances. P O I N C A R E partitions the parameter interval into subintervals of equivalent systems delimited by bifurcation points, picks a representative parameter value in each subinterval, and constructs the portion of its phase diagram in the bounding box. POINCARE constructs phase diagrams by identifying fixed points, saddle manifolds, and limit cycles in the bounding box and partitioning the remaining
bifurcation point •
sa
Fig. 5. Hopf bifurcation. bifurcation point
• sa
Fig. 6. L i m i t cycle bifurcation.
610
E. SACKS bifurcation point •
S • •
S
sr
Fig, 7. Saddle connection bifurcation. trajectories into basins and escape regions. These are the only possible asymptotic behaviors in a planar system by the Poincar6-Bendixson theorem. 2 It calculates the fixed points by setting all derivatives to zero and solving the resulting equations symbolically or numerically. It constructs saddle manifolds and limit cycles by numerical integration. The stable and unstable manifolds of a saddle are tangent to the eigenvectors of its negative and positive eigenvalues respectively. P O I N C A R E approximates the two trajectories in the stable manifold by integrating backward from two nearby points on the negative eigenvector, one on either side of the saddle. It approximates the unstable manifold by integrating forward from points on the positive eigenvector. Each limit cycle encloses a fixed point other than a saddle by an index theorem. P O I N C A R E scans out from the fixed points until it finds all enclosing limit cycles within the bounding box. F o r example, it finds three concentric limit cycles in the phase diagram shown in Fig. 4. P O I N C A R E generates this exact diagram as well as the other phase and bifurcation diagrams in the paper. Theoretical results provide P O I N C A R E with a framework for bifurcation analysis. One-parameter families of planar systems have a small set of possible bifurcations (except for rare nongeneric cases). Fixed points and limit cycles have three generic bifurcations. A fixed point undergoes a saddle node bifurcation at a parameter value where its Jacobian has a zero eigenvalue. As the parameter passes through that value, a saddle and a real attractor or repellor merge into a single degenerate fixed point then vanish. A fixed point undergoes a Hopf bifurcation when a conjugate pair of complex eigenvalues crosses the imaginary axis. Either a spiral attractor becomes a repellor and gives birth to an attracting limit cycle (Fig. 5), or a spiral repeUor becomes an attractor and gives birth to a repelling limit cycle. A limit cycle bifurcates at a parameter value where it attracts nearby trajectories from one side and repels them from the other. As the parameter passes through that value, a stable and unstable limit cycle collide and
co-destruct analogously to a saddle node bifurcation (Fig. 6). Fixed point and limit cycle bifurcations are called local bifurcations because they can be detected by tracking one phase space point through parameter space. Bifurcations that involve entire regions of phase space are called global. The most significant global bifurcations are ones where attractors and repellors appear and vanish. There are two generic possibilities in a planar system: a saddle connection, in which a saddle collides with and destroys a limit cycle (Fig. 7), and a saddle node infinite period (SNIPER), in which a saddle node bifurcation gives birth to a limit cycle (Fig. 8). POINCARE recognizes both types. It currently ignores secondary phenomena, such as changes in the orientation of the stable and unstable manifolds of saddles. POINCARE derives fixed point bifurcations directly from the input equations with a standard algebraic continuation package. 4 It detects the remaining bifurcations by means/ends analysis: it seeks bifurcations that explain the differences between the phase diagrams at every consecutive pair of "interesting" parameter values. It forms an initial list of interesting parameter values containing the endpoints of the parameter range and the midpoints between consecutive pairs of fixed point bifurcations. It compares the phase diagrams of each consecutive pair of parameter values on the list. If the diagrams are equivalent, it assumes that no further bifurcations occur at intermediate parameter values. If a single bifurcation explains the difference between the diagrams, it searches for such a bifurcation point. It looks up fixed point bifurcations and SNIPERs, which occur at fixed point bifurcations, in its precompiled list and finds other bifurcations by a binary search of the parameter interval. Otherwise, it assumes that multiple bifurcations occur between the parameter values and recursively searches the two subintervals formed by their midpoint. Comparing phase diagrams directly in means/ends analysis would be slow and complicated, since typical diagrams contain hundreds of points. Instead, the
bifurcation point
Fig. 8. SNIPER bifurcation.
Automatic phase space analysis m3: Phase diagram at Da=0.099596255: fixed-points: O. saddle at (0.40198523
611
1.9094298)
1. s p i r a l - s o u r c e at (0.8487833 4.031721) 2. s p i r a l - s i n k at (0.2238762 1.063412) limit cycles:
none
gl: Forward saddle-connection b i f u r c a t i o n at Da=0.10236387. nl: Phase diagram at Da:0.10257676: f i x e d - p o i n t s : O. saddle at (0.3358984 1.5955175) 1. s p i r a l - s o u r c e at (0.8578782 4.0749216) 2. s p i r a l - s i n k at (0.26993406 1.2821866) limit cycles:
stable enclosing
O, I, 2
Fig. 9. Symbolic descriptions produced by POINCARE.
Y
scriptors involves pattern-matching their common features and detecting discrepancies. 5. E X A M P L E
The continuously stirred tank reactor (CSTR) model from chemical engineering provides an example of complicated bifurcation patterns in a planar system. The equations are 2 = Da(1
j, = BDa(1
-~X
Fig. 10. Phase diagram for
D a = m 3.
comparison module employs compact symbolic descriptions of phase diagrams. It records the location and stability type of each fixed point, the stability type and enclosed fixed points of each limit cycle, and the forward and backward limit points of each saddle manifold. It finds the enclosed fixed points by standard computational geometry methods and obtains the other information directly from the phase diagram. Comparing two phase diagram de-
x)e' - x
--
-
x)e ) -
([3 +
l)y,
(4)
where x is the reactant concentration, y the reactor temperature, D a the Damk6hler number, B the heat of reaction, and [3 the heat transfer coefficient. Many research papers discuss the dynamics and bifurcations of this important system. POINCARE derives the complete bifurcation diagram for B = 19, 13= 3, parameter range 0.07 ~
/
Y
/ /
sr
,x Fig. 11. Phase diagram for
D a = n,
¢ after a saddle connection bifurcation.
612
E. SACKS
T
/
X
F(q)
Fig. 12. Trajectories of a T-periodic planar ODE and the corresponding return map F. bifurcation. Figures 10 and 11 show the corresponding phase diagrams. 6. LARGER SYSTEMS
P O I N C A R E demonstrates that a computer can match expert dynamicists in the analysis of autonomous planar ODEs, a problem that they solved (in principle) 40 years ago. 1° The next step is to automate the analysis of periodically driven planar equations, a problem at the frontier of dynamics research? ~''2 The Poincar6-Bendixson theorem breaks down, leading to a rich collection of additional asymptotic behaviors, including chaotic ones. The goal is to characterize trajectories and bifurcations, especially chaos and transition paths between order and chaos. Bifurcation analysis now involves a two-parameter space, since the behavior depends upon the frequency and amplitude of the drive function. The additional parameter introduces new generic bifurcations and complicates parameter-space search. We can reduce the analysis of a planar system with period T to that of a map on the state plane, called the return map. For each point p, there exists a unique trajectory whose value at time 0 is p. The return map carries p to the value of that trajectory at time T (Fig. 12). The return map F captures the essential dynamics of the periodic ODEs in a simpler setting. The set of iterates {F/(p)}, called the trajectory of p, describes the long-term behavior of the system given initial state p at time 0. Periodic trajectories of the ODEs with period T correspond to fixed points of the map where F ( p ) = p ; subharmonic periodic trajectories of order k correspond to periodic points where Fk(p) = p. Driven systems can also have quasiperiodic trajectories that contain two incommensurate periodic components and chaotic trajectories that wander erratically through phase space. These appear as closed invariant curves and fractal trajectories of the return map. The theory of planar maps closely resembles that of autonomous planar ODEs. The top-level analysis tasks in both cases are constructing phase diagrams and bifurcation diagrams. Although the maps have a
far richer set of behaviors than do the ODEs, the basic POINCARE analysis algorithm still applies. Finding and classifying fixed points works much as before. A fixed point is an attractor when the absolute values of both eigenvalues of the Jacobian are less than l, a repellor when both are greater than l, and a saddle when one is greater than and the other less than 1. Constructing saddle manifolds is harder because each manifold consists of a family of trajectories (Fig. 13), but still doable. 5 Applying these algorithms to F k yields the k periodic points and their manifolds. Another complication is that saddle manifolds can cross, yielding infinitely many periodic trajectories. In practice, we can safely ignore most of these trajectories because they attract tiny regions in phase space. We can find invariant curves by scanning outward from fixed and periodic points, much as POINCARE finds limit cycles. Once again, the invariant curve consists of a family of trajectories whose dynamics we must examine, whereas the limit cycle is a single trajectory with simple dynamics. Blind search is the only known way to find chaotic trajectories. POINCARE's basic bifurcation analysis strategy carries over to planar maps. Continuation of fixed and periodic points followed by means/ends analysis detects the generic bifurcations of one-parameter planar maps. These include all the ODE bifurcations
Fig. 13. Stable and unstable manifolds of a fixed point of a planar map.
Automatic phase space analysis as well as new bifurcations, such as a flip where an attracting period k point becomes repelling and gives birth to an attracting 2k point. The algorithm cannot find infinite sequences of bifurcations, which are common in maps. For example, an infinite sequence of flips collecting at a finite point produces chaos. We must modify the algorithm to postulate infinite sequences upon encountering a long, but finite initial segment. Adding a second parameter to the map changes the bifurcation analysis problem more dramatically. It increases the size of the search space quadratically and introduces new generic phenomena. The one-parameter bifurcations now occur along curves in parameter space, instead of at isolated points. New, more-complicated fixed point bifurcations occur at points in parameter space where the Jacobian has multiple eigenvalues with modulus one. New global phenomena also arise. We can still construct the curves of fixed point bifurcations by continuation, but finding the other bifurcation curves requires a planar search. I am currently extending the interval search algorithm of P O I N C A R E to the planar case. 7. RELATED WORK
A number of researchers have worked on automatic analysis of ordinary differential equations and maps. I developed a program that constructs piecewise linear approximations of autonomous planar differential equations and combines their local phase diagrams into a joint phase diagram. ~3 It chooses approximations that are accurate enough to reproduce the essential properties of their nonlinear prototypes, yet simple enough to be analyzed completely and efficiently. It derives the local phase diagrams from linear system theory and combines them by examining the boundaries between linear regions for possible transitions. This works well in systems with a small number of nonlinearities, but fails in essentially nonlinear systems, such as 2 = x 2. The program can sometimes deduce the existence of limit cycles, but can neither construct them explicitly nor detect their number or stability type. It also cannot construct attractor basins. It performs no bifurcation analysis. Yip's K A M program TM constructs phase diagrams for one-parameter area-preserving planar maps by numeric experimentation. A map is area-preserving if the area of any measurable set in phase space equals the area of its image. This models conservative phenomena, such as frictionless bouncing balls. KAM treats the equations as a black box for generating trajectories. It classifies trajectories according to their geometry and picks initial points for simulation according to the geometry of existing trajectories. It need not search for attractors and repellors, since these cannot exist in area-preserving maps. It finds local bifurcations by matching changes in the phase diagram against theoretical predictions. It does not find global bifurcations.
613
Abelson's bifurcation interpreter t5 finds periodic orbits and local bifurcations of one-parameter planar maps and one-parameter periodic planar ODEs. It calculates an initial periodic point, traces it in parameter space to find bifurcations, and interprets and summarizes the results according to bifurcation theory. It does not construct phase diagrams or detect global bifurcations. The AUTO software package perlbrms the same tasks, except for extrapolating infinite sequences of bifurcations, such as period-doubling. Unlike AUTO, the bifurcation interpreter searches a bounding box in phase space for all periodic orbits up to a user-specified maximum period. It also tests its conclusions for global consistency, for example merging an adjacent pair of unexplained bifurcations into a single generic bifurcation. Lee ~6classifies the response of a driven circuit over a range of driving frequencies and amplitudes. The circuit can contain any number of components of all types, including nonlinear ones, subject to the measurement constraints of her equipment. Her program works directly with the circuit, rather than with equations. It samples the frequencies and amplitudes at regular intervals, drives the circuit at each sample point, classifies the response as periodic, quasiperiodic, or chaotic, and groups the sample points into regions where the response is qualitatively identical. It classifies a response according to patterns in its time-signal and power-spectrum. For example, a single peak in the spectrum indicates a periodic response at the driving frequency. It chooses an arbitrary initial charge for the capacitors in the circuit, making no attempt to explore the phase space. Although unnecessary for planar ODEs, powerspectrum analysis is a powerful tool for classifying the behavior of planar maps and of larger systems of ODEs. 8. C O N C L U S I O N S
In this paper, I describe a program that analyzes one-parameter families of planar differential equations at the level of human experts. It constructs phase diagrams, performs bifurcation analysis, and summarizes its output into concise symbolic and graphical representations. I also sketch an extension to two-parameter periodic planar equations. What next? The primary dimensions for expansion are system type, phase space topology, phase space dimension, and parameter space dimension. Other system types include partial differential equations, integral equations, and delay equations. The topology of phase space affects the set of likely bifurcations and asymptotic behaviors. For example, planar ODEs on the torus have bounded trajectories that wrap around forever without approaching a fixed point or a limit cycle. Toridal and cylindrical phase spaces arise naturally in systems with periodic state variables, such as angles. We have already seen the effect of adding driving to planar ODEs. Increasing
614
E. SACKS
the dimension of phase space introduces additional phenomena, such as saddles with two-dimensional stable or unstable manifolds. Each new parameter gives rise to new bifurcation types and complicates parameter space search. Continuation techniques do not yet scale up to bifurcation surfaces. Essentially no work has been done on automating any of these extensions, but I expect that to change very soon.
7. 8. 9. 10.
Acknowledgements--The author is very grateful to Hal Abelson, Joel Friedman, John Guckenheimer, Yannis Kevrekidis, Gerry Sussman, Mike Wellman, and Ken Yip for useful advice on dynamics and suggestions regarding this work.
11.
12. REFERENCES
1. J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1986. 2. M. W. Hirsch and S. Smale, Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press College Division, Orlando, FL, 1974. 3. J. R. Leis and M. A. Kramer, "ODESSA--an ordinary differential equation solver with explicit simultaneous sensitivity analysis," ACM Transactions on Mathematical Software 14(1) (1988). 4. E. Doedel, "AUTO 86 user manual: software for continuation and bifurcation problems in ordinary differential equations," Technical Report, Princeton University, NJ, 1986. 5. T. S. Parker and L. O. Chua, Practical Numerical Algorithms for Chaotic Systems, Springer, New York, 1989. 6. J. Guckenheimer and S. Kim, "Kaos dynamical systems
13. 14.
15.
16.
toolkit with interactive graphics interface" (to be submitted). J. H. Hubbard and B. H. West, "Using the disk MacMath 6.1," Department of Mathematics, Cornell University, Ithaca, NY, 1986. H~ Kgcak, Differential and Difference Equations through Computer Experiments, Springer, New York, 1986. A. Uppal, H. Ray and A. B. Poore, "On the dynamic behavior of continuous stirred tank reactors," Chemical Engineering Science 29, 967-985 (1974). S. Lefschetz, Differential Equations: Geometric Theory, John Wiley, New York, 1962. D. G. Aronson, M. A. Chory, G. R. Hall and R. P. McGehee, "Bifurcations from an invariant circle for two-parameter families of maps of the plane: a computer-assisted study," Communications in Mathematical Physics 83, 303-354, 1982. B. B. Peckham, "The closing of resonance horns for periodically forced oscillators," Ph.D. thesis, University of Minnesota, MN, 1988. E. P. Sacks, "Automatic qualitative analysis of dynamic systems using piecewise linear approximations," Artificial Intelligence 41, 313-364 (1990). K. M. Yip, "Generating global behaviors using deep knowledge of local dynamics," Proceedings of the National Conference on Artificial Intelligence, American Association of Artificial Intelligence, 1988, pp. 280-285. H. Abelson, "The bifurcation interpreter: a step towards the automatic analysis of dynamical systems," AIM 1174, Massachusetts Institute of Technology, Artificial Intelligence Laboratory, 545 Technology Square, Cambridge, MA 02139, 1989 (to appear in Computers and Mathematics). M. K. Lee, "Summarizing qualitative behavior from measurements of nonlinear circuits," AI-TR 1125, Massachusetts Institute of Technology, Artificial Intelligence Laboratory, 545 Technology Square, Cambridge, MA 02139, 1989.