Computational Method for Phase Space Transport with Applications to Lobe Dynamics and Rate of Escape

Lobe dynamics and escape from a potential well are general frameworks introduced to study phase space transport in chaotic dynamical systems. While the former approach studies how regions of phase space are transported by reducing the flow to a two-dimensional map, the latter approach studies the phase space structures that lead to critical events by crossing periodic orbit around saddles. Both of these frameworks require computation with curves represented by millions of points-computing intersection points between these curves and area bounded by the segments of these curves-for quantifying the transport and escape rate. We present a theory for computing these intersection points and the area bounded between the segments of these curves based on a classification of the intersection points using equivalence class. We also present an alternate theory for curves with nontransverse intersections and a method to increase the density of points on the curves for locating the intersection points accurately.The numerical implementation of the theory presented herein is available as an open source software called Lober. We used this package to demonstrate the application of the theory to lobe dynamics that arises in fluid mechanics, and rate of escape from a potential well that arises in ship dynamics.



There are no comments yet.


page 3

page 7

page 13

page 19

page 22

page 24

page 25

page 27


A Crossing Lemma for Jordan Curves

If two Jordan curves in the plane have precisely one point in common, an...

A Crossing Lemma for Families of Jordan Curves with a Bounded Intersection Number

A family of closed simple (i.e., Jordan) curves is m-intersecting if any...

Bezier Curves Intersection Using Relief Perspective

Presented paper describes the method for finding the intersection of cla...

Computing the intersection of two quadrics through projection and lifting

This paper is devoted to presenting a method to determine the intersecti...

Elementary exposition of realizing phase space structures relevant to chemical reaction dynamics

In this article, we review the analytical and numerical approaches for c...

Discriminants of complete intersection space curves

In this paper, we develop a new approach to the discrimi-nant of a compl...

Enumerating Hassett's wall and chamber decomposition of the moduli space of weighted stable curves

Hassett constructed a class of modular compactifications of the moduli s...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In chaotic dynamical systems, phase space transport is an approach for understanding and characterizing how regions of phase space move over time. This approach is used widely in studying escape and transition rate in classical mechanics and dynamical astronomy [16, 10, 44, 45], loss of global integrity in ship capsize [24], reaction and dissociation rate in chemical reactions [23, 13, 41], chaotic advection in fluid mechanics [2, 1], wake generation behind a cylinder in a fluid flow [12], and transport in geophysical flow [36]. Furthermore, phase space transport is also relevant for prediction, control, and design in a myriad of natural processes and engineering systems; see Ref. [35] for an overview. In this approach, lobe dynamics, introduced in Ref. [32]

, is a geometric framework for studying the global transport in 1 and 2 degree-of-freedom (DOF) Hamiltonian systems that can be reduced to two-dimensional maps. Moreover, there has been few attempts at extending lobe dynamics to higher dimensional Hamiltonian systems with small perturbation where the existence of normally hyperbolic invariant set and its associated codimension-1 stable and unstable manifolds can be guaranteed; see Refs. 

[42, 13, 8, 9, 19] for details. In this article, we present a theory and numerical methods which are motivated by lobe dynamics in 2-DOF Hamiltonian systems that can be reduced to two-dimensional map by using suitable Poincaré section, or return map, or time-T map. Hence, the methods developed herein are applicable to a diverse array of problems in physical and engineering sciences. Following the developments in [32, 34], lobe dynamics can be stated as a systematic study of a fates and histories of a set of initial conditions. Formally speaking, the two-dimensional phase space of a Poincaré map (in general, this can a return map or a time-T map) can be partitioned into regions with boundaries consisting of parts of the boundary of (which may be at infinity) and/or segments of stable and unstable invariant manifolds of hyperbolic fixed points, , as shown schematically in Fig. 1(a). When the manifolds and are followed out on a global scale, they often intersect in primary intersection points or pips, for example as in Fig. 1(b), and secondary intersection points or sips, for example in Fig. 1(b). More precisely, an intersection point is called a pip if the curves connecting the point and the hyperbolic points intersect only at that point, and is called a sip otherwise. These intersections allow one to define boundaries between regions , as illustrated in Fig. 1(b). Moreover, the transport between regions of phase space can be completely described by the dynamical evolution of parcels of phase space enclosed by segments of the stable and unstable manifolds called lobes.

(a) (b)
Figure 1: (a) Pieces of the local unstable and stable manifolds, (red) and (green) of saddle fixed points . (b) When the manifolds and are followed out on a global scale, they often intersect in primary intersection points . These intersections allow one to define boundaries between regions .

Classically, invariant manifolds of hyperbolic fixed points are computed for as long as possible and lobes are extracted from the two curves [32, 10]. However, it is possible to get more lobes, hence compute transport for a much longer period of time, by integrating the boundaries directly, particularly in the complex case of multiple, self-intersecting, lobes [10]. In this approach, the notion of pips and sips are not suitable, and thus we propose a generalization of definition of lobes for two intersecting closed curves. When applying lobe dynamics to transport problems, it is of eventual interest to quantify volume of phase space that crosses the boundary with higher iterates of the Poincarè map . This is typically achieved in multiple ways:

  • Distribute test points inside the phase space region of interest and computing the iterate/s required for each to escape. This will also need to take care of the re-entrainment due to the underlying turnstile (a pair of lobe that moves points across a boundary) mechanism [20, 32, 43, 33].

  • Constructing a functional of the nonlinear system of vector field which measures area between the manifolds as parametrized by time. This is a semi-analytical approach and only remains valid near small perturbations, for example Melnikov method for mixing, stirring, optimal phase space flux, and action-integral method 

    [25, 37, 4, 5, 29, 6].

  • Following the boundaries of separatrices/manifolds as it is evolved in time and compute set operations with its pre-images/images. This is a more general but surely a difficult approach due to the stretching and folding of the curves that are involved in such computations [26, 27, 28].

For concreteness, let us consider the lobes formed due to a heteroclinic tangle as shown in Fig. 1(b) or due to a homoclinic tangle as shown in Fig. 2. In both cases, the lobes are regions of phase space bounded by curves oriented in opposite directions, this is because the segments bounding the lobes are stable and unstable manifold of the hyperbolic fixed points, and (Fig. 1(b)), or hyperbolic fixed point (Fig. 2). Topologically, the lobe areas can be represented by the difference of the area bounded by the individual curve. For example, let us consider two closed curves and in , oriented in counter-clockwise and clockwise direction bounding the two-dimensional subsets of , and , respectively, as in Fig. 2, then the lobe area can be represented by the set difference of the area, that is as shown in Fig. 2. Thus, studying phase space transport using lobe dynamics boils down to the fundamental problem of determining the area and defined by these two curves. We note that the area can be computed once is known by using


So, we approach this problem by developing methods that can calculate the area bounded by closed curves which are typically generated by intersecting manifolds in phase space transport.

Figure 2:  LABEL: Shows schematically lobes, pips and sips that are formed in a homoclinic tangle.  LABEL: Shows an example of two closed curves, and in , oriented in counter-clockwise and clockwise direction which bound the two-dimensional subsets of , and , respectively. LABEL: Shows schematically the area which represents a lobe area.

The main purpose of the paper is to expound on mathematical techniques useful for transport in systems that can be reduced to two-dimensional maps. We present examples in the context of time-periodic and autonomous Hamiltonian systems, however, the methods are applicable in finite-time situations such as numerical simulation and experiments in which one wants to compute the amount of fluid transported from one region of phase space to another. In the § 2-3, we derive the theory and numerical methods to separate the intersection points between the two curves into equivalence classes and derive the area of the lobes defined by the two curves. In § 4, we present an alternative method that can be used when portions of the two curves have non-transverse intersections, which arises in the case of an iterated boundary parametrized using an intersection point. In § 5, the numerical methods implemented as a software package Lober are applied to two example problems of interest in the dynamical systems literature: chaotic fluid transport in oscillating vortex pair flow and escape from a potential well in capsize of a ship.

2 Curve Area

2.1 One-Dimensional Integrals for Areas

Let us define and as the area enclosed by and , respectively and denote this using the standard Lebesgue measure for in as . The area of each region can be computed as:


by applying Green’s theorem to the vector field


Eq. (2) allows us to reduce the computation of the area of a complicated region to a one-dimensional integral over its boundary. Notice that the sign of the integral is to be reversed if the curves are oriented clockwise. Our hypothesis stating that the curves are oriented in a counter-clockwise direction is equivalent to


2.2 Numerical Methods

We assume that the curves and are given in terms of a sequence of points . We want an exact evaluation of the integral in Eq. (2) for curves given by piecewise linear segments connecting and . By defining


as parametric form of the segment, we have


As a result the exact value of the contour integral for a polygon is given by


It is to be noted that for computing the contour integral the point is repeated as , so that the curve is closed, and Eqn. 7 represents a contour integral. This numerical method forms the core of the software package Lober 111E-mail: along with the implementation of the numerical method derived in § 2,3, and 4. 11footnotetext: This is available as an open-source repository in Github at with implementation in C and additional wrapper scripts in MATLAB.

3 Intersection Points and Lobe area

Lobe dynamics is based on the geometry of a stable manifold, , and an unstable manifold, , of hyperbolic fixed points, and , their intersection points, and areas enclosed by the segments of the invariant manifolds. We note that when , is called homoclinic point, and when , and are called heteroclinic points. Following the definitions in [32, 42], a point is called a primary intersection point (pip) if the segment on connecting and and the segment on connecting and intersects only at , other than the point . A point is a secondary intersection point (sip) if there are other intersection points on the segment and . If and are two adjacent pips, then the area enclosed by the segments and is called a lobe. This is shown in Fig. 2, where are pips, and are sips. Furthermore, Fig. 1(b) shows the case when the unstable and stable manifolds are associated with two different hyperbolic fixed points. Thus, computing lobe areas require the knowledge of intersection points as accurately as possible and also identifying them as either pips or sips. Specifically, for closed intersecting curves and (see Fig. 2) encountered in phase space transport problems, we separate the set of intersection points between the two curves into classes of equivalence. For two curves corresponding to the invariant manifolds of a hyperbolic fixed point, each class of equivalence corresponds exactly to the two pips and sips on the segment of the invariant manifold. Then, we can compute the lobe area using the contour integral form of the Eqn. (7) which should be close to the lobe area given a well-resolved manifold.

3.1 Intersection Points

In this section, we assume that there are only transverse intersections of the curves. A numerical algorithm for efficiently computing the intersection points is presented below. The two curves are closed, so the number of intersection points must be even. We compute the intersection points between the two curves and . The unit tangent vector to the curve at point is denoted . For each intersection point, , we define


where is the angle between the tangents and , as shown in Fig. 3.

Figure 3: Shows the parameters involved in computing the orientation of an intersection point, , when traversing along the curves in the direction shown in dashed arrows.

The quantity represents the orientation of a given point . We note that when the intersections are transverse, the denominator of Eq. (8) is non-zero and . Fig. 4 shows two curves, their intersection points, and the value of for each point . The closed curves shown in Fig. 4 can be obtained by extracting the lobes from the intersection of stable and unstable manifolds or intersection of a cylindrical manifold with a plane. It is to be noted that the orientation of the closed curves is related to the geometry of the invariant manifolds, and as such one might have a counter-clockwise and a clockwise oriented curve.

Figure 4: Example of closed curves oriented in counter-clockwise (cyan arrows) and clockwise direction (magenta arrows) with transverse intersections at . (a) Green intersection points have . Blue intersection points have . (b) Intersection points of the same color belong to the same equivalency class. The arrows represent the one-to-one and onto relationship between the intersection points and the points joining the same equivalency class can either be a knob (black arrows) or handle (red arrows).

3.2 Classes of Intersection Points

The segments of curve between intersection points are most important to our computation, so we define and as the counter-clockwise segments of respectively and between the points and . Similarily, and are the clockwise segments of respectively and between the points and . We define the positive and negative adjacency on by


and the adjacency on as


The geometric interpretation of adjacency is that when traversing a curve in counter-clockwise or clockwise (that is, positive or negative sense) direction, the point with adjacency value of 1 is the point next to the point . The objective of this section is to determine a generalized notion of a lobe. Intuitively, lobes are bounded by “segments of the curves and turning in opposite directions on each curve”. We formalize this idea by defining the signed adjacency


which becomes for and for . Table 12, and 3 show the value of the functions , , and , respectively, for the example of intersecting curves shown in Fig. 4 .

0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
1 0 0 0 0 0
Table 1: Function for the intersection points between the curves of Fig. 4. The lines and columns correspond respectively to the first and second argument of .
0 0 0 0 0 1
1 0 0 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 1 0 0 0 0
0 0 1 0 0 0
Table 2: Function for the intersection points between the curves of Fig. 4. The lines and columns correspond respectively to the first and second argument of .
0 1 0 0 0 0
1 0 0 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
0 0 1 0 0 0
Table 3: Function for the intersection points between the curves of Fig. 4. The lines and columns correspond respectively to the first and second argument of .

In order to separate the intersection points in disjoint classes, we define an equivalency relation between the intersection points by

Definition 3.1

iff there exists a sequence of intersection points such that


Thus, the sequence of intersection points bound a domain while traversing the curve in clockwise or counter-clockwise direction. We have the following properties

Theorem 3.1 (Reflexivity of )

Proof.  Letting and in Eq. (12) gives the desired result.

Theorem 3.2 (Symmetry of )

Proof.  We notice that for any intersection point , there is one and only one intersection point satisfying . This is a consequence of the fact that the intersections are transverse and that there are not any self-intersections. There must be at least one segment of the curve leaving (based on a counter-clockwise orientation of ). If there was more than one segment leaving this point, the curve would self-intersect. Similarly, there is also one and only one intersection point satisfying . Since , there is only one intersection point ( or ) such that . We define as the unique intersection point satisfying


Notice that


and is therefore invertible. We construct the sequence using


Notice that the sequence satisfies the third condition of Eq. (12) for any length. The number of intersection point is finite, so there must be an integer such that is identical to for . We cut the sequence at the smallest possible , so there is only one repeated element in the sequence. The repeated element must be because is invertible. As a result, the sequence is the unique path from to that does not contain only and does not contain repeated points except for the endpoints. To prove the relationship, we notice that the point must be included in the unique sequence . As a result, we can take a subsequence from to in and .

Table 4: Functions and for the intersection points between the curves of Fig. 4. There are two equivalency classes for this example.
Theorem 3.3 (Transitivity of )

Proof.  The hypothesis gives us one path from to and one path from to . We can combine these two paths to go from to and . As a result, is an equivalence relation for the set of intersection points.

3.3 Lobe Area

Thus, the lobe area can be computed by using the contour integral form Eqn. (7) for the intersection points that are equivalent under . Let us define


where has been defined in the proof of Theorem 3.2. Then, the lobe area can be computed using the following

Conjecture 3.1 (Lobe area)

where we define the classes of equivalency as corresponding to subsets of containing all the elements that are equivalent under (i.e., the quotient of by ). We define the number of sets in as and we note that the equivalency classes constitute a partition of . Thus, the are disjoint and the union of their elements is . For the example intersecting curves in Fig. 4, there are equivalency classes, and . It is to be noted that using and instead of and in the definition of (Eq. 11) gives another equivalence relationship. The equivalence classes of the latter gives the area in a form identical to Proposition 3.1. We also note that Proposition 3.1 could be written with a single sum over all the intersection points. However, we prefer to keep the sum of the equivalency classes because the sum over each intersection point in the same class correspond to the area of a lobe, that is, an individual protuberance of outside .

3.4 Numerical Implementation

3.4.1 Intersection Points

The simplest approach of computing intersection points between two curves that are piecewise linear requires browsing every pair of linear segments (one on and one on ), and determining the possible intersection point between these two segments. This can become a computationally expensive operation for a manifold with million points since this approach requires solving a 2 2 linear system at each of the pair of piecewise linear segments, and hence requiring time, where is number of line segments. Furthermore, the Bentley–Ottmann algorithm uses a line to sweep across line segments with intersection points in time, and randomized algorithms can solve the problem in time. However, the application of these computational geometry algorithms is much more complicated, and in case of finite precision arithmetic leads to other challenges [39, 31]. Instead, we propose necessary and sufficient conditions to quickly check the existence of an intersection point between two piecewise linear segments.

Theorem 3.4

If the segments and intersect, then we must have


Proof.  The equation of the line passing through the two points and is


The theorem states that if the segment intersects the segment , then it must also intersect the line . As a result, the endpoints of the segment must be on opposite sides of the line and we have


Theorem 3.4 gives us a necessary condition for an intersection between two segments. Notice that we can reverse the role of each segment in Theorem 3.4 and get another necessary condition. In addition, we have

Theorem 3.5

The segments and intersect if and only if




Proof.  Theorem 3.4 directly implies one direction () of the equivalence. To prove , notice that if the two equations are satisfied, the endpoints of the first segment are on each side of the line containing the second segment. The endpoints of the second segment are also on each side of the line containing the first segment. As a result, the two segments must intersect in at least one point. The two theorems above allow for a very fast and efficient algorithm to detect intersection points. Each segment on curve is checked for intersections with each segment on curve . However, only the necessary condition given by Theorem 3.4 is checked. Only if this condition is satisfied is the second condition in Theorem 3.5 checked. If the necessary and sufficient condition is satisfied, then the intersection point is effectively computed.

3.4.2 Area bounded by a closed curve

We assume that the curves and are given in terms of a sequence of points . We want an exact evaluation of the integral in Eq. (19) for piecewise linear curves. By using Eq. (7) in the section above, Eq. (19) becomes


where and in a form suitable for numerical implementation.

4 Non-Transverse Intersections

In this section, in contrast with § 3, we relax the assumption that the two curves can only have transverse intersections. Thus, we allow non-transverse intersections of the curves, and the input curves can have common segments as shown in Fig. 5. This scenario arises in applying lobe dynamics on the Poincaré surface-of-section in 2 or more degrees of freedom, time-periodic 1-DOF system, 2 dimensional time-periodic fluid flow, and when multilobe, self-intersecting turnstiles are formed. This is shown in Fig. 6, and discussed in detail in Refs. [17, 10, 35]. In these case of multilobe turnstiles, transverse intersection cannot be guaranteed.

Figure 5: Schematic showing the boundary (green solid curve) and its pre-image (red dashed curve) with boundary intersection point (BIP) at and , respectively. When the boundaries have near-tangent intersections, computing intersection points requires expanding/shrinking the points on the curve to avoid false lobes formed by the segments of curve and .
Figure 6: (a) Schematic showing the geometry of multilobe, self-intersecting lobes. (b) An example of a multilobe turnstile that arises when considering the motion of a small body (for example, ejecta) in the field of a rotating asteroid; see Ref. [17] for details.

In the case of transverse intersections, we used an algorithm based on primary intersection points and secondary intersection points to identify the boundary of each lobe, which was then used to compute the area enclosed by the boundary. The near-tangency (when the determinant of the linear system resulting from the pair of linear segments is close to 0) of intersection of the two curves creates both a confusion in the ordering of the intersection points and the position of the lobes as well as an enormous computational difficulty to extract the actual position of the intersection. Thus, to get around the non-transverse intersection, we resort to computing a function that is amenable to the tangle geometry of manifolds and lobes.

4.1 Interior Function

Let us define the complex function


The integral of over a closed curve is equal to times the number of turns that the closed curve makes around the point (to prove that, use the fact that is analytic everywhere in the complex plane except at and use the residue theorem). We define




We note that is negative when the point is contained in the simple closed curve , and is positive otherwise.

4.2 Lobe Area

In order to extract and from the shape of the curves, we define the following quantities


A quick look at the different paths on Fig. 7 reveals that

Figure 7: Schematic view of the paths involved in the definition of , and .

Notice that the three equation above are not linearly independent and any of them can give us the expected results. However, we compute , and and use the redundancy to minimize the error on the integrals and provide an approximation of the computational error. Since




we have




And using


we find




Our algorithm combines these two results and provides the following final answer


where the last equation gives the approximate error on the computed area.

4.3 Numerical Methods

The only technical difficulty in this algorithm is the numerical computation of the functions in Eq. (28), and we present an algorithm for piecewise linear boundaries in § 4.3.1. The computation of the integrals in Eq. (30)), (31), and (32) is similar to the numerical method used for the integral in § 2.2. In addition, we present an algorithm to increase the number of points on the input curve close to the intersection points; we refer to the implementation of this algorithm in Lober as densifier, and its use is presented in § 4.3.2. The algorithm in § 4.3.2 has proved to increase the accuracy of this method without compromising the computational cost too much.

4.3.1 Computation of

In this section, we derive the exact value of the function defined in Eq. (28) when the boundary of is given as a sequence of piecewise linear segments. For each linear segment , we have


Using and , we obtain


where is the unit normal vector in the z-direction and , , .




We note that the discriminant of the denominator of Eq. (49) is


As a result, we have


We note that for , the increment to the integral is zero. This is consistent with the equation above where the right-hand term is continuous for and vanishes. This approach is similar to the point-in-polygon problem that is tackled in computational geometry (see Ref. [31]) by using the ray sweep method to check if the point is inside or outside a polygon. winding number approach. Although, our method is mathematically satisfying, it may show poor performance when compared to the more efficient method of ray passing. This is due to the integral and trigonometric function evaluations that are involved, but our approach is to adopt a method amenable to manifolds as input curves at the expense of computational performance.

4.3.2 Curve Densifier

There is a built-in densifier module in the light version of the package Lober, and which adds points on the curves close to the intersection points. The densifier can be activated by adding the parameters -DENS <nPass> <nDens> at the command line. The arguments <nPass> and <nDens> give the number of passes to be performed and the number of points to add near each intersection at each pass, respectively. The extra precision is always , where = <nDens> and = <nPass>. In other words, the precision is the same with than with . For a constant , the value of the two parameters and should minimize the computational time. Small means that fewer steps are necessary to densify the curve and can reduce the computational time. However, small usually implies large to maintain a constant . Since the extra length of the curve is , the number of points increases rapidly if is too small and lengthens the computation. So there is an optimal that minimizes computation time, and this optimal parameter depends on the local curvature of the manifold under consideration. We present a comparison of the densifier module in Fig. 8 by using this to compute the lobe area for two application problems. As demonstrated by the results, the accuracy in the lobe area does not improve dramatically beyond addition of points, and is comparable to the default densifier in the light version of Lober.

(a) OVP flow
(b) PCR3BP
Figure 8: Shows the numerically computed lobe area that used densifier option. (a) oscillating vortex pair (OVP) flow as discussed in Ref. [32], and (b) planar circular restricted three body problem (PCR3BP) as discussed in Ref. [18]. Along the x-axis, the indices denote the following densifier options: 1. -DENS 0, 2. -DENS 1 10, 3. -DENS 3 10, 4. -DENS 5 100, 5. without -DENS option which triggers the default densifier. The y-axis is the average of the entraining and detraining lobes for the parameters used in the references.

5 Application to lobe dynamics and rate of escape

In this section, we present applications of the computational approach to problems of chaotic transport in phase space that uses lobe dynamics in fluid flow and rate of escape from a potential well in a ship capsize model. In essence, using the numerical methods presented in § 3 and § 4, we demonstrate the qualitative and quantitative results that is relevant for these geometric methods of phase space transport.

5.1 Chaotic transport in fluid flow

In this section, we apply the numerical methods to quantify phase space transport in the oscillating vortex pair (OVP) flow that was introduced in Rom-Kedar et al. [32]. In this example from fluid dynamics, it is of interest to calculate the transport rate, in terms of the iterates of a two dimensional map, of species from the inner core to the outside region shown in Fig. 9(b). More generally and referring to the Fig. 1(b), this transport in phase space can be discussed using the language of lobe dynamics and transport in two dimensional maps developed in Refs. [32, 43]. We will briefly summarize the OVP flow and the transport quantities that are relevant for the numerical experiments in our study. In two-dimensional, incompressible, inviscid fluid flow, the stream-function is analogous to the Hamiltonian governing particle motion and the domain of the fluid flow is identified as the phase space. The OVP flow is generated by two counterclockwise rotating vortices under sinusoidal perturbation which can generate chaotic particle trajectories that causes regions of phase space to move across transport barriers. That is, time dependent velocity field, that is referred to as unsteady flow in fluid dynamics, can generate particle motion between the region with closed streamlines (shown as magenta in Fig. 9(a)) and the region outside the heteroclinic connection (shown as black in Fig. 9(a)). When the flow is time-dependent and has periodic forcing, a suitable approach is to construct a two-dimensional map of the flow by defining Poincaré section at a certain phase of the flow and given by