Painting the Phase Portrait of a Dynamical System with the Computational Tool of Lagrangian Descriptors

Introduction

“If one seeks to visualize the pattern formed by these two curves and their infinite number of intersections, each corresponding to a doubly asymptotic solution, these intersections form a kind of lattice-work, a weave, a chain-link network of infinitely fine mesh; each of the two curves can never cross itself, but it must fold back on itself in a very complicated way so as to recross all the chain-links an infinite number of times. One will be struck by the complexity of this figure, which I am not even attempting to draw. Nothing can give us a better idea of the intricacy of the three-body problem, and of all the problems of dynamics in general, …”

—H. Poincaré – New Methods of Celestial Mechanics (1899)

This prescient quote of Poincaré describing a fundamental structure of chaotic dynamics is at the heart of both the importance of invariant geometrical structures in phase space and the need for capabilities to visualize such structures. We will describe a technique, the method of Lagrangian descriptors, that offers sufficient flexibility and generality to apply to a wide range of applications in mathematics, science, and engineering.

As a preliminary motivation we describe the well-known and studied forced Duffing equation (or, a particle moving in a two-well potential subject to a time-periodic excitation) Wig03:

$$\begin{equation} \begin{cases} \dfrac{dx}{dt} = y \\[3.0pt] \dfrac{dy}{dt} = x - x^3 + A \cos (\omega t) \end{cases} \tag{1}{} \end{equation}$$

where $A$ is the amplitude of the time-dependent excitation and $\omega$ its angular frequency. In Fig. 1 we reveal the phase portrait of this dynamical system at time $t = 0$, using a driving force characterized by an amplitude $A = 0.3$ and an angular frequency $\omega = \pi$. For this calculation, we have applied the method of Lagrangian descriptors with an integration time $\tau = 13$ over a high resolution grid of initial conditions of size $1000\times 1000$. The Lagrangian descriptor field is shown in the left-hand panel from which the stable and unstable manifolds of the hyperbolic trajectory near the origin can be extracted, which are shown in the right-hand panel. The method requires only simple mathematical procedures that can be carried out by beginning graduate students, but the possible applications are extensive, as we will describe.

Figure 1.

(a) Lagrangian descriptor scalar field at time $t = 0$ for the Duffing system in Eq. 1. The integration time used is $\tau = 13$, the amplitude of the forcing is $A = 0.3$, and the angular frequency is selected as $\omega = \pi$. (b) Intricate geometry of the homoclinic tangle formed by the stable (blue) and unstable (red) invariant manifolds of the hyperbolic trajectory for the Duffing oscillator. These structures have been extracted from the Lagrangian descriptor values depicted in the left panel.

Computational tools are often essential for extracting quantitative predictions from complicated nonlinear dynamical systems. Insights from computations may also develop intuition and point the way to the discovery of new theorems. However, nonlinear dynamical systems have characteristics that make efforts to carry out geometrical analyses in phase space through computational explorations ambiguous and impractical. In particular, high dimensionality and strong instabilities can make the computation of geometric structures, such as normally hyperbolic invariant manifolds and their stable and unstable manifolds or Kolmogorov-Arnold-Moser tori, computationally intractable. Even if such structures could be computed, how one might visualize the structures in a manner that yields dynamical insights is unclear. High dimensionality and instability are not the only issues. Nonlinear dynamical systems may also have complicated, even stochastic, time dependence that complicates efforts to reveal time-varying phase space structures that organize the dynamics. Finally, the constantly growing availability of data sets describing a variety of phenomena across a wide range of time and length scales is resulting in dynamical models formulated as data sets (rather than the traditional equation based models) which, in addition to the issues described above, gives rise to a host of new challenges. After giving a description of the method of Lagrangian descriptors we will show how it meets these challenges in the setting of three examples.

(1)

Detection of dynamical matching. This is a dynamical phenomenon discovered in the context of reactions of organic molecules. Its understanding involves the simultaneous interaction of multiple phase space structures, i.e., unstable periodic orbits and their stable and unstable manifolds, and is ideal for study by Lagrangian descriptors since the method describes the interaction of all phase space structures in a selected region of phase space in a single calculation.

(2)

The Lagrangian descriptor as a stochastic process. The method of Lagrangian descriptors applies to systems with an arbitrary time dependence, which includes stochastic time dependence. The geometry of phase space structures and invariant manifolds in stochastic systems has not received extensive study in applications and we will apply this method to a familiar system: the Duffing oscillator with stochastic excitation.

(3)

Optimal path planning for ocean gliders. This is an example where the method is applied to a vector field defined as a data set. The data set describes the flow field in the North Atlantic ocean and is produced by the Copernicus Marine Environmental Monitoring Service (CMEMS). It was used to plan the journey of a trans-Atlantic glider in a way that took advantage of Lagrangian structures in the flow in order to optimize the speed and direction of the glider.

The Method of Lagrangian Descriptors

The method of Lagrangian descriptors is a mathematical and computational technique of dynamical systems theory that was originally developed in the context of geophysical fluid dynamics to analyze Lagrangian transport and mixing processes by means of identifying hyperbolic trajectories and their stable and unstable manifolds in flows with a general time dependence JMM09 (hence the word “Lagrangian” in its name).

Such invariant manifolds having dimension less than the ambient phase space act as transport barriers and therefore partition the phase space of the system into regions where trajectories have distinct dynamical behavior. This provides us with a geometrical template (or skeleton) that qualitatively describes the solutions of the differential equations governing the underlying dynamics. The idea of not looking for specific solutions through a given initial condition but instead seeking geometrical relationships between all solutions is the qualitative, geometrical approach to the study of dynamical systems theory developed by Henri Poincaré in the course of his studies of the three body problem of celestial mechanics BG97 and has had a profound impact on our study and understanding of the nonlinear character of natural phenomena ever since. It is this approach that is significantly facilitated by the method of Lagrangian descriptors, as we will show. Before showcasing several applications we want to establish the setting and provide a simple, yet general, description of the method.

Suppose that we are interested in obtaining a global picture of the phase portrait associated to the dynamical system described by the differential equation:

$$\begin{equation} \dfrac{d \mathbf{x}}{dt} = \mathbf{f}(\mathbf{x}(t),t) \,,\;\; \mathbf{x} \in \mathcal{M} \;,\; t \in \mathbb{R} \,, \tag{2}{} \end{equation}$$

where $\mathcal{M}$ is an $n$ dimensional phase space (set of all possible states), $\mathbf{x}(t)$ represents the state of the system at time $t$, and $\mathbf{f} = (f_1,\ldots ,f_n)$ corresponds to the vector field determining the time evolution of trajectories. We highlight the fact here that there are no restrictions on the time dependence of Eq. 2. It may be continuous, discrete, stochastic, or even time independent. Examples of each of these can be found in ASGGK$^{+}$21AASGG$^{+}$20. In applications, the vector field that defines the system, such as those obtained by remote sensing techniques, may take two forms: the function $\mathbf{f}$ may be represented by known analytical functions. On the other hand, the vector field may be specified by a discrete spatio-temporal data set acquired from the numerical simulation of a mathematical model or physical measurements.

Once the problem is defined in terms of a system of nonlinear ordinary differential equations, the goal is to understand the global dynamics in phase space in terms of the qualitatively distinct trajectories. This will require extracting the relevant data from large collections of trajectories. Lagrangian descriptors provide a tool for achieving this goal as it provides a simple and effective approach for classifying distinct dynamical behavior that focuses on the initial conditions of trajectories, rather than their futures and pasts after a specified amount of time. This implies that phase space structures are encoded in the initial conditions of the trajectories, which constitute the “primitive objects” used to explore the underlying dynamics. It is this property that provides an approach for analyzing high dimensional systems, since it allows for the computation of high resolution and low dimensional slices of the phase space, and it is the high resolution aspect that may reveal geometrical structures in high dimensional phase space. Specific examples can be found in AASGG$^{+}$20KGGW21ASGGK$^{+}$21, which we touch on in the conclusions.

The idea at the heart of Lagrangian descriptors is very simple and relies on assigning a scalar value to each initial condition in the phase space of the system. This trajectory-based diagnostic tool produces a scalar field that is calculated as follows. Consider a region of the phase space where one would like to reveal structure at time $t = t_0$ and define on it a uniformly spaced grid of initial conditions $\mathbf{x}_0 = \mathbf{x}(t_0)$. Once this setup is in place, the next step is to integrate all the initial conditions on the grid which will accumulate along their trajectories the values attained by a non-negative scalar function, $\mathcal{F}(\mathbf{x},t)$, which depends on the phase space variables. Trajectories are propagated both forward and backward for some time interval $[t_0-\tau ,t_0+\tau ]$, where $\tau$ defines the integration time (see Fig. 2). This computation yields a scalar field which highlights the location of the invariant stable and unstable manifolds. Mathematically, the method of Lagrangian descriptors is defined as follows:

$$\begin{equation} \mathcal{L}\left(\mathbf{x}_0,t_0,\tau \right) = \mathcal{L}_{f}\left(\mathbf{x}_0,t_0,\tau \right) + \mathcal{L}_{b}\left(\mathbf{x}_0,t_0,\tau \right) \,, \tag{3} \end{equation}$$

where the forward $(\mathcal{L}_{f})$ and backward $(\mathcal{L}_{b})$ components of the Lagrangian descriptor function have the expressions:

$$\begin{equation} \begin{split} \mathcal{L}_{f}\left(\mathbf{x}_0,t_0,\tau \right) &= \int _{t_0}^{t_0+\tau } \mathcal{F}(\mathbf{x}(t;\mathbf{x}_0),t) \, dt \,, \\[3.0pt] \mathcal{L}_{b}\left(\mathbf{x}_0,t_0,\tau \right) &= \int _{t_0-\tau }^{t_0} \mathcal{F}(\mathbf{x}(t;\mathbf{x}_0),t) \, dt \,. \end{split} \tag{4}{} \end{equation}$$

The simple formulation of Lagrangian descriptors highlights why it is so appealing, since its application only requires elementary skills such as the propagation of classical trajectories. Therefore, this technique is readily accessible to all practitioners of dynamical systems, from beginners to experts.

As a mathematical object, the Lagrangian descriptor is a scalar field of the form:

$$\begin{equation} \begin{array}{ccccc} \mathcal{L} & : & \mathcal{D} \times \mathcal{I} & \longrightarrow & [0,+\infty ) \\[3.0pt] & & (\mathbf{x}_0, t_0) & \longmapsto & {\cal L} (\mathbf{x}_0, t_0, \tau ) \end{array} \tag{5} \end{equation}$$

where ${\cal D} \subset {\cal M}$ is a chosen subset of the phase space $\cal M$, ${\cal I} \subset \mathbb{R}$ is a chosen subset of the time domain, and $\tau$ is regarded as fixed. The choice of $\cal D$ is made that is appropriate for a particular application, as we will see in the examples, e.g., it could be the entire phase space, a region of the phase space whose dynamical significance is relevant for the application, or a lower dimensional “slice” of the phase space which may be appropriate for analyzing high dimensional problems.

Figure 2.

Initial condition $\mathbf{x}_0 = \mathbf{x}(t_0)$ at time $t = t_0$ propagated forward and backward in time for a given integration time $\tau$. A value for the Lagrangian descriptor function is assigned to the phase space point $\mathbf{x}_0$ by accumulating a non-negative function along the trajectory that results from its evolution.

The implementation of Lagrangian descriptors is straightforward since, thanks to the Fundamental Theorem of Calculus, one can directly incorporate the computation of the integrals into any integrator used to propagate trajectories of the dynamical system. Indeed, we just need to add one more component to the vector field accounting for the integrand $\mathcal{F}$ that defines the Lagrangian descriptor function:

$$\begin{equation} \begin{cases} \dfrac{d \mathcal{L}_{f}}{d t} = \mathcal{F}(\mathbf{x}(t_0 + t;\mathbf{x}_0),t_0+t) \\[3.0pt] \dfrac{d \mathcal{L}_{b}}{d t} = \mathcal{F}(\mathbf{x}(t_0 - t;\mathbf{x}_0),t_0-t) \end{cases} ,\;\; t \in [0,\tau ]. \tag{6} \end{equation}$$

Therefore, Lagrangian descriptors can be calculated on the fly simultaneously as trajectories are propagated forward and backward in time, without becoming a computational burden. The scalar field obtained from the method, which synthesizes all the dynamical information concealed in the evolution of an ensemble of trajectories starting from a given set of initial conditions, has the capability of constructing the phase portrait of the underlying dynamical system. A detailed description of the algorithm for computing Lagrangian descriptors, along with software and examples, can be found in AASGG$^{+}$20.

In the literature, several definitions of LDs have been introduced depending on the non-negative scalar-valued function $\mathcal{F}$ used to define the method. In its origins JMM09MWCM13, the arclength of trajectories in phase space was applied, that is, $\mathcal{F}(\mathbf{x},t) = \left\lVert \mathbf{f}(\mathbf{x},t) \right\rVert$, and this provided a clear geometrical interpretation of the output generated by Lagrangian descriptors. However, other diagnostics have been widely used such as the action integral for Hamiltonian systems AASGG$^{+}$20 or the function:

$$\begin{equation} \mathcal{F}(\mathbf{x},t) = \sum _{i=1}^n |f_{i}(\mathbf{x},t)|^p \;,\;\; p \in (0,1] \tag{7}{} \end{equation}$$

inspired by the $p$-norm that arises in functional analysis LBIGG$^{+}$17. All these different types of Lagrangian descriptors are capable of highlighting the stable and unstable manifolds at locations where the scalar field values exhibit an abrupt change. These characteristics in the field, which align with the manifolds, are referred to as “singular features” and they partition the phase space into regions where trajectories have distinct dynamical behaviors. In particular, the computation of the forward component of the Lagrangian descriptor $\mathcal{L}_{f}$ detects the locations of the stable manifolds, whereas the backward Lagrangian descriptor, $\mathcal{L}_{b}$, does the same but in this case for the unstable manifolds.

The sharp transitions obtained in the Lagrangian descriptor scalar values across the stable and unstable manifolds produce large changes in its gradient that are very easy to detect when plotting the output provided by the method. In LBIGG$^{+}$17, a rigorous mathematical connection between the “singular features” displayed by the Lagrangian descriptor output and the stable and unstable manifolds was established for example systems. This was carried out for the Lagrangian descriptor defined from Eq. 7. With this approach, the Lagrangian descriptor scalar field becomes in fact non-differentiable at the phase space points belonging to a stable or unstable manifold. This property is crucial in many ways since it allows us to easily extract the location of the stable and unstable manifolds in the Lagrangian descriptor plot by applying edge and ridge detection algorithms such as those widely used in image processing. In addition to invariant manifolds, the method also provides information about other phase space structures, such as Kolmogorov-Arnold-Moser tori, but in this case they appear in regions where the Lagrangian descriptor function varies smoothly. In fact, tori structures can be visualized by means of computing long-term time averages of Lagrangian descriptors as discussed in LBIGG$^{+}$17.

To complete this section, we demonstrate how the method of Lagrangian descriptors is applied to analyze the phase portrait of the following two dimensional, autonomous dynamical system where the stable and unstable manifolds of a hyperbolic equilibrium point can be computed analytically:

$$\begin{equation} \begin{cases} \dot{x} = \dfrac{\partial \mathcal{H}}{\partial y} = y \\[3.0pt] \dot{y} = - \dfrac{\partial \mathcal{H}}{\partial x} = \dfrac{x}{(1+x^2)^2} \end{cases} . \tag{8}{} \end{equation}$$

This system is described by the Hamiltonian function:

$$\begin{equation} \mathcal{H}(x,y) = \dfrac{y^2}{2} + \dfrac{1}{2(1+x^2)} \,, \tag{9} \end{equation}$$

and it is a simple exercise to show that this system has a hyperbolic equilibrium point at the origin, and that the equations of its stable ($\mathcal{W}^s$) and unstable ($\mathcal{W}^u$) manifolds are given by:

$$\begin{equation} \begin{split} \mathcal{W}^s(0,0) &= \left\lbrace (x,y) \in \mathbb{R}^2 \; : \; y = \dfrac{-x}{\sqrt {1+x^2}} \right\rbrace \\[3.0pt] \mathcal{W}^u(0,0) &= \left\lbrace (x,y) \in \mathbb{R}^2 \; : \; y = \dfrac{x}{\sqrt {1+x^2}} \right\rbrace . \end{split} \tag{10} \end{equation}$$

Consider a uniformly spaced grid of $500\times 500$ initial conditions in the rectangle $[-2,2]\times [-1.25,1.25]$ and integrate them forward and backward, starting at $t_0 = 0$, for a time $\tau = 3$. Along their evolution, we accumulate the phase space arclength of the trajectories, and display the total Lagrangian descriptor scalar field obtained from this computation, that is $\mathcal{L} = \mathcal{L}_f+\mathcal{L}_b$, in Fig. 3. Notice that the value attained by Lagrangian descriptors at the origin is zero, since it is an equilibrium point and hence it is not moving. Moreover, if we depict the forward and backward components of Lagrangian descriptors separately, they will reveal the locations of the stable and unstable manifolds respectively. We illustrate this property in Fig. 4 where we have depicted only the forward LD, given by $\mathcal{L}_f$.

Figure 3.

Total Lagrangian descriptor ($\mathcal{L}_f+\mathcal{L}_b$) calculated for the dynamical system in Eq. 8 with an integration time $\tau = 3$. The Lagrangian descriptor used to reveal the underlying geometry is that defined from the arclength of trajectories in phase space. The output has been normalized by subtracting the minimum value of the Lagrangian descriptor function, and then dividing by the difference between its maximum and minimum values.

Figure 4.

Forward Lagrangian descriptor function $(\mathcal{L}_f)$ calculated for the dynamical system in Eq. 8 with an integration time $\tau = 3$. The Lagrangian descriptor definition used to visualize the phase portrait is that of the arclength of trajectories in phase space.

However, although the manifolds are visible in Figs. 3 and 4, they are poorly resolved. The reason this occurs is that the integration time used for this system is small and hence the output provided by Lagrangian descriptors cannot distinguish clearly between distinct dynamical behaviors of the trajectories. The detection of phase space structure as revealed by the method will improve when the integration time gets larger, since this implies incorporating more information about the dynamical history of the trajectories. If we repeat the calculation of Lagrangian descriptors, but using an integration time $\tau = 10$, the manifolds appear as sharp features in the scalar field, see Fig. 5(a). This becomes evident when we look at the values attained by the Lagrangian descriptor along a line in the phase space, say $y = 0.5$, and plot them as we have done in Fig. 5(b). The gradient exhibits an abrupt change in its values at the points where the manifolds are located. Furthermore, the Lagrangian descriptor displays a local minimum at the manifolds. This characteristic behavior of the method on the manifolds can be used to extract these objects from the scalar field by means of applying a ridge/edge detection algorithm, such as those developed in image processing. We have done so in Fig. 5(c), calculating the norm of the gradient of the scalar field, $\left\lVert \nabla \mathcal{L}\right\rVert$, and then discarding those points below a certain threshold. This operation is carried out separately for the forward and backward components of the Lagrangian descriptor, and the stable and unstable manifolds are drawn in blue and red, respectively.

Figure 5.

(a) Total Lagrangian descriptor function (forward$+$backward) calculated for the dynamical system in Eq. 8 with an integration time $\tau = 10$. (b) Value of Lagrangian descriptors along the line $y = 0.5$ depicted in magenta in the first panel. (c) Geometry of the stable (blue) and unstable (red) manifolds extracted from the sharp ridges (singular features) of the Lagrangian descriptor scalar field displayed in the left panel.

But why does the method reveal the manifolds? To provide an intuitive explanation to this question, we focus for example on the forward component of Lagrangian descriptors, $\mathcal{L}_f$. Take two neighboring initial conditions, one lying on the curve that corresponds to the stable manifold and another slightly off it. If we integrate them forward for a small time $\tau$, the segments of both trajectories are comparable in length, so that the value of Lagrangian descriptors obtained for both initial conditions is almost equal. However, if the integration time $\tau$ is larger, the arclengths of the two trajectories become very different. Indeed, since the initial condition on the stable manifold converges to the origin as $t \to \infty$, the value attained by $\mathcal{L}_f$ at that point is bounded, since the curve segment of the stable manifold between this point and the origin has a finite arclength. On the other hand, the value of the Lagrangian descriptor function at the initial condition that does not belong to the manifold can be made as large as we desire if we set a very large integration time, because the trajectory does not converge to the equilibrium point at the origin. Moreover, if we consider a curve of initial conditions that crosses transversely the stable manifold, the Lagrangian descriptor value along it attains a minimum on the manifold. Notice also that, by following this same argument, but constructing the backward time evolution term of Lagrangian descriptors, $\mathcal{L}_b$, we can arrive at the conclusion that backward integration of initial conditions will highlight the unstable manifold of the hyperbolic equilibrium point at the origin. This heuristic argument justifies why the Lagrangian descriptor values vary significantly near the invariant manifolds in comparison to those elsewhere, and as a consequence the method successfully reveals these mathematical structures in the resulting scalar field.

Detection of Dynamical Matching

In this section we describe a phase space dynamical phenomenon that occurs in a variety of reactions of organic molecules. This phenomenon was discovered by Carpenter Car95 and is referred to as dynamical matching. The existence of dynamical matching has been extensively confirmed both through molecular dynamics simulations carried out for real chemical systems, and also by means of simple mathematical models built with the purpose of developing a deeper theoretical understanding of this phenomenon and shed some light on the reasons why it occurs. As we will describe, its understanding involves the interaction of several phase space structures simultaneously, which makes its study particularly appropriate for the method of Lagrangian descriptors.

The prototypical setting for the study of dynamical matching is carried out in the context of a two degree-of-freedom Hamiltonian system defined by a caldera potential energy surface (PES). The name derives from the fact that the topography of the PES resembles the depression that is formed in a terrain when a volcano erupts and then collapses. An analytical form for a caldera PES is derived in CKC$^{+}$14 and a graphical illustration of this PES is shown in Fig. 6.

Figure 6.

An example of a caldera-type potential energy surface.

We see from the figure that the caldera PES is characterized by a local minimum in a relatively flat region, surrounded by four saddle points that control access/exit to/from the central area where the minimum is located. Typically, two of the saddles have low energy values, while the other two are higher in energy.

Computation of trajectories of Hamilton’s equations defined by the caldera PES for different total energies CKC$^{+}$14KW19 have revealed two main types of dynamical behavior for trajectories in such systems: dynamical matching and trajectory trapping (transitory or permanent) in the central region of the caldera. The phenomenon of dynamical matching can be understood as the tendency to conserve the direction of momentum that some trajectories display along their evolution. What this means is that for a caldera PES with the topographical features we have described above, it is expected that the directionality associated with some of the trajectories initialized at a high energy saddle point of the system is preserved. The consequence of this dynamical feature is that trajectories entering the caldera by passing over the upper right-hand saddle, exit by passing over the lower left-hand saddle. “Dynamical matching” means that entering trajectories with this dynamical property are “matched” with a particular manner of exiting. This gives a strategy for predicting trajectory behavior. The goal here is to understand the phase space structures that “guide” this trajectory behavior.

We will show how Lagrangian descriptors can be used to reveal the phase space structures responsible for dynamical matching by considering the two degree-of-freedom Hamiltonian system developed in CKC$^{+}$14KW19, where the Hamiltonian to explore dynamical matching is expressed as the sum of the kinetic and potential energies:

$$\begin{equation} \mathcal{H}(x,y,p_x,p_y) = \dfrac{p_x^2 + p_y^2}{2} + V(x,y) \,. \tag{11} \end{equation}$$

In these works, the caldera-type PES is described by the function:

$$\begin{equation} V(x,y) = c_1 \, ((\lambda x)^2+y^2) + c_2 \, y - c_3 \, ((\lambda x)^4 + y^4 - 6 (\lambda x)^2 y^2) \,, \tag{12} \end{equation}$$

which depends on the parameters $c_1,c_2,c_3$, and $\lambda \in (0,1]$. The latter parameter allows one to probe the influence of symmetry of the PES by stretching it in the $x$-direction. Hamilton’s equations of motion have the form:

$$\begin{equation} \begin{cases} \dot{x} = \dfrac{\partial \mathcal{H}}{\partial p_x} = p_x \\[3.0pt] \dot{y} = \dfrac{\partial \mathcal{H}}{\partial p_y} = p_y \\[3.0pt] \dot{p}_x = -\dfrac{\partial \mathcal{H}}{\partial x} = 2 \lambda ^2 x \left[2 \, c_3 \left(\lambda ^2x^2 - 3y^2\right) - c_1\right] \\[3.0pt] \dot{p}_y = -\dfrac{\partial \mathcal{H}}{\partial y} = 2 \, y \left[2 \, c_3\left(y^2-3\lambda ^2x^2\right) - c_1\right] - c_2. \end{cases} \tag{13}{} \end{equation}$$

This dynamical system has five equilibrium points: a minimum close to the origin of center type stability and four saddle points with saddle$\times$center stability that are related by a reflection symmetry about the $y$ axis. Moreover, the saddles located in the upper-half of the OXY plane ($y > 0$) have higher energy than those in the lower-half plane ($y < 0$).

In CKC$^{+}$14KW19, a parametric study of the dynamics of this system was conducted in terms of the total energy and the stretching coefficient, and simulations were performed using the model parameter values $c_1 = 5$, $c_2 = 3$ and $c_3 = -3/10$. These are the same parameters used, as well as $\lambda =1$, for Fig. 6.

Ensembles of trajectories were initialized near the upper right-hand saddle for energy $\mathcal{H} = 29$ and with initial momentum directed towards the caldera. It was observed that for the unstretched caldera PES ($\lambda = 1$) all the trajectories exhibited dynamical matching. Hence, they evolved across the central region of the caldera, preserving the direction of their momentum vector, and escape the PES through the diagonally opposite exit channel associated to the lower left-hand saddle. Interestingly, as the stretching parameter $\lambda$ was varied, it was found that there exists a critical value below which full dynamical matching behavior is destroyed and some trajectories became temporarily trapped in the central area of the caldera, before they exited.

While the saddle points certainly have an influence on the dynamics, they do not exist on the energy surface $\mathcal{H} = 29$. Instead, the Lyapunov subcenter theorem states that saddle type periodic orbits exist for energies above that of the saddles Wig03. From the point of view of dynamical systems theory it is therefore natural to ask if heteroclinic connections between these periodic orbits are responsible for guiding trajectories in the dynamical matching mechanism. This was recently shown to be the case in KGGW20. In order to describe this in more detail, we recall that the stable and unstable manifolds of unstable periodic orbits of two degree-of-freedom Hamiltonian systems have the topology of cylinders, i.e. $S^{1}\times \mathbb{R}$. Moreover, when they intersect a two dimensional section in the three dimensional energy surface they give rise to closed curves which have been referred to as “reactive islands” OdADLMM90. We give a geometrical representation of this situation in Fig. 7.

Figure 7.

Stable or unstable manifold of a 2 DoF Hamiltonian system intersecting a transverse surface of section, giving rise to a closed curve known as a reactive island.

In the case of the unstretched PES ($\lambda = 1$), where all trajectories initialized near the upper right-hand unstable periodic orbit exhibit dynamical matching behavior, it has been shown that this occurs because the unstable manifold of the upper-right unstable periodic orbit is contained inside the region bounded by the stable manifold of the lower-left unstable periodic orbit. Therefore, we would have the geometrical situation where one cylinder lies inside the other. However, when the PES is stretched there exists a critical value of $\lambda$ (depending on the energy) for which a heteroclinic intersection between these manifolds occurs and dynamical matching begins to break down, giving rise to trajectories that get diverted and become temporarily trapped in the central region of the caldera.

We will show next that the method of Lagrangian descriptors can be applied to detect the existence of dynamical matching in the system described by Eq. 13, and provide an estimation for the critical value of the stretching parameter at which the heteroclinic intersection develops. To do so, we calculate Lagrangian descriptors on the following two dimensional section:

$$\begin{equation} \Sigma = \left\lbrace (x,y,p_x,p_y) \in \mathbb{R}^4 \; : \; y = 0 \;,\; p_y < 0 \right\rbrace , \tag{14}{} \end{equation}$$

where the total energy of the system is fixed at $\mathcal{H} = 29$. In Fig. 8 we display the output of Lagrangian descriptors calculated for the unstretched system ($\lambda = 1$) with an integration time $\tau = 4$. From this plot we can clearly see that the unstable manifold of the upper-right unstable periodic orbit is inside the stable manifold of the lower-left unstable periodic orbit, and this geometrical configuration gives rise to full dynamical matching.

Figure 8.

Lagrangian descriptor scalar field (using the $p$-norm definition with $p = 1/2$) calculated on the two dimensional section in Eq. 14 for the dynamical system in Eq. 13 with an integration time $\tau = 4$. The total energy is $\mathcal{H} = 29$ and the stretching parameter is chosen as $\lambda = 1$ (unstretched PES). UPO stands for unstable periodic orbit.

On the other hand, if we set $\lambda = 0.72$, Lagrangian descriptors reveal, as displayed in Fig. 9(a), the intersection between the manifolds responsible for reducing the effect of dynamical matching in the system. This produces transitory trapping of the trajectories located inside the lobe formed by the connection between the invariant manifolds. Hence, Lagrangian descriptors allow us to easily identify the location of these special trajectories, and distinguish between distinct dynamical behaviors. We illustrate this fact by selecting two initial conditions marked as green and cyan dots in Fig. 9(a) and show their evolution in forward time, projected onto configuration space, in Fig. 9(b). While the green initial condition displays dynamical matching, the cyan one is partially trapped in the central area of the caldera.

Figure 9.

(a) Lagrangian descriptor scalar field (using the $p$-norm definition with $p = 1/2$) calculated on the two dimensional section in Eq. 14 for the dynamical system in Eq. 13 with $\tau = 4$. The total energy is set to $\mathcal{H} = 29$ and the stretching parameter is $\lambda = 0.72$. (b) Projection onto the OXY plane of the trajectories of two initial conditions selected in the left panel. Evolution in forward time is depicted in blue and red whereas backward propagation is represented in black. The green and cyan initial conditions display dynamical matching and partial trapping behavior, respectively.

We end this section by demonstrating how Lagrangian descriptors reveal the intricate geometry of the stable and unstable manifolds with an ever increasing level of detail, as the integration time $\tau$ used for its computation gets larger. This is expected because the propagation of trajectories for longer times provides more information about their dynamical fates, leaving a distinctive fingerprint on the Lagrangian descriptor scalar field values that allows for a clear detection of the invariant manifolds. We illustrate this property in Fig. 10, where the locations of the stable and unstable manifolds of all the existing UPOs in the system have been calculated with Lagrangian descriptors using $\tau = 15$. Therefore, Lagrangian descriptors identify all these dynamical objects in the same calculation without the need to apply other classical methods from nonlinear dynamical systems that rely on locating first the periodic orbits, and then globalizing the invariant manifolds from them by numerically propagating trajectories along the appropriate eigendirections.

Figure 10.

Stable (blue) and unstable (red) manifolds extracted from the output of Lagrangian descriptors calculated on the two dimensional section defined in Eq. 14 for the system in Eq. 13 with $\tau = 15$. The total energy is set to $\mathcal{H} = 29$ and the stretching parameter is $\lambda = 0.72$.

The Lagrangian Descriptor as a Stochastic Process

Dynamical systems that occur in Nature are typically subjected to the effects of random influences or noise that come in the form of external forcing, internal agitation due to thermal fluctuations, uncertainties arising from the measurements of initial conditions, or as a result of our lack of knowledge on the values attained by the parameters describing the mathematical model under study, etc. Probably, the most classical example of a stochastic dynamical system is that of Brownian motion, which is characterized by the Langevin equation modeling the action produced by a solvent on the random motion of colloidal particles (such as large molecules or proteins). Randomness is thus an important ingredient to take into account when modeling complex phenomena that occur in chemical, biological, physical, and other type of systems that evolve in time, and this is done through the use of stochastic differential equations (SDEs) Dua15.

We show here how to construct a stochastic Lagrangian descriptors that can reveal the phase space structure of a continuous random dynamical system (RDS) with general time dependence described by a set of stochastic differential equations in the form:

$$\begin{equation} d \mathbf{X}_t = \mathbf{b} \left(\mathbf{X}_t,t\right) \, dt + \boldsymbol{\sigma }(\mathbf{X}_t,t) \, d\boldsymbol{\xi }_t \;, \tag{15}{} \end{equation}$$

where $\mathbf{b} \in \mathbb{R}^n$ is the deterministic vector field, known as the drift, while $\boldsymbol{\sigma } \in \mathbb{R}^{n \times m}$ is called the random forcing intensity. The stochastic process $\mathbf{X}_t \in \mathbb{R}^n$ is the solution of the SDE in Eq. 15 satisfying the initial condition $\mathbf{X}_{t_0} = \mathbf{x}_0$ at time $t = t_0$, and $\boldsymbol{\xi }_t \in \mathbb{R}^m$ is a stochastic process defined by a Brownian motion (a Gaussian process) or Lévy motion (a non-Gaussian process). The noise term $\boldsymbol{\sigma } d\boldsymbol{\xi }_t$ provides us with a measure of the uncertainty acting on the system, which can be caused by the fluctuations in the parameters of the mathematical model, or the influence of external random factors.

Consider a probability space $(\Omega ,\mathcal{F},\mathbb{P})$, where $\Omega$ denotes the sample space, $\mathcal{F}$ is a $\sigma$-field and $\mathbb{P}$ is the probability measure. If we take a random experiment $\omega \in \Omega$, this gives a particular realization of the stochastic processes involved in the RDS, and hence the solution to Eq. 15 can be formally written in closed-form as:

$$\begin{equation} \mathbf{X}_t(\omega ) = \mathbf{x}_0(\omega ) + \int _{t_0}^t \mathbf{b} \left(\mathbf{X}_s(\omega ),s\right) \, ds + \int _{t_0}^t \boldsymbol{\sigma }(\mathbf{X}_s(\omega ),s) \, d\boldsymbol{\xi }_s(\omega ) \;. \tag{16}{} \end{equation}$$

Notice that in this expression, the second integral is a stochastic integral, and thus it has to be interpreted in terms of It$\hat{o}$ or Stratonovich calculus. We will focus here on the case where $\boldsymbol{\xi }_t = \mathbf{W}_t \in \mathbb{R}^m$ is a Wiener process, in particular a two-sided Brownian motion. The sample space in this context is termed a Wiener space and is defined by $\Omega = C_0(\mathbb{R};\mathbb{R}^m)$, that is, the class of continuous functions that vanish at the origin. Under such assumptions, the existence and uniqueness of solutions for Eq. 15 requires that the functions $\mathbf{b}$ and $\boldsymbol{\sigma }$ are locally Lipschitz.

In order to introduce a Lagrangian descriptor to reveal the phase space of the RDS in Eq. 15, it is important to allow for forward and backward time evolution, and therefore, we need to use a two-sided Brownian motion $\boldsymbol{W}_t = \left(W^{1}_t,\ldots ,W^{m}_t\right) \in \mathbb{R}^m$. This stochastic process satisfies the following properties:

(1)

$\mathbf{W}_0 = \mathbf{0}$ almost surely.

(2)

$\mathbf{W}_t$ has independent increments, that is, if $0 < t_1 < t_2 < \ldots < t_{n-1} < t_n$, then the random variables $\mathbf{W}_{t_1}$, $\mathbf{W}_{t_2} - \mathbf{W}_{t_1},\ldots ,\mathbf{W}_{t_n} - \mathbf{W}_{t_{n-1}}$ are independent.

(3)

For any $t,s \in \mathbb{R}$, the increment $\mathbf{W}_{t} - \mathbf{W}_{s}$ follows a normal distribution with zero mean and variance $|t-s| \mathbb{I}_m$, that is, $\mathbf{W}_{t} - \mathbf{W}_{s} \sim \mathcal{N}(0,|t-s| \mathbb{I}_m)$, where $\mathbb{I}_m$ is the $m \times m$ identity matrix.

To construct the solution $\mathbf{X}_t = \left(X^{1}_t,\ldots ,X^{n}_t\right) \in \mathbb{R}^n$ of the RDS in Eq. 15 for a time interval $\mathbb{T} = [t_0-\tau ,t_0+\tau ]$, we start by selecting an initial condition $\mathbf{X}_{t_0} = \mathbf{x}_0 \in \mathbb{R}^n$ at time $t = t_0$. The next step is to apply any numerical method to solve an SDE, such as for example the Euler-Maruyama (EM) scheme. With the purpose of integrating forward in time, define a regular partition of $\mathbb{T}^+ = [t_0,t_0 + \tau ]$ into $N$ subintervals of sufficiently small length $\Delta t = \tau /N > 0$. This gives rise to the ordered sequence of times $t_k = t_0 + k \Delta t$ with $k \in \lbrace 0,\ldots ,N \rbrace$. Applying EM in forward time, we get that the $i$-th component of the solution evolves in discrete steps as follows:

$$\begin{equation} X^{i}_{t_{k+1}}(\omega ) = X^{i}_{t_{k}}(\omega ) + b^i\left(\mathbf{X}_{t_k}(\omega ),t_k\right) \Delta t + \sum _{j = 1}^{m} \sigma ^{i}_{j}\left(\mathbf{X}_{t_k}(\omega ),t_k\right) \Delta W^{j}_{t_k}(\omega ) \,, \tag{17}{} \end{equation}$$

where $k \in \lbrace 0,\ldots ,N-1\rbrace$ and $\omega \in \Omega$ is a random experiment that determines a realization of the stochastic processes $\mathbf{X}_t$ and $\mathbf{W}_t$. The increments of the Wiener process follow the Gaussian distribution $\Delta W^{j}_{t_k} = W^{j}_{t_{k+1}} - W^{j}_{t_{k}} \sim \mathcal{N}(0,\Delta t)$ for each $j \in \lbrace 1,\ldots ,m\rbrace$. Similarly, this process can also be carried out in backward time. To do so, we partition the interval $\mathbb{T}^- = [t_0-\tau ,t_0]$ into $N$ subintervals of equal length $\Delta t = \tau /N > 0$, producing the sequence of times $t_k = t_0 + k \Delta t$ with $k \in \lbrace -N,\ldots ,0\rbrace$. Applying EM in backward time yields in this case:

$$\begin{equation} X^{i}_{t_{k}}(\omega ) = X^{i}_{t_{k+1}}(\omega ) - b^i\left(\mathbf{X}_{t_{k+1}}(\omega ),t_{k+1}\right) \Delta t - \sum _{j = 1}^{m} \sigma ^{i}_{j}\left(\mathbf{X}_{t_{k+1}}(\omega ),t_{k+1}\right) \Delta W^{j}_{t_{k+1}}(\omega ) \,, \tag{18}{} \end{equation}$$

where $k \in \lbrace -1,\ldots ,-N\rbrace$, and the Brownian increments $\Delta W^{j}_{t_{k+1}} = W^{j}_{t_{k+1}} - W^{j}_{t_{k}} \sim \mathcal{N}(0,\Delta t)$ for each component $j \in \lbrace 1,\ldots ,m\rbrace$.

Once the solution to the SDE is known on a discrete set of points in time for a grid of initial conditions, we can define a Lagrangian descriptor function from this data set by selecting a value $p \in (0,1]$ and using the expression introduced in BILWM16, which can be written as:

$$\begin{equation} \mathcal{L}\left(\mathbf{x}_{0},t_0,\tau ,\omega \right) = \sum _{k = -N}^{-1} \sum _{i = 1}^{n} \big |X^{i}_{t_{k+1}}(\omega ) - X^{i}_{t_{k}}(\omega )\big |^p + \sum _{k = 0}^{N-1} \sum _{i = 1}^{n} \big |X^{i}_{t_{k+1}}(\omega ) - X^{i}_{t_{k}}(\omega )\big |^p \,. \tag{19}{} \end{equation}$$

It is clear that the value of this scalar field, and hence of the phase space structure revealed by it, depends fundamentally on the realization of the random experiment $\omega \in \Omega$. Therefore, the phase portrait of the RDS in Eq. 15 is unpredictable and we can only address questions about the locations of hyperbolic trajectories and their stable and unstable manifolds in terms of probabilities. If we run $M$ experiments $\mathcal{S} = \lbrace \omega _k\rbrace _{k=1}^{M} \subset \Omega$, giving rise to a collection of outcomes:

$$\begin{equation} \mathcal{L}_k = \mathcal{L}\left(\mathbf{x}_{0},t_0,\tau ,\omega _k\right) \,,\;\; k = 1,\ldots ,M \tag{20} \end{equation}$$

of the Lagrangian descriptor scalar field, one can understand the geometry of transport for an SDE by means of calculating the expectation values of the Lagrangian descriptor considered as a uniformly distributed random variable, that is:

$$\begin{equation} \left\langle \mathcal{L} \right\rangle _{\mathcal{S}} = \mathbb{E}\left[\mathcal{L}\left(\mathbf{x}_{0},t_0,\tau ,\cdot \right)\right] = \dfrac{1}{M} \sum _{k=1}^{M} \mathcal{L}_k \,. \tag{21}{} \end{equation}$$

For a sufficiently large number of experiments, this quantity will provide us with an approximate geometrical picture of the average phase portrait of the system, which is a very useful diagnostic tool for exploring the influence of random phenomena on dynamics which is discussed in more detail in Dua15.

For the purpose of illustrating how the method of Lagrangian descriptors can be used to analyze the phase space structure of a RDS, we will take a look at the dynamics of the stochastically forced Duffing oscillator, given by the system:

$$\begin{equation} \begin{cases} d X_t = Y_t \, dt \\[3.0pt] d Y_t = \left(X_t - X_t^3\right) dt + \sigma \, dW_t \end{cases} \,, \tag{22}{} \end{equation}$$

where $W_t$ is a two-sided Wiener process and $\sigma$ represents the stochastic forcing strength which we consider to be additive, that is, it does not depend on the solution. In particular, we will give it a constant value $\sigma = 0.2$. To study the phase portrait of this SDE, we apply Lagrangian descriptors at an initial time $t = 0$ to a grid of $400\times 400$ initial conditions on the rectangle $[-1.8,1.8]\times [-1.5,1.5]$, integrating trajectories forward and backward for a time $\tau = 20$. To do so, we have implemented the Euler-Maruyama scheme discussed in Eqs. 17 and 18 using a time step of $\Delta t = 0.005$. First, we carry out only one experiment and the results of this numerical simulation are shown in Fig. 11. Our next goal is to reconstruct the most likely phase portrait for the Duffing oscillator by means of performing $50$ numerical experiments and average the resulting LD scalar fields as we discussed in Eq. 21. The outcome of this calculation is displayed in Fig. 12. Notice how the method nicely highlights the homoclinic tangle formed by the stable and unstable manifolds of the RDS. These can be compared with the analogous structures in the periodically forced Duffing equation shown in Fig. 1.

Figure 11.

Phase portrait of the stochastic Duffing system in Eq. 22, as revealed by Lagrangian descriptors, for a particular random experiment.

Figure 12.

Expected phase portrait of the stochastic Duffing system in Eq. 22 for $50$ random experiments. (a) Averaged Lagrangian descriptor scalar field. (b) Approximate locations of the stable (blue) and unstable (red) invariant manifolds extracted from the values of the Lagrangian descriptor function represented in the first panel.

Optimal Path Planning for Ocean Gliders

Autonomous underwater vehicles (AUV) can travel the ocean without being controlled by an operator. An underwater glider is a particular type of AUV that generates propulsion by modulating their buoyancy at specified depths. This procedure transfers a component of the induced vertical acceleration forward by means of swept wings attached to the glider. In this section we discuss the mission of the underwater glider Silbo and describe how Lagrangian descriptors were used to plan its mission in real time. All of this work was originally reported in RGGM$^{+}$18.

Silbo was deployed on the 13th of April 2016 near Falmouth, Massachusetts and was recovered on March 9th 2017 just south of Ireland, a 6506 km journey. An illustration of the path of Silbo is given in Fig. 13. Monitoring and navigation of the glider along a selected route is carried out at a collection of locations along the route that are selected ahead of time. These locations are referred to as waypoints. Lagrangian descriptors played the key role in selecting the waypoints that led to the remarkable success of the journey of Silbo in a way that we will now describe.

Silbo began its journey in the Gulf Stream so as to take advantage of this well-known current. However, gliders typically have a low horizontal speed (0.2–0.4 m/s) and this makes them extremely sensitive to the current fields that they experience. As a result, they are not capable of following an arbitrarily specified trajectory to reach a required location. Consequently, it would be a significant advantage to utilize the inherent geometrical patterns of the ocean current field in a way that would enable us to design trajectories that optimize the mission of the glider. This idea is similar in spirit to previous efforts that exploited the underlying dynamics generated by the gravitational fields of planetary systems for planning spacecraft missions with low thrust capabilities KLMR01MR06.

Copernicus Marine Environment Monitoring Service (CMEMS) provides daily velocity fields for the region of the North Atlantic where Silbo flew. Two dimensional velocity fields were provided in 50 vertical geopotential levels ranging from 0 to 5500 meters. In particular, to perform the Lagrangian path planning simulations, daily operational velocity fields have been derived from the data set by vertically averaging the velocity fields over the water column that extends from 0 to 902 meters depth, which was Silbo’s diving depth. The velocity fields had a 10 day forecast window. From these daily operational velocity fields Lagrangian descriptors were computed every 2 or 3 days using an 8 day forecast window which corresponded to setting $\tau = 8$, both forward and backward.

Using the Lagrangian descriptor fields, hyperbolic trajectories and their stable and unstable manifolds could be identified along the route of Silbo and these hyperbolic trajectories were chosen as waypoints. The stable and unstable manifolds of the hyperbolic trajectories provided natural navigation routes, as the stable manifolds provided optimal paths towards the waypoint and the unstable manifolds provided optimal paths away from the waypoints. In this way the Lagrangian descriptor maps were used to monitor the trajectory of Silbo and decide where the next waypoint should be placed.

In Fig. 14 we show Lagrangian descriptor maps in grayscale on the 20th and 23rd June 2016. We can see in the figure that the singular features of the Lagrangian descriptor field intersect in a hyperbolic trajectory, and correspond to the stable and unstable manifolds of the hyperbolic trajectory. The magenta arrows indicate the velocity field and confirm the direction of motion of fluid parcels and the saddle-like structure corresponding to the hyperbolic trajectory and its stable and unstable manifolds. Silbo’s path is denoted by the green curve and the speed of Silbo is indicated at the end of the curve.

Figure 13.

The path of Silbo as it crossed the northern Atlantic. This figure is from RGGM$^{+}$18 where details of the creation of this figure can be found.

Figure 14(a) shows that the speed of Silbo is reduced as it approaches the hyperbolic trajectory along its stable manifold on 20th June 2016. Figure 14(b) displays Silbo moving away from the hyperbolic point along its unstable manifold on the 23rd June 2016 and gaining speed. Remarkably, on that day Silbo achieved a record speed (1.04 m/s), which is unprecedented for this type of glider mission since typical speeds for this type of gliders are below 0.5 m/s.

Figure 14.

Glider path and Eulerian velocity fields in the neighborhood of a hyperbolic trajectory highlighted by the Lagrangian descriptor. (a) 20th June 2016. (b) 23rd June 2016. This figure is from RGGM$^{+}$18 where details about its creation can be found.

These results illustrate that stable manifolds (SM) are optimal paths towards the hyperbolic trajectory (HT), i.e., Silbo approaches the HT most closely along this direction. Similarly, unstable manifolds (UM) are the optimal paths for moving away from the HT. In the neighborhood of the HT, the glider slows down. Putting these observations together we conclude that the optimal path for navigating near a HT would be to remain close to the sequence of dynamical objects SM-HT-UM. These observations also indicate why HTs are good choices for waypoints. They are points where we can control the direction, and to some extent speed, of the glider, with one caveat.

In order to avoid a significant decrease in speed as the glider approaches the HT, a new waypoint is chosen at a different hyperbolic trajectory. This new HT is selected so that its SM is aligned with the UM of the previous HT. In this way, the new waypoint forces the glider to leave the neighborhood of the previous HT along the direction of the UM which, in general, would lead to an increase in speed of the glider. An appropriate navigation sequence thus would concatenate the dynamical objects SM-HT-UM/SM-HT-UM. This will result in a wave-like path with the glider moving alternatively from stable to unstable manifolds at the chosen waypoints.

Conclusions and Future Challenges

We have shown that Lagrangian descriptors can be applied in broad and diverse settings that are continually expanding. They can be applied to both deterministic and stochastic dynamical systems, low and high dimensional systems, and dynamical systems defined as data sets. For example, recently Lagrangian descriptors have been at the center of a new approach to understand the geometrical structure of uncertainty quantification GSMW22, in the area of machine learning they have been used as a feature for training support vector classifiers NKW21, and in chemical reaction dynamics they have been central to achieving a long sought after theory of reactive islands in Hamiltonian systems with three degrees-of-freedom KGGW21. These latter two examples highlight how Lagrangian descriptors can be applied to provide insight into high dimensional dynamics. The flexible, easily implementable and interpretable nature of this computational tool results in a very low learning curve for applied scientists outside of mathematics, and this has undoubtedly inspired the growth and breadth of applications that use this tool. We hope to encourage this further by providing an open source software package ASGGK$^{+}$21, complete with tutorial examples, that show how to analyze specific problems in the spirit of all of the topics described above.

https://github.com/champsproject/ldds