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 degreeoffreedom (DOF) Hamiltonian systems that can be reduced to twodimensional 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 codimension1 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 2DOF Hamiltonian systems that can be reduced to twodimensional map by using suitable Poincaré section, or return map, or timeT 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 twodimensional phase space of a Poincaré map (in general, this can a return map or a timeT 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) 
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, selfintersecting, 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:

Constructing a functional of the nonlinear system of vector field which measures area between the manifolds as parametrized by time. This is a semianalytical approach and only remains valid near small perturbations, for example Melnikov method for mixing, stirring, optimal phase space flux, and actionintegral method
[25, 37, 4, 5, 29, 6].
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 counterclockwise and clockwise direction bounding the twodimensional 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
(1) 
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.
The main purpose of the paper is to expound on mathematical techniques useful for transport in systems that can be reduced to twodimensional maps. We present examples in the context of timeperiodic and autonomous Hamiltonian systems, however, the methods are applicable in finitetime 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 § 23, 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 nontransverse 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 OneDimensional 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:
(2) 
by applying Green’s theorem to the vector field
(3) 
Eq. (2) allows us to reduce the computation of the area of a complicated region to a onedimensional 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 counterclockwise direction is equivalent to
(4) 
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
(5) 
as parametric form of the segment, we have
(6) 
As a result the exact value of the contour integral for a polygon is given by
(7) 
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 ^{1}^{1}1Email: shiba@vt.edu along with the implementation of the numerical method derived in § 2,3, and 4. ^{1}^{1}footnotetext: This is available as an opensource repository in Github at https://github.com/Shibabrat/curve_densifier 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 wellresolved 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
(8) 
where is the angle between the tangents and , as shown in Fig. 3.
The quantity represents the orientation of a given point . We note that when the intersections are transverse, the denominator of Eq. (8) is nonzero 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 counterclockwise and a clockwise oriented curve.
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 counterclockwise 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
(9) 
and the adjacency on as
(10) 
The geometric interpretation of adjacency is that when traversing a curve in counterclockwise 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
(11) 
which becomes for and for . Table 1, 2, 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 
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 
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 
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
(12) 
Thus, the sequence of intersection points bound a domain while traversing the curve in clockwise or counterclockwise direction. We have the following properties
Theorem 3.1 (Reflexivity of )
(13) 
Proof. Letting and in Eq. (12) gives the desired result.
Theorem 3.2 (Symmetry of )
(14) 
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 selfintersections. There must be at least one segment of the curve leaving (based on a counterclockwise orientation of ). If there was more than one segment leaving this point, the curve would selfintersect. 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
(15) 
Notice that
(16) 
and is therefore invertible. We construct the sequence using
(17) 
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 .
Class  

1  
1  
2  
2  
2  
2 
Theorem 3.3 (Transitivity of )
(18) 
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
(19) 
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)
(20) 
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
(21) 
Proof. The equation of the line passing through the two points and is
(22) 
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
(23) 
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
(24) 
and
(25) 
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
4 NonTransverse Intersections
In this section, in contrast with § 3, we relax the assumption that the two curves can only have transverse intersections. Thus, we allow nontransverse 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é surfaceofsection in 2 or more degrees of freedom, timeperiodic 1DOF system, 2 dimensional timeperiodic fluid flow, and when multilobe, selfintersecting 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.
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 neartangency (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 nontransverse 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
(27) 
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
(28) 
and
(29) 
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
(30) 
(31) 
(32) 
A quick look at the different paths on Fig. 7 reveals that
(33) 
(34) 
(35) 
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
(36) 
and
(37) 
we have
(38) 
and
(39) 
And using
(40) 
we find
(41) 
and
(42) 
Our algorithm combines these two results and provides the following final answer
(43) 
(44) 
(45) 
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
(46) 
Using and , we obtain
(47)  
(48) 
where is the unit normal vector in the zdirection and , , .
(49) 
where
(50) 
We note that the discriminant of the denominator of Eq. (49) is
(51) 
As a result, we have
(52)  
(53)  
(54) 
We note that for , the increment to the integral is zero. This is consistent with the equation above where the righthand term is continuous for and vanishes. This approach is similar to the pointinpolygon 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 builtin 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.
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 RomKedar 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 twodimensional, incompressible, inviscid fluid flow, the streamfunction 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 timedependent and has periodic forcing, a suitable approach is to construct a twodimensional map of the flow by defining Poincaré section at a certain phase of the flow and given by
Comments
There are no comments yet.