Efficient manifolds tracing for planar maps Free

In dynamical systems, the determination of unstable periodic orbits and their invariant manifolds is a problem on its own and is fundamental to understanding the underlying structure of the chaotic orbits of planar maps.¹ These manifolds provide the skeleton for chaotic dynamics and are relevant in physical applications for the interpretation of the transport induced by non-integrability.^2,3

Physically, invariant manifolds determine the geometrical features of advection and delimit the extension of irregular flows. For instance, in magnetically confined plasmas, a small non-axisymmetric magnetic field induces Hamiltonian chaos near the plasma edge of single-null tokamak discharges.^4–8 In this situation, the invariant manifolds of the magnetic saddle can be related to the heat flux patterns where the confinement chamber interacts with the plasma.^9,10 In analogous systems, time-dependent perturbations lead to orbit transfers between otherwise isolated regions of bistable non-linear oscillators.¹¹ Here, invariant manifolds determine the transport channels between the phase space domains.

From the dynamics point of view, manifold tracing is required for determining the geometry of the basins of attraction of chaotic sets,¹² and can be used to identify homoclinic and heteroclinic intersections which are contained in non-attracting invariant sets that play a fundamental role on the description of the dynamics of typical chaotic orbits in the phase space.¹³ Although global quantities like the Lyapunov exponents allow us to determine numerically the bifurcation parameters, they do not, in general, provide details on the type of bifurcation that takes place. In this situation, precise invariant manifolds tracing determines clearly the topological mechanisms that take place during the bifurcations and accounts for physically relevant features, like regions of enhanced transport or the destruction of transport barriers.¹⁴ 

In some situations, the geometry of the manifolds can be estimated through the mapping of a large collection of orbits close to the saddle, a procedure similar to the one used to determine chaotic saddles.^15–18 However, without an ordering scheme and refinement, this method is computationally expensive and limited in resolution.

In this work, we introduce an exact and an approximation method to calculate the invariant manifolds of periodic saddles of planar maps. The maps in question can be explicit or induced by three-dimensional flows. Both calculations are based on the manifold decomposition in primary segments introduced by Hobson.¹⁹ The exact method relies on the efficient calculation of primary segments from a single seed segment near the saddle, while the approximated method is based on the determination of reliable interpolants for arbitrary primary segments. The latter is specially relevant for maps induced by flows, where lengthy numerical integrations are required to produce a single iteration of the Poincaré map. This allows us to deal with the increase in length of the manifold segments due to the stretching and folding mechanism inherent to the chaotic dynamics.

The approximation method involving the primary segments interpolation has been discussed before by Hobson¹⁹ and Krauskopf,²⁰ and more detailed implementations were made by Goodman and Wróbel²¹ and references therein. In this work, we introduce a versatile interpolation approach based on the curve decomposition in bare and fine details, which allow us to focus on different aspects of the manifold separately. The bare details of the curve are determined by a discretization procedure that determines an appropriate set of nodes containing most of the curve information, while the fine details are contained in a normal displacement function or shape function which ensures the smoothness of the manifold. The presented normal displacement description is new to our knowledge and can be easily implemented, containing a single set of adjustable parameters which can be fixed from intuitive geometrical conditions, complementing other approaches based on Catmull-Rom splines.²¹ The resulting approximation method gives the invariant manifold as a continuous parametric curve from which we can measure the distance to the discretized manifold obtained from the exact calculation. This allows us to estimate the approximation error that, in general, is well below an estimated baseline. The approximation procedure is numerically stable and more efficient than the exact one, making it ideal for fast and precise calculations.

The methods presented in this work are specific for one-dimensional manifolds of planar maps. Other sophisticated methods are available for tracing manifolds for continuous dynamical systems,²² which can then be used to intersect the surface of section instead of using a Poincaré map. Also, parametrization methods^23,24 are available for obtaining a patched manifold in discrete systems with arbitrary precision. In contrast to these methods, our simpler approach can be easily implemented by a non-specialist and is robust enough to deal with fairly complex manifolds. This work is intended to enable the reader to trace exact and approximated manifolds with a reliable numerical procedure, without dipping into the subject of numerical interpolation or parametrization functions.

The paper is organized as follows: in Sec. II, we present a review on invariant manifolds and their decomposition in primary segments. In Sec. III, we introduce an approximation method to obtain the zeroth segment of the manifold to machine precision. In Sec. V, we present an exact calculation method and in Sec. VI, we introduce its approximated version, where the interpolation procedure is detailed. In Sec. VII, we show an example application of the manifold tracing for bi-stable systems where the transition from homoclinic to heteroclinic chaos is used to explain the transition from alternating chaos to global chaos. In Sec. VIII, we present our conclusions and perspectives.

In this review section, we want to introduce the concept of invariant manifolds for planar maps and show their decomposition in primary segments. We also show, formally, the method to obtain any primary segment from a single known segment which will be discussed in the following section.

Without loss of generality, we will develop our discussion around the Poincaré maps of three-dimensional flows, since they represent a situation in which we do not have access to a closed-form expression for the planar map of interest. Consider the dynamical system

where $x∈R3$ and $f:R3→R3$⁠. The solution $ϕt(x0)$ of this differential equation takes an initial condition $x0∈R3$ and maps it to a new point $x∈R3$ for a given $t∈R$⁠. Now, consider an orientable two-dimensional surface $Σ$ which is transverse to the flow everywhere, so that orbits starting in $Σ$ must eventually return to $Σ$⁠. Then, the Poincaré map $T$ is a diffeomorphism $T:Σ→Σ$ such that $y=T(x)$⁠, where $x,y∈Σ$ and $y=ϕτ(x)$ for the smallest $τ>0$ such that the flow traverses $Σ$ always in the same direction (Fig. 1). Analogously, the inverse map$T−1$ is such that $x=T−1(y)$⁠, where $x=ϕτ(y)$ for the largest $τ<0$⁠.

A fixed point $x∗∈Σ$ satisfies $T(x∗)=x∗$⁠, which means that the orbit starting in $x∗$ is closed or periodic. For the case of interest, $x∗$ is a saddle point, i.e., the Jacobian matrix of $T$ has real eigenvalues ${λs,λu}$⁠, such that $|λs|<1$ and $|λu|>1$⁠.

The unstable manifold of $x∗$⁠, $Wu(x∗)$⁠, is the collection of points that converge to $x∗$ under $T−1$⁠, and the stable one $Ws(x∗)$ is the collection that converge to $x∗$ under $T$⁠. Formally, we can write¹³ 

For the numerical calculation of the invariant manifolds, it is useful to develop an ordering scheme that allows us to calculate them recursively, for this we resort to the concept of primary segments.¹⁹ A primary segment based on $x$ is a continuous subset of the manifold $Ps,u(x)⊂Ws,u(x∗)$⁠, containing every point in $Ws,u(x∗)$ between $x$ and $T(x)$⁠, including $x$ but not $T(x)$⁠. Notice that, if $y∈Pu(x)$⁠, then $T(y)∉Pu(x)$⁠, but clearly $T(y)∈Pu[T(x)]$⁠, i.e., every point inside $Pu(x)$ has an image in the next primary segment and a pre-image in the previous one, and this is also valid for the base point $x$ (Fig. 2). Without loss of generality, we will continue this discussion in terms of the unstable manifold, since the definitions and methods are analogous for the stable one.

The unstable manifold can be written as the union of all primary segments based on the images of an arbitrary $x∈Wu(x∗)$⁠,

Since every point in $Pu(x)$ has an image in $Pu[T(x)]$⁠, we can say that the primary segment based on $T(x)$ is given by

where $T[A]$ is the set obtained after mapping every element of $A$⁠. This implies that the manifold can be obtained by joining subsequent mappings of an arbitrary primary segment

Consequently, if we know every point in any primary segment, we can map them forward and backward repeatedly to obtain the desired extent of the manifold. This is the theoretical basis of interpolation methods,^19,21 and it is valid for both conservative and dissipative diffeomorphisms.

Now, we need to determine at least one primary segment of the manifold. This can be done to any desired precision in the region near the saddle point and, here, we illustrate this with a third order parametric form. Consider an arbitrarily small displacement vector $δx$ around $x∗$⁠. To a second order, we can write

where $JT(x∗)$ and $HT(x∗)$ are the Jacobian matrix and Hessian tensor of the map $T(x)$ at the saddle point. The idea here is to determine a function $f(u)$ that determines the normal displacement of the manifold at position $u$ along the tangent line $lu$⁠, touching the manifold at $x∗$⁠. In this representation, the parametric form of the local manifold takes the vector form

where $u^$ is the unit unstable eigenvector of $JT(x∗)$⁠, satisfying

with $λu>1$⁠, and $v^$ is perpendicular to $u^$⁠. Using Eq. (9), it can be shown that a vector of magnitude $u$ aligned with the unstable direction $u^$ gets stretched by a factor $λu+κu$⁠, in the unstable direction, with the parameter $κ$ given by

where $ui$ are the coordinates of $u^$ in some set of orthogonal coordinates, and $Ti,jk=∂2∂xj∂xkTi(x1,x2)x∗$ are the map derivatives with respect to these coordinates.

If we represent the phase space points by their displacements $(u,v)$ along $u^$ and $v^$ with respect to the fixed point $p∗$⁠, it can be shown that the normal displacement function satisfies

This comes because the effect of $T$ over a point in the manifold is moving it by a factor $λu+κu$ along the unstable direction and displacing it normally to a new distance given by $f(λuu+κu2)$⁠. Provided that we are close to the saddle, we can use (7) in the new coordinates to write-down the LHS of (11), which leads to a consistency relation for $f$⁠,

In practice, we can use a polynomial expansion for $f$⁠,

where $1$ and $u$ are absent because of the chosen coordinates. Provided that (12) is valid for any small $u$⁠, we have that the coefficients resulting from replacing $f(u)$ in (12) must vanish independently. This allows us to relate the coefficients of $f$ with the local properties of the map, i.e.,

In the  Appendix, we provide additional information on how to obtain the Hessian elements in (14) for an arbitrary planar map in arbitrary coordinates. With this, we have the information required to write explicitly a parametric form of the manifold near the saddle point $x∗$ and, consequently, we are able to determine, to a leading order, one primary segment of the manifold that we will call the zeroth segment,

for some $u0<<1$⁠. Numerically speaking, a representation of the local manifold can be determined to machine precision by making the highest order term of the manifold expansion to be below the roundoff error. In practice, we can check the performance of the local approximation by applying the map repeatedly over an initial condition obtained by $Wuloc(x∗,u0)$ for some small $u0$⁠. In Fig. 3, we compare the performance of a linear and cubic expansions of the local manifold for one saddle of the Henon map²⁵ with the classical parameters $aH=1.4$ and $bH=0.3$⁠.

In regular situations, a linear approximation of the local manifold may be sufficient,²⁰ but having higher orders in the expansion of $f(t)$ have important advantages in practical applications, given that a larger primary segment can store more preimages required to obtain numerically the global manifold. For instance, the primary segment between the $4th$ and $5th$ iterates of the cubic Henon manifold can generate preimages to machine precision along a portion of size $10−3$⁠, while the linear manifold can only generate preimages with precision $10−13$ between the initial condition and its first iterate, in a portion of size $10−6$⁠. Clearly, there is an important gain when using higher order approximations of the manifold, in this example, a $1000$-fold increase in available space with $1000$-fold increase in precision (assuming machine precision around $10−16$⁠).

Numerical indistinguishability between the actual manifold and its functional approximation guarantees that they will be equally affected by roundoff errors due to the application of the map itself, i.e., approximation errors will not propagate numerically if they were numerically nonexistent. In this context, general methods based on parametrization functions valid for high-dimensional systems can be used as well.²⁶ Generalizations for continuous systems are also available,²³ allowing us to approximate the manifold along the flow instead of in a cutting surface (the Poincaré section).²⁴ 

Provided that we identified an appropriate parametric form of the zeroth segment, the unstable manifold can then be written as the following semi-infinite union:

where $Pi+1=T[Pi]$ and $x∗x~$ is the union of inverse mappings of $P0$ which can be approximated by the parametric form $Wuloc(x∗,t)$ for $t∈(0,t0)$⁠.

With this simple definition, we can compute a finite representation of the invariant manifold to the desired extension, provided we know how to calculate $T[Pi]$ and we know one segment $P0$⁠. Formally speaking, if we know how to compute $T(y)$ for any $y∈Pi$⁠, then we are able to compute $T[Pi]$⁠, but mapping a set of uniformly distributed points in $Pi$ will not, in general, lead to evenly distributed set of points in $Pi+1$⁠.

As we do not know a priori which distribution of points in $P0$ will lead to a uniform discretization of $Pn$⁠, an interim solution consists in discretizing $P0$ uniformly with a large number of nodes, although, eventually, the stretching and folding mechanism will overcome the excess points, and the segments will become poorly resolved.

In Fig. 4, we show a short portion of the unstable manifold of the period-one saddle of the Chirikov-Taylor map,²⁷ 

for $k=1.5$⁠. Although the resulting curve appears smooth, close inspection shows the undesired effects of uniform discretization in $P0$⁠. The segment stretching leads to a complex mix of over-resolved and under-resolved portions of the manifold. In Sec. V, we introduce an adaptive method that allows us to determine the discretization of $P0$ that leads to a well resolved discretization of an arbitrary $Pn$⁠.

In this section, we introduce a refinement method for producing a suitable exact discretization of the manifold which satisfies some predefined resolution criteria by limiting the distance between subsequent nodes and the angles between consecutive secant lines. The word exact is used here to refer to the procedure by which nodes in an arbitrary segment are obtained from the zeroth segment, given that it can be determined to machine precision. More precisely, the only source of error in this case is the application of the mapping itself. For explicit mappings, it is the accumulated roundoff error, while for Poincaré maps, it is the inherent error of the numerical integration routine.

The exact method is then a procedure to obtain the most reliable discretization of the manifold for a prescribed finite precision, while the actual curve between each pair of nodes must be approximated by some suitable method.

Consider the following initial discretization of $P0$⁠:

where the “$→$” symbol is used instead of “=” to indicate that we are only representing the continuous set $P0$ by a discrete set of points. The initial set of points can be obtained by evaluating

for $i=1,2,…,M$⁠, $M$ is a small number, usually $<10$⁠. Then, the induced discretization of the segment $Pn$ is

where

As mentioned before, the resulting points of $Pn$ will be unevenly spaced due to the non-uniform stretching of the primary segments, but this can be solved by adding new points to $Pn$ at specific locations. For instance, if some resolution criterion is violated between $xn,k$ and $xn,k+1$⁠, we need to insert a new manifold point $xn′$ between these points. Since we do not have a parametric representation for $Pn$⁠, we need to obtain the new points from initial conditions in the first segment $P0$⁠, where we do have a parametric form given by (8). Then, we calculate $x0′$ between $x0,k$ and $x0,k+1$ in $P0$ and calculate $xn′=Tn(x0′)$⁠, which is the required refinement point (Fig. 5). This procedure can be iterated until the discretization of $Pn$ satisfies the desired resolution requirements.

The refinement procedure must result in evenly spaced points for low curvature regions and a higher density of points in high curvature regions. Such a scale-invariant resolution is desired if we require a reliable description of densely packed portions of the curve, or if we need to study the intersection between manifolds. In Fig. 6, we show the quantities involved in resolution conditions. In general, the length of the chord $li$ connecting the nodes $xi$ and $xi+1$ is required to be below some small value $lc$⁠. This chord belongs to an infinite secant line that forms and angle $θi$ with the secant line containing the chord $li−1$⁠. Analogously, the inter-secant angle $θi$ must be below some predefined value $θc$⁠. Provided that $lc$ and $θc$ are sufficiently small, the arcs connecting the points $xi−1,xi$⁠, and $xi+1$ will not differ much from straight lines. With this, we establish that most of the geometric information is contained in the nodes instead of the arcs connecting them. Then, the small details contained in the parameters that determine the arcs are used to guarantee the manifold differentiability.

Upon violation of the resolution criteria, a refinement procedure must be carried out as illustrated in Fig. 7. In short, a new point is introduced between $xi$ and $xi+1$ if the chord $li$ is larger that the critic distance $lc$⁠, or if $li$ is the longest chord in an angle that exceeds some small critical value $θc$⁠. This simple procedure and the rules of motion in Fig. 7 prevent clustering of new nodes and distribute the new entries evenly along the manifold. The procedure is iterated until the criteria are satisfied for all nodes in the segment $Pn$⁠, then these nodes are mapped to obtain an initial discretization of $Pn+1$ and so on.

This method differs from others in that we do not try to extrapolate the location of new nodes based on the information of previous nodes,^19,20 we start with an initial set of nodes generated by mapping the nodes of the previous segment, then perform the refinement operations following Fig. 7. A similar method is mentioned by Goodman and Wróbel,²¹ but it adds two points instead of one and does not discuss on movement rules, which in our experience have important consequences for the efficiency and reliability of the refinement.

In Fig. 8, we show the same portion of Fig. 4 calculated with the mapping-refinement method. The calculation was done so that the number of nodes in $P0$ is the same for the uniform discretization example, and the number of calls to the map function is the same for both calculations. In the inset region of Fig. 8, it is clear that the resulting nodes are evenly spaced in regions of low curvature and accumulate in the high curvature regions as required. Although the results from the mapping-refinement method have no error intrinsic to the manifold tracing, they have some inherent error sources, as the roundoff error when applying explicit maps or the error of numerical integration routine when using a Poincaré map.

Before introducing the interpolation-mapping approach, let us discuss a relevant limitation of the presented method. Since all the points in the manifold are seeded in $P0$⁠, the mapping-refinement method is numerically intensive; more precisely, if $Mn$ is the number of points required to resolve the $n$th segment, the total number of calls to the map function when tracing $n$ segments of a manifold is

Assuming that the stretching mechanism leads on average to an exponential growth in the primary segments, we obtain

for an appropriate $λ>0$⁠. Consequently, the presented method can be computationally expensive for manifold calculations involving a large number of primary segments, although in most situations a few tens of segments are sufficient to obtain sufficient information on the manifold geometry. A more fundamental problem lies in the finite-precision of the numerical representation, which limits the number of different initial conditions that can be represented in $P0$⁠, and can be insufficient to resolve a given $Pn$ due to the exponential growth of $Mn$⁠. This is why it is important to use high order interpolants for the zeroth segment, as this allows us to work with a large primary segment. In practical situations, the stretching mechanism can be sufficiently large for two neighbor initial conditions in a linear representation of $P0$ separated by the minimum representable numerical difference in coordinates to become separated after a few mappings by a distance that exceeds the resolution limit $lc$ in the phase space.

In Sec. V, we introduced a refinement method to produce a suitable discretization of the exact manifold to any desired resolution and extent. Despite our control over the calls to the map function and a precise definition of the refinement regions, the resulting procedure can be computationally intensive (but far less intensive than a non-refining procedure). To reduce the computational cost we can move to an approximation paradigm,¹⁹ where the refinement of segment $Pn+1$ is based on the approximation of initial conditions in segment $Pn$⁠, instead of initial conditions in $P0$⁠.

In this section, we concentrate on producing an approximated continuous representation for $Pn+1$ based on an approximated representation of $Pn$⁠. This is formally closer to the original composition presented in (16), but, in practice, we use the refinement methods presented in Sec. V to define an iterative scheme that allows us to produce reliable interpolants from a previous segment in the manifold.

Assume that we have a planar parametric curve $P¯n(s):R→R2$ for $s∈[sl,sr]$ that approximates the points in the $n$th segment $Pn$ to a given precision $ε$⁠. This is

In other words, the points in $P¯n(s)$ differ from those in $Pn$ by at most a distance $ε$⁠. To approximate the segment $Pn+1$⁠, we can infer its behavior from an appropriate set of nodes

where $s0$ and $sN$ are the parameters corresponding to the end-points of the segment $Pn+1$⁠. The set $s0,s1,…,sN$ is chosen so that the curve between any pair ${xi,xi+1}$ cannot differ much from a straight line. The procedure to do this is analogous to that presented in Fig. 7, where a new node between $xi$ and $xi+1$ is added when the resolution criteria are not met. The new point is obtained by parametric bisection

After determining an appropriate set of nodes for the whole domain of $Pn+1$⁠, the interpolant curve to approximate it is a piecewise vector function

where $fi:R→R2$ and $γi$’s are adjustable parameters that determine the local properties of the curve and can be determined by imposing geometrical constrains on the interpolant pieces $fi$⁠.

Consider the following expression for the position vector along the interpolant of the $i$th arc:

where $x→i$ is the position vector of node $xi$⁠, $l→i=x→i+1−x→i$⁠, and $t∈[0,1]$ is the normalized position along the chord line joining $xi$ and $xi+1$ (see Fig. 9). The function $hi(γi,t):R→R$ is proportional to the normal distance between the interpolant curve $fi(s)$ and the chord $xi,xi+1¯$⁠, and satisfies $hi(0)=hi(1)=0$⁠, which guarantees that the curve passes through the nodes $xi$ and $xi+1$⁠. Notice also that the L.H.S. in (29) depends on the global parameters $s$⁠, and the R.H.S. depends on the normalized parameter $t$ so that there is a function $s=g(t)$ that allows us to make the transformation from the local parameter $t$ to the global one $s$⁠. The interesting part of this arc representation is that the bare geometrical aspects of the curve are contained in the first half of the expression $x→i+tl→i$⁠, and the fine tuning is carried by the shape function$hi(γi,t)$⁠, in the second part, which contains the error adjustable parameters $γi$⁠. Notice also that a good choice of the partition nodes results in $hi(γi,t)<<1$⁠, everywhere in the interval, guaranteeing that the main geometry is well described by the nodes themselves and only small (but extremely important) details are described by $hi(γi,t)$⁠.

By definition, the arc representation (29) is continuous across the nodes $xi$⁠, but we also require the curve model to be differentiable everywhere. Since the interpolant is smooth inside the arcs, we concentrate on the nodes and require

for all $si$ in the partitioning. In terms of the local representations of the interpolant, this becomes

where $t$ and $t′$ are the corresponding local parameters of the arcs $i−1$ and $i$⁠. Because the global parameter $s$ has the same meaning across the node, the differential $ds$ is the same in both sides, and we are left with a simpler relation that does not depend on the global parameter

The relation $dt′/dt$ can be determined by the ratio of the projections of a line element $dr→$ along the secant lines $xi−1xi¯$ and $xixi+1¯$⁠, leading to

where $αi$ and $βi−1$ are the angles between the tangent line at $xi$ and the secant lines. In other words, $αi$ and $βi$ are the angles between the $i$th arc and its corresponding chord line at the endpoints (see Fig. 9). Writing $βi−1$ as the difference between $αi$ and the inter-secant angle $θi$⁠, we obtain a consistent equation for the differentiability at the node $i$⁠,

where the angles $θi$ can be determined from the position of the nodes alone. For a suitable choice of the shape function$hi(γi,t)$⁠, the chord-tangent angles $αi$ determine the overall behavior of the model curve and can be considered as the fundamental parameters of our interpolant.

Using the explicit arc representation in (29) and the differentiability condition (34), we can relate the angles $αi$ with the derivatives of the shape function at the end-points

Requiring the shape function $h(γi,t)$ to be the lowest order polynomial, vanishing for $t=0$ and $t=1$⁠, and allowing for different left and right derivatives, we obtain a cubic function as the simplest expression for the shape function

where $ai=hi′(0)$ and $bi=hi′(1)$ are the end-point derivatives, and $bi$ can in turn be related to the parameters of the next shape function using the conditions obtained from differentiability in (35)

where $mi+1=tan⁡θi+1$⁠. Given that the inter-chord angles $θi$ are the resulting features of the discretization procedure, they are not adjustable parameters and are kept fixed. This leaves us with the set $a→={a0,a1,…,aN}$ which can, in principle, be adjusted to obtain an interpolant $P¯n+1$ which is closest to the primary segment $Pn+1$⁠. However, to adjust this set, an error functional must be minimized, and this requires obtaining more information about the curve being modeled. For instance, we can calculate a set of intermediate points $yi=T[P¯n(si′)]$ where $si′∈(si,si+1)$ and minimize the functional

by performing subsequent variations $δa→$ obtained from a Levenberg-Marquardt procedure.²⁸ As expected, this results in an interpolant that passes through every node and is closest to the intermediate points $y→i$⁠. However, there can be an arbitrary distance between the interpolant and the primary segment $P¯n+1$ in every other location without affecting the error functional. In fact it was observed that an appropriate guess of the set ${a0,a1,…,aN}$ performed better at arbitrary locations than the numerically optimized set, even when additional optimization constraints were applied, like requiring the interpolant arc-length to be small or limiting the domain of variation of the $αi$’s.

A suitable guess consists in requiring the tangent line on each node to bisect the inter-secant angle. In other words, $αi=θi/2$⁠, resulting in the following relations between parameters:

with $mi=tan⁡θi$⁠. This results in a very smooth curve passing through all the nodes, but as discussed before, the success of this approximation relies on leaving a minimum amount of information on the shape function. Most of the information must be contained in the node positions, and the interpolant between two nodes must differ very little from a straight line, i.e., the values of $ai$ must be very small. Notice that this is automatically satisfied for a good set of nodes because all inter-secant angles satisfy $θi<θc$ with $θc<<1$⁠.

We can estimate the maximum distance from the interpolant curve to the corresponding secant line in terms of the resolution parameters $lc$ and $θc$⁠. Provided that $θc<<1$⁠, we have

where $θc$ is measured in radians and $hc=h(1/2)$ for $a=b=tan⁡θc$⁠. This gives us a baseline to measure the interpolation error, which must be small in units of $Hc$⁠.

In Fig. 10, we compare the approximated manifold and the exact one in a $1×1$ region of the phase space for the Chirikov-Taylor map with $k=1.5$⁠. For this comparison, we traced $26$ primary segments with $dc=0.01$ and $θc=3.0°$ for both the exact and approximated calculation. This gives us a distance baseline of $Hc≈6.5×10−5$ so that the local interpolation error must satisfy $ε<<Hc$⁠. In the scale of Fig. 10, differences between the manifolds are not observable, and in fact they are not observable in any low or moderate-curvature region, because the distance between consecutive nodes of the manifold is much larger than the separation between the curves. To depict any difference we must look into very sharp corners, where any tiny amount of compression along the manifold amplifies the differences between the exact and approximated manifold.

To have a better perspective of the approximated manifold error, we can measure the distance between an exact node and its expected position from the manifold interpolant. This can be done by finding the normal projection of the exact node along the secant line between two nodes of the approximated manifold. This gives us the value of $t$ in (29), and we can calculate the corresponding interpolant position.

In Fig. 11, we show the distance between the exact nodes and their corresponding estimation by the interpolants that compose the approximated manifold. Because the first segment of the approximated and exact methods is the same, we expect smaller errors for the first segments. However, the correlation between the nodes of the exact calculation and approximated method is quickly lost because of the difference in the construction methods, then the distance from the exact nodes to their estimated positions gives us a good estimate of the local error. An important feature of Fig. 11 is that the interpolant error does not grow along the manifold as a consequence of mapping approximated new nodes based on an interpolant curve. This indicates that the approximation method is stable, i.e., small errors introduced by the interpolation of new nodes get damped due to the compression perpendicular to the manifold.

For this example, the manifold error is quite small, as can be seen in Fig. 12, where the statistical distribution of distances of Fig. 11 is shown. The error of the interpolants that compound the approximated manifold concentrate around $10−9$ and below $10−17$⁠, and error values close to the baseline $Hmax≈6.5×10−5$ show a very small probability. Remarkably, these results were obtained after setting the adjustable parameters ${α0,α1,…}$ only by geometrical requirements. It is important to realize that errors below $10−17$ must not be considered to be meaningful, since they are below machine precision, and their numerical representation is not reliable. However, it is important to keep them in the error histogram because they account for the events where the approximated and the exact manifolds were numerically identical.

The resolution criteria led to $106$ calls to the map function for the exact calculation and $105$ calls for the approximated one. This factor ten reduction in the number of calculations is due to the seeding of a new initial condition in the segment previous to the segment under refinement, instead of the zeroth segment. The number of calculations involved in the approximated calculation can be approximated by

so that, for a large stretching rate $λ$ and large number of primary segments $n$⁠, the number of calls to the map approaches

then, for a finite number of segments, the ratio to the number of calls for the exact calculations satisfies

where $κ=(1−e−λ)−1>1$⁠. This rapidly decreasing function of the number of segments $n$ guarantees a sustained growth in efficiency as the number of primary segments increases.

Now that we have illustrated the approximated manifold calculation for an explicit map, we consider a situation where the computational efficiency is more relevant, a continuous time dynamical system. Consider, for instance, the conservative Duffing oscillator with the Hamiltonian function

where $q$ and $p$ are conjugate variables and satisfy the Hamilton equations. Physically, $q$ and $p$ are the position and momentum of a particle in a double-well potential subjected to a periodic force field. In the unforced situation, particles with energy $E<0$ are restricted to one side of the well, but in the forced situation, the external field enables the transition from one well to the other for energies moderately below zero. This problem has interest in modeling the increase of the transition rate over potential barriers induced by external monochromatic fields.¹¹ Moreover, the corresponding deterministic diffusion in the chaotic region competes with quantum-mechanical tunneling between wells.^29,30

In the subject of magnetically confined plasmas this problem is equivalent to that of an asymmetric single-null tokamak discharge.³¹ Here, small asymmetries cause some magnetic field lines to become chaotic around the plasma column and the external divertor coils, leading the charged particles to encounter the tokamak chamber.^4,5 The use of resonant magnetic perturbations to break the magnetic invariants in the separatrix region has important applications in the control of plasma edge instabilities and plays an important role in modern tokamak operation.

Broadly, the potential barrier penetration can be understood through the formation of a chaotic layer around the separatrix of the well, which overcomes the potential barrier of the integrable case, allowing chaotic orbits to wander in both sides of phase space. However, there are interesting situations where the transit between sides of the well is virtually suppressed for orbits in the chaotic layer during long extents of time, and manifold tracing provides relevant insight into the mechanisms involved in transport or lack of it.

In Fig. 13, we show $30000$ cycles of the stroboscopic position of the Duffing system for an arbitrary orbit in the chaotic layer with different forcing amplitudes. For $ϵ=0.025$⁠, the orbit spends a long time on each side of the well with shorter interludes of global motion. As we increase the amplitude, the times spent in local motions become shorter, and global motion begins to dominate. For instance, at $ϵ=0.05$⁠, periods of localized motion are very short compared to those of global motion.

To understand these transitions in a geometrical basis, we need to know which invariants control the motion in the chaotic region. From structural stability,¹³ it is expected for the perturbed system to contain a saddle $xc∗$⁠, which is a time-dependent version of the unperturbed saddle at $(0,0)$⁠. The invariant manifolds of this saddle are responsible for driving the chaotic orbits from one side of the well to the other, but this is not the only period-one saddle in the chaotic region. For $ω=1.5$⁠, there are two resonant tori near the separatrix of the double well, one for each side of the well. For $ϵ=0.025$⁠, they are destroyed by the perturbation creating period-one K.A.M. islands surrounded by a thin chaotic layer driven by the invariant manifolds of a helical saddle $xh∗$⁠. The chaotic layers of the central saddle and the islands are in contact, but their interaction depends on the perturbation strength.

In Fig. 14, we show a single orbit in the merged chaotic layers and an inset depicting the central and right period-one saddles with some of their invariant manifolds. The periodic orbits $xc∗$ and $xh∗$ were determined to $14$ significant digits using a modified Levenberg-Marquardt algorithm.²⁸ For $ϵ=0.025$⁠, a finite tracing of the invariants shows no entanglement between the manifolds of the central and the helical saddle. The lobes of the stable and unstable manifolds of $xh∗$ develop around $xh∗$ so that chaotic orbits near the helical saddle tend to remain so, leading to a slow diffusion from the chaotic layer around the island to the global layer responsible for transitions, as if they were only weakly connected.

When the forcing amplitude is increased, the lobes of the stable and unstable manifolds of $xh∗$ develop around $xc∗$⁠, and vice versa, causing a heteroclinic connection between these saddles and a strong entanglement between their manifolds which results in a full mixing of the chaotic layers (Fig. 15). This transition from homoclinic to heteroclinic connection explains the transition from the alternating localized chaotic motion to the global chaotic dynamics, and leads to the suppression of the local motion features of around the period-one islands. Provided that there are transitions for the $ϵ=0.025$ case, there must be some heteroclinic connection between the saddles, and such can be observed after a rather extensive manifold tracing; however, this is a second-order feature and the homoclinic connections determine the most relevant features of the weak forcing situation.

In this paper, we introduced a mapping-refinement approach to discretize the exact invariant manifolds of planar maps. This intensive method prevents redundant calculations and makes the best use of every new orbit resulting in a lower computational cost when compared to non-refining approaches. An intensive non-refinement technique with the same computational cost results in a non-smooth representation of the manifold with under-resolved and over-resolved portions.

Then, we introduced an interpolant-mapping approach to approximate the invariant manifolds. Our interpolation method was based on the curve decomposition in bare and fine details, where the fine details were contained in a normal displacement function or shape function. This curve decomposition is new to our knowledge and does not involve any parametric optimization stage, which eases the method implementation.

With the approximation method, the efficiency was greatly increased with a small precision penalty due to the interpolation of the primary segments for seeding new orbits. Comparison between the exact and approximated manifold resulted in errors well below the determined baseline, even when the shape functions of the interpolant curve were chosen to be the simplest polynomial function with independent end-point derivatives. This leaves space for further improvement if more elaborate shape functions are used, or even parametric optimization techniques, with the cautions mentioned in Sec. VI B.

To illustrate the different manifold tracing routines, we obtained the invariants for the Chirikov-Taylor map during the methods comparison, and used the approximation technique to study the chaotic transitions between local and global motions in the conservative Duffing oscillator, where the transition from homoclinic to heteroclinic connections between the central and helical saddles in the chaotic layer was evidenced in agreement with the observed intermittent dynamics.

This work was developed with financial support from the National Council for Scientific and Technological Development (CNPq, Brazil) Grant No. 433671/2016-5, the São Paulo Research Foundation (FAPESP, Brazil) Grant Nos. 2012/18073-1 and 2011/19296-1, and the U.S. Department of Energy under DE-FC02-04ER54698 and DE-SC0012706.

In Sec. III, we described with some detail how to determine the normal displacement function coefficients for an approximated local manifold near a saddle point of an arbitrary map. Any planar map $T:R2→R2$ consists of two functions $T1,2$⁠, one for each dimension

The unstable unit vector $u^$ satisfies

with $Ti,j=∂∂xjTi(x1,x2)x∗$⁠, and its determination is elementary in two dimensions. The entries of $u^$ are related to the angle $θ$⁠, between $u^$ and $x^1$⁠, by

Coordinates $(u,v)$ are such that $u$ is the displacement along the unstable subspace with respect to the saddle $x∗$⁠, and $v$ is the displacement perpendicular to $u$⁠. These coordinates are just a rotation by $−θ$ of coordinates $(x1,x2)$ and, analogously, the representation $(Tu,Tv)$ of map $T$ can be obtained by

Using elementary transformation rules, we can obtain the required elements for the calculations in (14). For the Jacobian elements, we have

and for the Hessian elements

where $Ti,jk=∂2∂xj∂xkTi(x1,x2)x∗$ can be obtained by elementary differentiation of the original map or can be obtained numerically by finite differences when an explicit form is not available.