Visual exploration of 2D autonomous dynamical systems
VVisual exploration of 2D autonomous dynamical systems
Thomas M ¨uller
Visualisierungsinstitut der Universit¨at Stuttgart (VISUS),Allmandring 19, 70569 Stuttgart, GermanyE-mail:
Filip Sadlo
Visualisierungsinstitut der Universit¨at Stuttgart (VISUS),Allmandring 19, 70569 Stuttgart, GermanyE-mail:
Abstract.
In an introductory course on dynamical systems or Hamiltonian mechanics, vectorfield diagrams are a central tool to show a system’s qualitative behaviour in a certain domain.Because of their low sampling rates and the involved issues of vector normalization, theseplots give only a coarse insight and are unable to convey the vector field behaviour atlocations with high variation, in particular in the neighborhood of critical points. Similarly,automatic generation of phase portraits based on traditional sampling cannot precisely captureseparatrices or limit cylces. In this paper, we present
ASysViewer , an application for theinteractive visual exploration of two-dimensional autonomous dynamical systems using lineintegral convolution techniques for visualization, and grid-based techniques to extract criticalpoints and separatrices.
ASysViewer is addressed to undergraduate students during their firstcourse in dynamical systems or Hamiltonian mechanics.(Some figures may appear in colour only in the online journal)PACS numbers: 01.50.H-,01.50.hv
1. Introduction
The qualitative behaviour of two-dimensional autonomous dynamical systems is traditionallyexplored visually by means of either a vector field diagram or as a phase portrait that showscritical points, several typical orbits, and separatrices. Both techniques, however, have thedisadvantage that they only coarsely sample the region of interest. Furthermore, traditionalautomatic generation of these plots, in general, cannot precisely capture the separatrices orreproduce limit cylces. A continuous visualization of the phase space can be achieved usingline integral convolution (LIC), which is a well-known technique for vector field visualizationintroduced by Cabral and Leedom [1] in 1993. Unfortunately, up to now, this technique hasnot found its way into lectures about dynamical systems because of the lack of software thatcan produce high-quality phase images interactively.The aim of this article is to present how visualization and numerical techniques can helpin the visual and interactive exploration of two-dimensional autonomous dynamical systemswhich might be helpful in undergraduate courses on classical Hamiltonian mechanics aswell as on dynamical systems in general. For that, we mainly use line integral convolution a r X i v : . [ phy s i c s . e d - ph ] M a r isual exploration of 2D autonomous dynamical systems ASysViewer is based on the cross-platform framework Qt [5] and the OpenGraphics Library (OpenGL) [6]. The source code and Windows binaries are available from http://go.visus.uni-stuttgart.de/asysviewer .
2. Autonomous dynamical systems
A two-dimensional autonomous dynamical system is defined by two coupled first-orderdifferential equations˙ x = P ( x , y ) , ˙ y = Q ( x , y ) , (1)where a dot represents differentiation with respect to time, t , and P and Q are two (arbitrary) C -differentiable functions of x and y in an open subset Ω ⊂ R . Both, x and y themselvesare functions of t . Because of the theorem of existence and uniqueness, solutions of Eq. (1)do not cross. Hence, we can interpret the system of equations (1) geometrically as a two-dimensional vector field with a vector (cid:126) f ( x , y ) = ( P ( x , y ) , Q ( x , y )) T attached at each point ( x , y ) ∈ Ω . The trace of a solution (cid:126) σ ( t , x , y ) = ( σ x ( t , x , y ) , σ y ( t , x , y )) T of the coupledsystem (1) with initial values ( x , y ) and parameter t is called integral curve (trajectory, orbit,flow line, characteristic line, or streamline), where (cid:126) σ is always tangential to the vector (cid:126) f atthe current point ( x , y ) . Contour lines for which d y / d x = ˙ y / ˙ x = const are called isoclines .To obtain a qualitative view of the vector field, a phase portrait that consists of criticalpoints and several characteristic lines and separatrices is used. Stationary points (alsocalled fixed or equilibrium points) are defined by ˙ x = ˙ y = Critical points are isolatedstationary points, i.e., stationary points where the determinant of the Jacobian does not vanish,det ( J ) (cid:54) =
0. The behaviour of the system (1) in a close neighborhood of a critical point ( x c , y c ) is characterized by the eigenvalues λ , of the Jacobian J = (cid:18) J J J J (cid:19) = (cid:18) ∂ x P ∂ y P ∂ x Q ∂ y Q (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) x c , y c , (2)which follows from linearizing the system (1) with a Taylor series expansion,˙ X = X ∂ P ∂ x (cid:12)(cid:12)(cid:12)(cid:12) x c , y c + Y ∂ P ∂ y (cid:12)(cid:12)(cid:12)(cid:12) x c , y c , ˙ Y = X ∂ Q ∂ x (cid:12)(cid:12)(cid:12)(cid:12) x c , y c + Y ∂ Q ∂ y (cid:12)(cid:12)(cid:12)(cid:12) x c , y c , (3)and X = x − x c , Y = y − y c . Then, the eigenvalues λ , can be determined by means of thecharacteristic equation det ( J − λ ) = λ , = (cid:20) tr ( J ) ± (cid:113) ( tr ( J )) − ( J ) (cid:21) (4)with trace, tr ( J ) = J + J , determinant det ( J ) = J J − J J , and being the identitymatrix in two dimensions. If the real part of the eigenvalues is nonzero, the critical point isual exploration of 2D autonomous dynamical systems ( J ) and det ( J ) . If det ( J ) <
0, the critical point canonly be a saddle point. If det ( J ) >
0, it depends on tr ( J ) whether the critical point is stable orunstable, i.e., in which way orbits approach or move away from the critical point, respectively. tr( J )det( J )Saddle point Unstable nodeStable node Unstable focusStable focus C e n t e r U S N o r U D N SS N o r S D N Figure 1.
Classification of critical points depending on the trace and the determinant of theJacobian J (see also Lynch[2]). The parabola is given by det ( J ) = ( tr ( J )) /
4, see Eq. (4).
Figure 2 shows examples of the six main types of critical points which are located at thecenters of the vector field diagrams. The arrows are normalized to unit length to better see thecharacteristics of the field. The red curves represent some typical orbits.To further characterize a vector field, one has to determine the characteristic lines thatmark the boundary between different regions. These lines are also called separatrices, ingeneral. A characteristic line is called homoclinic orbit if the limitslim t → ∞ (cid:126) σ ( t ,(cid:126) x ) = lim t →− ∞ (cid:126) σ ( t ,(cid:126) x ) = (cid:126) x c (5)have the same critical point (cid:126) x c . If the limits have different critical points, the characteristicline is called heteroclinic orbit ; otherwise, it is a generic orbit .
3. Visualization techniques
The most straightforward visualization of the two-dimensional system (1) is a vector fielddiagram with arrows pointing in the direction (cid:126) f = ( P , Q ) T as already shown in figure 2.Figures 3(a),(b) show a vector field diagram where the arrows are either normalized or not.Without normalization (Fig. 3(a)), the arrow lengths are too diverse so that the overall structureof the vector field gets lost if the diagram cannot be scaled interactively.While in figure 3(b) the critical points at (cid:126) c = ( − , ) T and (cid:126) c = ( , ) T are obviously acenter and a saddle, respectively, the location and the types of the critical points in figure 3(c), (cid:126) c = ( , ) T , (cid:126) c = ( . , − . ) T , (cid:126) c = ( . , − . ) T , are not as obvious. Even morecomplicate is the retrieval of the separatrices in figure 3(c). Line integral convolution is a well-known texture synthesis technique in flow visualization. Ithas the advantage that its spatial resolution is much higher than that of a vector field diagram. isual exploration of 2D autonomous dynamical systems -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (a) unstable node -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (b) stable node -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (c) saddle -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (d) unstable focus -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (e) center -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (f) stable focus Figure 2.
Vector plots for critical points at ( x c = , y c = ) . The arrows are normalized to unitlength. The red lines show some example orbits. -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (a) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2-2-1.5-1-0.500.511.52 xy (b) -1 -0.5 0 0.5 1 1.5 2 2.5 3-2-1.5-1-0.500.511.52 xy (c) Figure 3.
Vector field diagram for P ( x , y ) = y , Q ( x , y ) = x + x with critical points (cid:126) c =( − , ) T and (cid:126) c = ( , ) T ; (a) real length, (b) normalized. (c) Normalized vector fielddiagram for P ( x , y ) = x − x − y , Q ( x , y ) = x ( + xy ) with critical points (cid:126) c = ( , ) T , (cid:126) c = ( . , − . ) T , and (cid:126) c = ( . , − . ) T . The idea is to locally smear an image (usually an image consisting of white noise) along thedirection of a vector field. Thereby, pixel values along streamlines are strongly correlatedwhile orthogonal to the vector field, there is virtually no correlation.The basic algorithm works as follows. First, we have to generate a window-filling noise isual exploration of 2D autonomous dynamical systems T of size N x × N y pixels with random numbers in the range [ , ] for each pixel thatcovers the domain Ω , see figure 4. x min x max y min y max Ω (a) (b) Figure 4. (a) 2 D domain of interest Ω = { ( x min , y min ) , ( x max , y max ) } . (b) Noise texture/image T with resolution N x × N y pixels, slightly smoothed by a Gaussian kernel to reduce highfrequencies. Next, we perform a discrete line integral convolution of the noise image for each pixel p ∈ T using the trajectory (cid:126) σ ( t ,(cid:126) x ) . The resulting intensity I ( p ) at the pixel p follows from I ( p ) = k ( ) T ( p ) + N ∑ n = [ k ( n ) T ( u n ) + k ( − n ) T ( v n )] , (6)where k is a filter kernel and 2 N + T ( p ) ∈ [ , ] of the noise texture at pixel p weighted by k ( ) and then sum up all the intensityvalues of T along the positive and negative parts of the trajectory weighted by the kernel value k at the corresponding sampling points u n and v n , respectively.” For that, we first need therelation between pixel coordinates q = ( q x , q y ) ∈ T and domain coordinates (cid:126) x = ( x , y ) T ∈ Ω which is defined by (cid:126) x (cid:55)→ q with q x = N x x − x min x max − x min , q y = N y y − y min y max − y min . (7)For the sake of convenience, we employ the Euler method to integrate the trajectories (cid:126) σ ( t ,(cid:126) u = (cid:126) v ) in positive and negative direction. The initial pixel p defines the initial positions (cid:126) u = (cid:126) v which follow from the inverse of equation (7). Then, Euler iteration delivers thesampling points (cid:126) u n ,(cid:126) v n ∈ Ω with (cid:126) u n + = (cid:126) u n + h (cid:126) f n , (cid:126) v n + = (cid:126) v n − h (cid:126) f n , n = , . . . , N − . (8)According to equation (7), we obtain the sampling points u n , v n ∈ T from (cid:126) u n and (cid:126) v n .The filter kernel for standard LIC can be any symmetric function. Here, we use a Hannfilter k ( n ) = N (cid:16) + cos π nN (cid:17) , n = {− N , . . . , N } , N ∑ n = − N k ( n ) = . (9)The line integral convolution is computationally quite expensive. But with the highparallelism of GPUs, we can accomplish the convolution at interactive rates. At the end,we obtain the LIC image composed of all the intensity values I ( p ) calculated above, whichcan afterwards be scaled arbitrarily to enhance the field structures, see, e.g., figure 6.Unfortunately, the standard LIC loses the directional information of the vector field, i.e.,forward and reverse direction are indiscernible. To resolve this disadvantage, we can use an isual exploration of 2D autonomous dynamical systems oriented LIC ), see, e.g., Wegenkittl etal. [7]. The most straightforward kernel for that is a sawtooth. Furthermore, for oriented LIC,we have to replace the noise texture of figure 4 by a texture with a finite number of randomlydistributed spots and black background. Figure 5 shows the vector fields of figure 3 usingoriented LIC. Additionally, we use a repeating colour map for the vector magnitudes (cid:107) (cid:126) f (cid:107) . Figure 5.
Oriented line integral convolution images of the vector fields of figure 3, generatedby
ASysViewer using the scripts paperExp1.js and paperExp2.js . Here, colour encodesvector magnitude (cid:107) (cid:126) f (cid:107) . As these values are very diverse, we repeat the colour map with n ∈ N . Another possibility to maintain the direction information of the vector field is to animatethe line integral convolution ( animated LIC ). For that, we change the frequency of the Hannfilter and modulate the phase by animation time τ . Additionally, we multiply the modifiedHann filter with a Hann window. Thus, k ( n ) from equation (9) is replaced by h ( n ) = (cid:20) + cos (cid:18) π nN + τ (cid:19)(cid:21) (cid:16) + . π nN (cid:17) . (10) So far, we have demonstrated the benefits of dense vector field visualizations by the LICfamily of algorithms as compared to the discrete sparse visualizations by vector plots: dueto their dense nature, they provide detailed information about the direction of (cid:126) f , and to someextent also about its magnitude and orientation. However, it is still difficult to spot the criticalpoints and, in particular, the separatrices emanating therefrom (figure 5). In other words, it isstill difficult to identify regions of qualitatively different behaviour.Topology extraction usually starts with the extraction of critical points. A major difficultywith this step is, however, the lack of a theory on the existence and location of critical pointsin arbitrary (higher-order) vector fields. Nevertheless, in the special class of low-order vectorfields arising from 2D (bilinear) tensor-product linear interpolation, their number is restrictedto two. Hence, in vector fields where Ω is discretized into cells and (cid:126) f is discretized at the nodes(the corners of the cells), such a field is present in each cell due to bilinear interpolation.The algorithm to detect the critical points works as follows:(i) Detect cells with different signs at their nodes both in terms of P and Q . isual exploration of 2D autonomous dynamical systems P and Q intersect.(iii) Analytical derivation of these intersections is possible, but subject to numericalinstability. Thus, we follow a subdivision approach:(iv) Each cell with different signs in P and Q is subdivided, and the resulting four childrencells are again tested for different signs in P and Q . This procedure is repeated until nocell exhibits different signs in P and Q or until a maximum subdivision level (a prescribedaccuracy) is achieved. The centers of the remaining cells with different signs representstationary points.(v) Those stationary points that have det ( J ) (cid:54) = ( J ) <
4. The viewer
A screenshot of
ASysViewer ’s graphical user interface is given in figure 6. The main windowshows standard LIC with critical points and separatrices of the “mathematical pendulum withfriction” example, see Eq. (11). In the “
System ” window, the functions P ( x , y ) and Q ( x , y ) can be given as mathematical expressions, and the domain Ω can be controlled either usingmouse interaction (zooming: middle button; panning: left button) or by typing the valuesin the corresponding boxes. With the right mouse button a trajectory can be shown whichstarts at the current mouse position. The parameters for this trajectory can be set within the“ Trajectory ” window. The “
FuncPlot ” window shows the values of P and Q along thewhite intersection line drawn in the main window. This line can be set with the right mousebutton and “ intersection ” selected in the “ Show ” window. After defining the dynamicalsystem, the critical points as well as the separatrices are determined by pressing the “ calc ”button in the “
Topology ” window.There are also two windows to control the parameters a and eps which can be used as freeparameters to define the functions P and Q . is based on the cross-platform application and user interface framework Qt [5],the graphics library OpenGL, and the OpenGL Shading Language (GLSL) [6]. In the mainwindow, we draw a single quad that covers the whole 2 D domain of interest, see also figure 4.For each fragment (pixel) of this quad, the line integral convolution is calculated within afragment program/shader that works as a SIMD (single instruction multiple data) architecture.The fragment program is written in a C-like programming language (GLSL) and is compiledand transferred to the GPU at runtime. When the functions P and Q are modified, the functionstrings are inserted in the fragment code and the fragment program is recompiled. Dueto limitations of mathematical functions/operators in GLSL, expressions like “ x ” must beformulated as “ x ∗ x ∗ x ∗ x ”. isual exploration of 2D autonomous dynamical systems Figure 6.
Screenshot of
ASysViewer ’s graphical user interface showing standard LIC withseparatrices and critical points. In the “
FuncPlot ” window, the blue and red lines correspondto the functions P and Q evaluated along the interactively defined white line in the mainwindow. The integration of the trajectories to visualize the separatrices or individual orbitsselected via mouse control is done using a step-controlled Runge-Kutta-Fehlberg integrator.For that, the functions P and Q are parsed using the muparser [8] library by Ingo Berg. Forbetter reproducibility, the parameter settings of ASysViewer are also scriptable by means ofQtScript which is based on the ECMAScript [9] standard. The necessary scripting functionscan be checked up from the source code.
5. Examples
As a first example, we consider the damped oscillation of the mathematical pendulum of unitlength which is described by the equation of motion ¨ ϕ + ω sin ϕ + α ˙ ϕ = ϕ , angular frequency ω , and a positive frictional coefficient α ≥
0, see figure 6.Substituting ϕ = x , ˙ ϕ = y converts the second order ordinary differential equation into thestandard form of equation (1),˙ x = P ( x , y ) = y , ˙ y = Q ( x , y ) = − ω sin x − α y . (11)The positions of the corresponding critical points immediately follow from ˙ x = = ˙ y . Thus, y c = x c = sin ( n π ) with n ∈ Z . The Jacobian yields tr ( J ) = α and det ( J ) = ω cos ( n π ) .If n is odd, the determinant is negative and the critical point is a saddle. If n is even, on theother hand, the critical point is either a center (undamped oscillator, α = α < ω cos ( n π ) ), or an unstable node. In the damped case, each unstable focus/nodeis connected to the two neighbouring saddle points via heteroclinic orbits, whereas in theundamped case, each saddle point is connected to its two neighbouring saddle points. isual exploration of 2D autonomous dynamical systems A limit cylce is an isolated periodic solution which is defined by the set of all points (cid:126) y with (cid:126) σ ( t ,(cid:126) x ) → (cid:126) y for t → ∞ . An example can be found in the phase space for the differential equationrelated to the oscillation of a violin string as derived by Rayleigh [10],¨ v − µ (cid:18) −
13 ˙ v (cid:19) ˙ v + v = , µ > , (12)which can be brought into the standard form (Eq. (1)) via v = x and ˙ v = y . (a) -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5-2.5-2-1.5-1-0.500.511.522.5 xy (b) Figure 7. (a) Phase space image of the Rayleigh system, Eq. (12), visualized using orientedLIC. The trajectory (yellow line) starting at x = . , y = . Ω = [ − . , . ] × [ − . , . ] and µ =
1. (b) Vector field diagram.
As shown in the phase space image, figure 7(a), the oriented LIC representation of theRayleigh system indicates that there must exist a boundary curve separating the outflowingand inflowing vector field. Every integral curve (cid:126) σ ( t ,(cid:126) x i ) approaches this boundary curve (limitcycle) irrespective of the initial point (cid:126) x i . The
Poincar´e index i P of a closed curve Γ is determined by integrating the change in the angleof the vectors at each point of Γ along Γ . For the numerical calculation, we set Γ to a circleof radius r and replace the integration by the sum i p = π N ∑ n = arcsin (cid:107) (cid:126) f n × (cid:126) f n + (cid:107)(cid:107) (cid:126) f n (cid:107) · (cid:107) (cid:126) f n + (cid:107) , (13)where (cid:126) f n = (cid:126) f ( r cos ϕ n , r sin ϕ n ) T and ϕ n = π n / N . For a sink, source, or center, we have i P = +
1. A saddle point has i P = −
1. For details to the Poincar´e index, we refer the readerto Wiggins [3]. In
ASysViewer , the Poincar´e index of an arbitrary point (cid:126) x can be determinedinteractively by constructing a circle around (cid:126) x via mouse handling (right button). For that,select “ poincare idx ” in the mouse combo box of the “ Show ” window. isual exploration of 2D autonomous dynamical systems Consider the system˙ x = ax − x , ˙ y = − y , (14)with free parameter a ∈ R . If a is continuously varied from negative to positive values, thissystem undergoes a pitchfork bifurcation . As long as a is negative, there is only a singlecritical point at (cid:126) x neg = ( , ) T . For a =
0, there is one non-hyperbolic critical point at theorigin. If a >
0, we obtain three critical points (cid:126) x = ( , ) T (saddle), (cid:126) x = ( −√ a , ) T (stablenode), and (cid:126) x = ( √ a , ) T (stable node).Figure 8 shows the “pos,neg P , Q ” image for the pitchfork bifurcation system with a = P and Q . Hence, a point where all fourcolours hit is a stationary point, P = Q = (a) “pos,neg P , Q ” image for a = ax x = √ ax = −√ a (b) Bifurcation diagram Figure 8.
The system (Eq. (14)) shows a pitchfork bifurcation at a =
0. If a >
0, the criticalpoint at the origin is a saddle whereas the other two are stable nodes. The colour coding in(a) is as follows: yellow ( P > , Q > P > , Q < P < , Q < P < , Q >
6. Summary
In this work we presented a framework for interactive visual exploration of 2 D autonomousdynamical systems. For that, we use different types of line integral convolution techniques fora dense representation of the vector field. We also described a robust algorithm to numericallydetect critical points. As the source code of ASysViewer is freely available, the interestedstudent is encouraged to conduct his or her own experiments and to look also at the codingdetails.
Acknowledgments
This work was partially funded by Deutsche Forschungsgemeinschaft (DFG) as part ofthe Collaborative Research Centre SFB 716 and the Cluster of Excellence in Simulation isual exploration of 2D autonomous dynamical systems
References [1] B. Cabral and L. C. Leedom. Imaging vector fields using line integral convolution. In
Proc. 20th ann. conf.Comp. graph. and inter. techn. , SIGGRAPH ’93, pages 263–270, New York, USA, 1993.[2] S. Lynch.
Dynamical Systems with Applications using Mathematica . Birkh¨auser, Boston, 2007.[3] S. Wiggins.
Introduction to Applied Nonlinear Dynamical Systems and Chaos . Springer, 2003.[4] R. S. Laramee, H. Hauser, H. Doleisch andB. Vrolijk, F. H. Post, and D. Weiskopf. The state of the art in flowvisualization: Dense and texture-based techniques.
Computer Graphics Forum , 23:2004, 2004.[5] Qt – A cross-platform application and UI framework. An LGPL version is freely available at qt-project.org .[6] Details about OpenGL can be found at .[7] R. Wegenkittl, E. Gr¨oller, and W. Purgathofer. Animating Flowfields: Rendering of Oriented Line IntegralConvolution. In
Computer Animation ’97 Proceedings , pages 15–21, 1996.[8] “muparser – fast math parser library” can be found at http://muparser.beltoforion.de . Here, we use theversion 2.2.4 from 2014-12-10.[9] ECMA-262 language specification: .[10] J. Rayleigh.