${x}_{r}=\frac{i-1}{628}2\text{\pi}$ (3)

${y}_{r}=\frac{j-1}{628}2\text{\pi}$ (4)

To identify the hyperbolic trajectories, we use a filter to apply a threshold to our phase space and highlight regions that undergo large dispersion. For a given grid point, $\left(i,j\right)$ , and threshold value, $\tau $ , the filtered matrix of dispersion values is given by,

$\varphi \left({R}_{ij}\right)=(\begin{array}{cc}1& {R}_{ij}\ge \tau \\ 0& {R}_{ij}<\tau \end{array}$ (5)

Figure 1. Phase portrait of the Standard map after 12 iterations for $K=0.971635$ .

The resulting binary, filtered matrix consists of ones in the $\left(i,j\right)$ locations that correspond to relative dispersion values above a given threshold and zeroes elsewhere. The threshold is determined relative to the maximum possible dispersion value. Converting the resulting grid points from $\varphi \left({R}_{ij}\right)$ intensity space to $\left({x}_{r}\mathrm{,}{y}_{r}\right)$ space yields the path of maximum transport relative to neighboring points for the map. Based on distance arguments, relative dispersion on the $\left[\mathrm{0,2}\right)$ torus is bound. For the point located at $\left(x\mathrm{,}y\right)$ , there are four neighboring points that are followed in iteration and used to calculate relative dispersion. The locations of the neighboring points are provided in Figure 2 and Figure 3. On the torus, separation distance in the horizontal direction is bound by the modulo value. Given initial particle locations $\left({x}_{i-\mathrm{1,}j}^{n}\mathrm{,}{y}_{i-\mathrm{1,}j}^{n}\right)\mathrm{,}\left({x}_{i+\mathrm{1,}j}^{n}\mathrm{,}{y}_{i+\mathrm{1,}j}^{n}\right)$ (Figure 2), each particle is advected with iteration to step, N. The resulting dispersion of the particles that were initialized in the horizontal direction relative to the center point is bound as follows:

${R}_{h}=\sqrt{{\left({x}_{i+1,j}^{N}-{x}_{i-1,j}^{N}\right)}^{2}+{\left({y}_{i+1,j}^{N}-{y}_{i-1,j}^{N}\right)}^{2}}$ (6)

Similarly, for points $\left({x}_{i\mathrm{,}j-1}^{N}\mathrm{,}{y}_{i\mathrm{,}j-1}^{N}\right)$ and $\left({x}_{i\mathrm{,}j+1}^{N}\mathrm{,}{y}_{i\mathrm{,}j+1}^{N}\right)$ , the separation distance after advection is bound in the vertical direction by:

${R}_{v}=\sqrt{{\left({x}_{i,j+1}^{N}-{x}_{i,j-1}^{N}\right)}^{2}+{\left({y}_{i,j+1}^{N}-{y}_{i,j-1}^{N}\right)}^{2}}$ (7)

For the Standard map:

${R}_{h}<2\text{\pi}\sqrt{2},\text{\hspace{1em}}{R}_{v}<2\text{\pi}\sqrt{2}$ (8)

$R={R}_{h}+{R}_{v}$ (9)

In other versions of this method, the dispersion value is defined relative to the initial dispersion of the particles in the x and y directions. To apply the method to the Standard map, we use a uniformly spaced, square grid to generate the initial particle locations. Setting the initial dispersion in the horizontal and vertical directions to ${10}^{-2}$ , we can consider only the dispersion after iteration. Using Equations 8 and 9 the relative dispersion response can be bound by

$R<\sqrt{32{\text{\pi}}^{2}}\approx 17.771$ (10)

We compare the R value at each location, $\left(i\mathrm{,}j\right)$ , with the upper bound of R values to generate our threshold as a percentage:

$Thresh=\frac{R}{\mathrm{max}\left(R\right)}$ (11)

In applying the method of relative dispersion to the Standard map, we see that points that are initialized within one of the chaotic islands of the Standard map (below the line $y=Kx$ in the interval $x\in \left(\mathrm{0,3}\right)$ ) remain entrained within the circulation that travels between the upper and lower halves of the phase space. The particles that begin above the line $y=Kx$ are trapped outside of the chaotic islands and travel more as a packet, undergoing less dispersion relative to each other.

For initial positions near the hyperbolic trajectory as in Figure 3, the particles may follow different travel paths. Consider the particles on the horizontal axis, the particle on the left starts inside of the chaotic region, while the particle on the right is largely restricted to travel within the non-chaotic island. Similarly, in the vertical direction, the particle on the bottom will travel within the non- chaotic region while the corresponding particle at the top follows a path within the chaotic region. The result is a larger relative dispersion between points that begin in the vicinity of the central grid point (the lighter blue regions of Figure 4). In instances where the quartet of points is initialized within the chaotic regions (or between the chaotic regions) dispersion is minimized, resulting in lower relative dispersion values (the darker blue regions of Figure 4).

3. Relative Dispersion in a Rossby Wave

A system representing an autonomous Rossby wave is used to further assess how well relative dispersion identifies hyperbolic trajectories. The use of the Rossby wave helps provide insight into the method’s impact and validity relative to dynamical systems. The general formulation of the Rossby wave under consideration:

$\frac{\text{d}y}{\text{d}t}={a}_{0}\mathrm{sin}\left(kx\right)+{b}_{0}\mathrm{sin}\left(\omega t\right)$ (12)

$\frac{\text{d}x}{\text{d}t}=y$ (13)

The equation is often used to simulate gyre motion. The system as used can be considered a simple model for a double gyre flow (Figure 5). The Rossby wave is selected as an example of a physical, Hamiltonian system, providing some additional tools to aid in evaluation. Like the Standard map, the differential

Figure 2. Relative locations of grid points used for calculating Relative Dispersion at $\left(x\mathrm{,}y\right)$ .

Figure 3. Example of initial (squares) and final (triangles) positions of the four points used to find ${R}_{ij}$ . In the vicinity of the hyperbolic trajectory, some points may become trapped in a circulation while others are transported longer distances.

equation uses a parameter, ${a}_{0}$ , to control the non-linear “kick” for the system. The parameter, ${b}_{0}$ controls the time dependence [3] . The time independent case can be investigated by setting ${b}_{0}=0$ . The flow field and relative dispersion derived response for the time independent case are provided (Figure 6). We use the time independent case to assess the underlying mechanics and the ability of relative dispersion to identify the hyperbolic trajectories, as well as some elliptic orbits. A benefit of the selection of the representation is the existence of a closed form solution for the time-independent case:

${\psi}_{0}=\frac{1}{2}{y}^{2}-\frac{{a}_{0}}{k}\mathrm{cos}\left(kx\right)+c$ (14)

4. Comparisons

We investigate the distribution of the relative dispersion values to determine optimal threshold values for isolating specific trajectories. The histograms in Figure 7, and Figure 8 depict the frequency distributions of dispersion values for the 395,641 grid points used. The majority of the points that are used generate a dispersion value, $R<2$ . Values of R in this range denote grid points whose neighbors do not travel large distances apart, and whose travel is more in the form of a packet. In focusing only on particles with larger dispersion values individual trajectories can be identified. An additional perk of the process allows for additional features to be extracted. Isolating mid-ranges of the R-values makes it possible to identify other orbits within the map. For the Standard map, the range of possible dispersion values is bound based on the periodicity of the

Figure 4. Relative dispersion derived portrait of the Standard map for $K=0.7$ (sub- chaotic), 0.971 (chaos threshold), and 1.3 (chaotic).

Figure 5. The intensity phase portrait and isoclines of the example Rossby wave for ${a}_{0}=1$ , ${b}_{0}=0$ , and $k=1$ generated by the analytic solution to the system.

map. There are also a large number of grid points whose relative dispersion value falls in the interval $R\in \left(\mathrm{2,7}\right)$ . While the number of grid points that result in an $R>7$ decreases, the cumulative effect of relative dispersion for all values greater than $R=7$ results in poor isolation of the trajectory. Other orbits are well defined at that value while values of $R>12$ result in visually and

Figure 6. The flow field for a simple, time-independent Rossby wave and the relative dispersion response of the system.

quantitatively closer approximations to the stable and unstable manifolds of the Standard map. The distribution of the $628\times 628$ grid of relative dispersion values is binned in Figure 7. Each spectrum of blue to dark red represents the 628 columns of the relative dispersion intensity matrix. The horizontal axis uses the minimum and maximum dispersion values, $\left(0,\sqrt{32{\text{\pi}}^{2}}\right)$ (10) to partition the range of dispersion values into ten intervals. The vertical axis shows how many grid points in each column returned a dispersion value in the given range. Each column is represented by the same color in each interval. We find that the majority of the particles that are initialized at a given grid point undergo dispersion of less than 2 units. We also find that a small subset of particles undergoes large scale dispersion, returning a dispersion assessment of $R\ge 14$ .

To help evaluate how well the relative dispersion derived trajectories approximate the traditionally derived trajectories, Jacobian analyses are per- formed on each system. Using the results classifies and provides the orientations for the trajectories corresponding to the fixed point.

In the Standard map, there are fixed points at $\left(2\text{\pi}n,0\right)$ and at $\left(\left(2n-1\right)\text{\pi}\mathrm{,0}\right)$ . A Jacobian analysis classifies the points as hyperbolic fixed points and centers, respectively.

$J=\left(\begin{array}{cc}0& 1\\ k\mathrm{cos}x& 0\end{array}\right)$ (15)

$J\left(0,0\right):\to {\lambda}_{1}=\sqrt{k},\text{\hspace{0.17em}}{\lambda}_{2}=-\sqrt{k}$ (16)

$J\left(\left(2n-1\right)\text{\pi},0\right):\to {\lambda}_{1}=i\sqrt{k},\text{\hspace{0.17em}}{\lambda}_{2}=-i\sqrt{k}$ (17)

In the Rossby wave, a Jacobian analysis of the fixed points for $k>0$ on the interval $x\in \left(-2\text{\pi}\mathrm{,2}\text{\pi}\right)$ returns 2 real eigenvalues at $\left(\mathrm{0,0}\right)$ and 2 purely imaginary eigenvalues at $\pm \left(\text{\pi}\mathrm{,0}\right)$ . The fixed points at $\left(2\text{\pi}n,0\right)$ are hyperbolic, while the fixed points corresponding to $\left(\left(2n-1\right)\text{\pi}\mathrm{,0}\right)$ result in centers.

$J=\left(\begin{array}{cc}0& 1\\ {a}_{0}k\mathrm{cos}x& 0\end{array}\right)$ (18)

$J\left(0,0\right):\to {\lambda}_{1}=\sqrt{{a}_{0}k},\text{\hspace{0.17em}}{\lambda}_{2}=-\sqrt{{a}_{0}k}$ (19)

Figure 7. The top image displays the values returned by relative dispersion applied to the Standard map, uniformly parsed into 10 bins. The majority of grid points return a value less than 7 (40% of $\mathrm{max}\left(R\right)$ ). The bottom image shows the 394 ,384 points grouped based on the dispersion value that is generated for their initial locations.

$J\left(\left(2n-1\right)\text{\pi},0\right):\to {\lambda}_{1}=i\sqrt{{a}_{0}k},\text{\hspace{0.17em}}{\lambda}_{2}=-i\sqrt{{a}_{0}k}$ (20)

With a cursory glance, one can anticipate lower dispersion values to be generated when the four points used are contained inside of either of the two constrained regions (Figure 5). Similarly, if all four points are located outside of the regions, the expectation is that they travel more as packet than not. When some of the points are entrained within one of the regions centered at $\left(-\text{\pi}\mathrm{,0}\right)$ or $\left(\text{\pi}\mathrm{,0}\right)$ , and the others are not, larger dispersion values are generated.

Though a closed form solution exists for the system, the solution for all systems may not be easily found, if it is possible at all. In the case of observed data, the equations generating the motion are often unknown. The accuracy of the method must be investigated relative to data that is typically available. As in the case of the Standard map, the assessment is conducted by identifying a comparison hyperbolic trajectory, densely discretizing along it, and iterating in time. To generate the comparison trajectory, points are initialized along the hyperbolic trajectory for $x\in \left(\mathrm{0,0.1}\right)$ , requiring the sampling rate $\delta x\ll 0.1$ and iterated forward in time. In the trial, $\delta x=0.00001$ . After iteration, the distance between neighboring points along the trajectory is calculated, $\delta {x}^{\prime}$ . The distance between any two successive points after iteration is required to be less than the grid discretization that is used in calculating the relative dispersion values $\left(0.01\right)$ .

To compare the data, a rigid transformation from $\left(i\mathrm{,}j\right)$ space to $\left(x\mathrm{,}y\right)$ space must again be applied to the data. Due to the difference in the domain size (as compared to the Standard map), a different scaling is required. The system is analyzed on the interval $x\in \left(-2\text{\pi}\mathrm{,2}\text{\pi}\right)$ and $y\in \left(-\text{\pi}\mathrm{,}\text{\pi}\right)$ . The conversion from $\left(i\mathrm{,}j\right)$ to $\left(x\mathrm{,}y\right)$ space is

$x=4\text{\pi}\frac{i-1}{{N}_{x}}-2\text{\pi}y=2\text{\pi}\frac{j-1}{{N}_{y}}-\text{\pi}$ (21)

where i and j are both integers and $i\in \left[\mathrm{1,}{N}_{x}\right]$ and $j\in \left[\mathrm{1,}{N}_{y}\right]$ .

5. Results

We performed a relative dispersion assessment for several values of K for the Standard map. Recall the parameter K controls the chaos in the map. Chaotic behavior in the map is traditionally associated with values of K above 0.971. We present the relative dispersion assessments for $K=0.971635$ . Each map was iterated for $N=16$ steps to normalize the comparison. The trajectory of focus can be seen starting near $\left(\mathrm{0,0}\right)$ , peaking near $\left(4.4,2\right)$ and oscillating as it approaches $\left(\mathrm{2,0}\right)$ . A full discussion of the features and dynamics of the Standard map can be found in numerous articles and books, and is omitted here. Examples of the intensity images returned from relative dispersion are pre- sented in Figure 9. For comparison, the same images are overlaid with the corresponding hyperbolic trajectories that are used for validation (in red). Per initial observations, relative dispersion seems to provide very good alignment with the manifolds.

From our assessment of R distribution (Figure 7), we found that values of $R>12$ , are the optimal focal points for efforts to delineate the key trajectories. We also found that values of $R>16$ , provided accurate but very sparse data for approximating the desired orbits. While values closer to $\mathrm{max}\left(R\right)$ would seem ideal, the accuracy relative to location that is gained from a higher threshold

Figure 8. The top image displays values returned by relative dispersion for a Rossby wave, uniformly parsed into 10 bins. The majority (95%) of grid points return a value less than 5 ( $\stackrel{\u02dc}{1}\mathrm{8\%}$ of $\mathrm{max}\left(R\right)$ ). The bottom image shows the total number of grid points returning a given dispersion value. The 31,752 points are grouped based on the dispersion value that is generated for their initial locations.

value sacrifices density of data points used to delineate the trajectory. The loss of additional data points introduces room for interpolation errors. Therefore we have found that there is an inherent balancing act that must be performed by the investigator based on the desired outcome and intended use of the data. Values of $R\in \left(\mathrm{14,16}\right)$ were found to provide the best balance between accuracy of approximation and density of data.

In an effort to generalize R thresholds for other maps and flows, the threshold has also been investigated as a percentage relative to the maximum dispersion value. We found values of $R\in \left(75\%\mathrm{max}\left(R\right),90\%\mathrm{max}\left(R\right)\right)$ provide a general form for the optimal R threshold values for delineating both the trajectories of the Standard map. These values are equivalent to an interval of $R\in \left(\mathrm{12.8147,15.377}\right)$ .

Figure 9. The Relative Dispersion generated intensity field for the Standard map.

Figure 10. The Relative Dispersion generated intensity field, the locations that return a value above 80% of the maximum value, and the locations returning 80% of the maximum value superimposed on the orbits generated by a 4th order Runge- Kutta method.

In Figure 9, the dark blue regions correspond to small dispersion values. Lighter blue and white correspond to larger dispersion values. The grid points that satisfy a given threshold value are highlighted in green. In each image, the corresponding trajectory (stable or unstable) is superimposed in red. In many of the images, the green is difficult to observe due to its (desired) coincidence with the superimposed, comparison trajectory, in red.

A similar assessment is applied to the Rossby wave example. The Rossby example did not require as high a threshold as the Standard map (Figure 8). A threshold that isolates particles that have traveled in the top 50% of total

Figure 11. The Relative Dispersion generated intensity field, the locations that return a value above $R=15$ , and the locations returning a value above $R=15$ superimposed on the hyperbolic trajectories generated by the map.

Figure 12. The sum of registration error values from the approximation to the hyperbolic trajectories of a Standard map as a function of iteration. The locations of highest dispersion value are compared with the values derived from the map definition.

distances eliminates many of the excess paths, while increasing to 80% refines the results. Increasing to 90% does not improve the results in the unstable directions, while losing some information in the stable directions (Figure 10).

6. Discussion

The relative dispersion response for the Standard map is presented for dispersion

Figure 13. The Relative Dispersion generated intensity field for the Rossby wave.

Figure 14. The sum of registration errors for all points generated by the relative dispersion approximation to the hyperbolic trajectories of a Rossby wave, as a function of time. The locations of highest dispersion value are compared with the locations derived from a fourth order Runge-Kutta implementation. As the system evolves, the relative dispersion response converges to expected trajectories.

values greater than 15 (the maximum possible is approximately 17.77) (Figure 11). The relative dispersion approximation for the Standard map example is approximately second order accurate, after about 30 iterations of the map. In early steps of the iteration, the evolution of the system is in its infancy resulting in less dispersion of neighboring points, whether near the hyperbolic trajectory or not. The result is higher error in locating the desired orbit. As time progresses, there is an averaging out effect, as points near the trajectory undergo larger relative dispersion than point quartets that are wholly located either inside or outside of the chaotic islands of the Standard map. The mean error is second order accurate. The dispersion error is calculated as the system evolves and the relative dispersion response is calculated, which contributes to the error convergence (Figure 12). The results may allow one to determine how long in time or iteration is necessary to allow a system to run to return sufficient representation of its behavior.

When only the final time step is considered, and the per point error is plotted, additional information can be isolated. Initial locations near fixed points in the stable direction immediately undergo dispersion, and it grows with each time step. Values that begin near the unstable trajectory travel as a packet in the short term and undergo increasing amounts of dispersion as they near the next fixed point.

In the Rossby wave example, the oscillations in the stable direction that occur in the Standard map are not present (Figure 13). Upon inspection, the relative dispersion generated approximation appears to be better. The mean error is improved, but more importantly is also on the scale of the discretization. The method for choosing the threshold used for delineation of the desired trajectory is updated.

Instead of using a threshold involving an arbitrarily chosen dispersion value, the threshold is defined in terms of a percentage of the maximum dispersion value. All of the dispersion values are considered and only central locations that return 80% of the maximum dispersion value over the time period are used. The response from relative dispersion is similar to that of the Standard map, though with better delineation in the unstable directions. Using a filter of 80%, there is good consistency in the localization of the points. Superimposing the responses from the relative dispersion implementation, the Runge-Kutta derived trajectory, and the flow field generated by the equations results in good alignment upon inspection (Figure 10). The relative dispersion response shows better delinea- tion of the desired trajectories than in the Standard map example.

To gauge the degree of accuracy of the approximation, an error calculation is performed. For each time step, the total registration error between the points generated by relative dispersion and Runge-Kutta is computed. The error decreases with time. There is very quick convergence of the relative dispersion response to the Runge-Kutta response, which can aid investigators in deter- mining the appropriate length of time to run the model (Figure 14). Depending on the needs of the project, the model may be able to be run for as few as 10 time steps rather than 45 depending on desired accuracy.

In the future, an application to other time independent systems would provide further insight into the viability of the method for identifying stable and unstable trajectories and other orbits. While relative dispersion can and has been applied to time dependent systems, additional methods are needed for evaluating the results and may be informative. Modifying the parameters of the systems to induce additional dynamics will give insight into the robustness of the method. In the larger scheme, expansion to three dimensional systems and an assessment of the method’s ability to accurately identify surfaces is a major goal.

Cite this paper

Redd, T. , Moghraby, A. and Tillman, T. (2017) Assessing Relative Dispersion.*Applied Mathematics*, **8**, 1572-1589. doi: 10.4236/am.2017.811115.

Redd, T. , Moghraby, A. and Tillman, T. (2017) Assessing Relative Dispersion.

References

[1] Haller, G. and Poje, A. (1999) Geometry of Cross-Stream Mixing in a Double-Gyre Model. Journal of Physical Oceanography, 29, 1649-1665.

https://doi.org/10.1175/1520-0485(1999)029<1649:GOCSMI>2.0.CO;2

[2] Branicki, M. and Wiggins, S. (2010) Finite-Time Lagrangian Transport Analysis: Stable and Unstable Manifolds of Hyperbolic Trajectories and Finite-Time Lyapunov Exponents. Nonlinear Processes in Geophysics, 17, 1-36.
https://doi.org/10.5194/npg-17-1-2010

[3] Malhotra, N. and Wiggins, S. (1998) Geometric Structures, Lobe Dynamics, and Lagrangian Transport in Flows with Aperiodic Time-Dependence, with Applications to Rossby Wave Flow. Journal of Nonlinear Science, 8, 401-456.
https://doi.org/10.1007/s003329900057

[4] Ottino, J. (1989) The Kinematics of Mixing: Stretching, Chaos, and Transport. Cambridge UP, Cambridge.

[5] Rom-Kedar, V., Leonard, A. and Wiggins, S. (1999) An Analytical Study of Transport, Mixing, and Chaos in Unsteady Vortical Flow. Journal of Fluid Mechanics, 214, 347-394.

https://doi.org/10.1017/S0022112090000167

[6] Miller, P., Jones, C., Rogerson, A. and Pratt, L. (1997) Quantifying Transport in Numerically-Generated Velocity Fields. Physica D, 110, 105-122.

[7] Madrid, J. and Mancho, A. (2009) Distinguished Trajectories in Time Dependent Vector Fields. Chaos, 19, Article ID: 013111. https://doi.org/10.1063/1.3056050

[8] Hamilton, P., Larsen, J., Leaman, K., Lee, T. and Waddell, E. (2004) Transports through the Straits of Florida. Journal of Physical Oceanography, 35, 308-322.
https://doi.org/10.1175/JPO-2688.1

[9] Abraham, E. and Bowen, M. (2002) Chaotic Stirring by a Mesoscale Surface-Ocean Flow. Chaos, 12, 373-381. https://doi.org/10.1063/1.1481615

[10] Poje, A., Toner, M., Kirwan, A. and Jones, C. (2002) Drifter Launch Strategies Based on Lagrangian Templates. Journal of Physical Oceanography, 32, 1855-1869.

https://doi.org/10.1175/1520-0485(2002)032<1855:DLSBOL>2.0.CO;2

[11] Wiggins, S. (1992) Chaotic Transport in Dynamical Systems. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4757-3896-4

[12] Kuznetsov, L., Jones, C., Toner, M. and Kirwan, A. (2004) Assessing Coherent Feature Kinematics in Ocean Models. Physica D: Nonlinear Phenomena, 191, 81-105.

[13] Bowman, K. (1999) Manifold Geometry and Mixing in Observed Atmospheric Flows.

[14] Haller, G. (2000) Finding Finite-Time Invariant Manifolds in Two Dimensional Velocity Fields. Chaos, 10, 99-108. https://doi.org/10.1063/1.166479

[15] Von Hardenberg, J., Fraedrich, K., Lunkeit, F. and Provenzale, A. (2000) Transient Chaotic Mixing during a Baroclinic Life Cycle. Chaos, 10, 122-134.
https://doi.org/10.1063/1.166491

[16] Jones, C. and Winkler, S. (2002) Invariant Manifolds and Lagrangian Dynamics in the Ocean and Atmosphere. Handbook of Dynamical Systems, 2, 55-92.

[17] Froyland, G. and Padberg, K. (2009) Almost-Invariant Sets and Invariant Manifolds-Connecting Probabilistic and Geometric Descriptions of Coherent Structures in Flows. Physica D, 238, 1507-1523.

[18] Klocek, D. (2005) Weather and Chaos Theory: The Standard Map.

http://docweather.com/4/show/124/

[19] Gelfreich, V. (1999) A Proof of the Exponentially Small Transversality of the Separatrices for the Standard Map. Communications in Mathematical Physics, 201, 155-216.

https://doi.org/10.1007/s002200050553

[20] MacKay, R. and Meiss, J. (1987) Hamiltonian Dynamical Systems: A Reprint Selection. Taylor & Francis, New York.