# Local Trajectory Variation Exponent (LTVE) for
Visualizing Dynamical Systems

Yunchen Tsai Department of Mathematics, the Hong Kong University of Science
and Technology, Clear Water Bay, Hong Kong. Email<EMAIL_ADDRESS>Shingyu Leung Department of Mathematics, the Hong Kong University of Science
and Technology, Clear Water Bay, Hong Kong. Email<EMAIL_ADDRESS>
###### Abstract

The identification and visualization of Lagrangian structures in flows plays a
crucial role in the study of dynamic systems and fluid dynamics. The Finite
Time Lyapunov Exponent (FTLE) has been widely used for this purpose; however,
it only approximates the flow by considering the positions of particles at the
initial and final times, ignoring the actual trajectory of the particle. To
overcome this limitation, we propose a novel quantity that extends and
generalizes the FTLE by incorporating trajectory metrics as a measure of
similarity between trajectories. Our proposed method utilizes trajectory
metrics to quantify the distance between trajectories, providing a more robust
and accurate measure of the LCS. By incorporating trajectory metrics, we can
capture the actual path of the particle and account for its behavior over
time, resulting in a more comprehensive analysis of the flow. Our approach
extends the traditional FTLE approach to include trajectory metrics as a means
of capturing the complexity of the flow.

##  1 Introduction

The study of coherent structures has been a fundamental area of research in
dynamic systems and computational fluid dynamics. Among the various types of
coherent structures, Lagrangian coherent structures (LCS) are particularly
noteworthy. LCS focuses on identifying regions within a flow that display
strong attraction and repulsion of particles over a specific time span. It
serves as a valuable tool for analyzing, visualizing, and extracting
information from complex dynamical systems. The primary objective is to
delineate surfaces of trajectories and separate regions in the phase space
that exert significant influence on neighboring trajectories during a defined
time interval [15]. LCS has found applications in diverse fields, including
oceanography [18, 31], meteorology [29], flight mechanics [4, 33], gravity
waves [35], and bio-inspired fluid flows [23, 11, 24]. Although the underlying
physical intuition behind LCS is straightforward, an accurate and formal
mathematical definition capturing this behavior is still a subject of ongoing
development. Numerous studies have attempted to establish robust frameworks
for identifying LCS structures [14, 31], but these frameworks often face
computational challenges or suffer from certain limitations that affect their
precision.

The finite-time Lyapunov exponent (FTLE) is a widely used measure for locating
Lagrangian coherent structures (LCS) [16, 12, 13, 31, 19]. It quantifies the
rate of change in distance between neighboring particles over a finite time
interval, considering an infinitesimal perturbation in the initial position.
Computing the FTLE field involves several steps. Firstly, the flow map is
computed, which establishes the connection between the initial and arrival
positions of particles along the characteristic line. The FTLE is then defined
based on the largest eigenvalue of the deformation matrix derived from the
Jacobian of the flow map. Various Eulerian approaches have been developed to
numerically compute the FTLE on a fixed Cartesian mesh in recent studies . For
more in-depth discussions on this topic, interested readers can refer to [20,
21, 47, 42, 43, 22, 27, 44, 45, 46] and related references.

Another approach to analyzing trajectories involves trajectory analysis, which
aims to develop techniques for understanding and classifying trajectory
characteristics. While trajectory clustering is not a new topic, previous
methods primarily focused on visualizing or constructing dictionaries and
summaries to reveal hidden patterns or predict future routes [17, 38, 10, 9,
39, 8, 40]. These approaches were not specifically designed to discover
coherent structures within underlying flows and dynamical systems. However,
recent works have explored clustering-based approaches for identifying LCS.
Some of these methods involve constructing networks of trajectories and
determining edge weights based on trajectory similarity [30, 25, 28]. In [41],
a coherent ergodic partition method was developed to separate trajectories
into clusters by integrating and computing long-time averages of functions
along the trajectories. This effectively projects high-dimensional data onto a
low-dimensional manifold, allowing for partitioning based on these function
averages. In [6], a clustering-based approach was proposed to identify
coherent flow structures in continuous dynamical systems. This method treats
particle trajectories over a finite time interval as high-dimensional data
points and clusters them from different initial locations into groups using
normalized standard deviation or mean absolute deviation to quantify
deformation.

The aim of this paper is to propose a general framework based on the
trajectory analysis technique to the study of LCS. Specifically, we seek to
develop a new metric for efficiently extracting and visualizing LCS using
trajectory analysis. Instead of clustering all provided trajectories into
groups, we investigate the variation of the trajectories in a small
neighborhood of each trajectory. We define a local quantity that measures the
similarity of these trajectories from a small local neighborhood based on
various possible trajectory metrics. This metric is a commonly used tool for
quantifying the similarity between trajectories, and its metric properties may
prove useful in developing the theoretical framework [32, 36, 26]. Our
proposed metric considers the effect of a small perturbation on the entire
trajectory over a fixed time period, rather than just the final position,
which we believe provides a more robust method for extracting LCS over an
extended period.

This paper is organized as follows: Section 2 presents the related works and
background, while Section 3 details the model and algorithm for our proposed
quantity. Section 3.6 discusses several possible choices of the trajectory
metric. Finally, Section 4 provides several numerical results.

##  2 Background

In this section, we provide several definitions and assumptions on our problem
that is crucial for the discussion in the later section.

###  2.1 Trajectory Metric

Consider a closed and bounded domain $\Omega\subset\mathbb{R}^{n}$, a finite
time domain $[t_{0},t_{0}+T]$ and a Lipchitz velocity field
$\mathbf{v}:\Omega\times\mathbb{R}\mapsto\mathbb{R}^{n}$. The movement of any
particle in $\Omega$ can be described by the following ordinary differential
equation (ODE),

$\displaystyle\dot{\mathbf{x}}(t)$ $\displaystyle=$
$\displaystyle\mathbf{v}(\mathbf{x}(t),t)\quad\forall\mathbf{x}\in\Omega,t\in[t_{0},t_{0}+T]$
(1)

for any initial condition $\mathbf{x}(t_{0})=\mathbf{x}_{0}\in\Omega$.

######  Definition 2.1 (Trajectory of Particle).

def_traj Given $\mathbf{x}_{0}\in\mathbb{R}^{n}$, a trajectory
$\mathbf{x}:\mathbb{R}\mapsto\Omega$ is defined as the solution of the ODE in
(1) with initial condition of $\mathbf{x}(t_{0})=\mathbf{x}_{0}$.

Instead of simply measuring the distance between two trajectories, we are more
interested in studying the intrinsic difference in the motion characteristics.
One idea is to consider only the so-called displacement trajectory, i.e. the
relative displacement from the initial location at the initial time $t=t_{0}$.

######  Definition 2.2 (Displacement Trajectory of Particle).

Given $\mathbf{x}_{0}\in\Omega$ with the trajectory
$\mathbf{x}:\mathbb{R}\mapsto\Omega$ as defined in def_traj, a displacement
trajectory $\Psi_{\mathbf{x}_{0}}:\mathbb{R}\mapsto\mathbb{R}^{n}$ is defined
as $\Psi_{\mathbf{x}_{0}}(t):=\mathbf{x}(t)-\mathbf{x}(t_{0})$.

Since we work with a discretized time domain, we consider the following
discrete displacement trajectory.

######  Definition 2.3 (Discrete Trajectory).

def_discrete Given a trajectory $\Psi$, then a $k$-discrete trajectory is
defined as $\Psi_{\Delta}:[0,k-1]\cap\mathbb{N}\mapsto\Omega$ such that
$\Psi_{\Delta}(j)=\Psi\left(t_{0}+{jT}/{(k-1)}\right)$ for $j=0,1,\cdots,k-1$.

We can map any $k$-discrete trajectory into a unique vector in the space of
$\mathbb{R}^{kn}$ by

$L(\Psi_{\Delta})=(\pi_{e_{1}}(\Psi_{\Delta}(0)),\pi_{e_{2}}(\Psi_{\Delta}(0)),...,\pi_{e_{d}}(\Psi_{\Delta}(k-1)))$

where $\pi_{e_{j}}(\mathbf{v})$ is the projection map that projects
$\mathbf{v}$ to $e_{j}$ with $e_{j}$ denotes the $j$-th standard basis in
$\mathbb{R}^{n}$. One can quickly verify that the map $L$ is linear and
bijective. Thus the space of all $k$-discrete trajectory is isomorphic to the
space $\mathbb{R}^{kn}$. For the simplicity of presentation, we will abuse the
notation of $\Psi_{\Delta}$ for the rest of the article to represent both the
mapping $[0,k-1]\cap\mathbb{N}\mapsto\Omega$ and also the vector
$L(\Psi_{\Delta})$ whenever it is clear from the context.

To study the similarity between two displacement trajectories, we follow a
common approach to define a trajectory metric. Such a method is standard in
trajectory analysis, especially in classifying and clustering a set of
trajectories [32]. The definition is straightforward. Denoting $\Phi$ as the
space of all trajectories, we obtain a trajectory metric
$d:\Phi\times\Phi\mapsto\mathbb{R}$ such that $(\Phi,d)$ forms a metric space.

######  Definition 2.4 (Trajectory Metric).

$\Phi$ is the space of $\mathbf{C}^{1}$ function $\mathbb{R}\mapsto\Omega$,
then $d:\Phi\times\Phi\mapsto\mathbb{R}$ is a trajectory metric $\iff$ it
satisfies

  * •

$d(\Psi_{1},\Psi_{2})=0\iff\Psi_{1}=\Psi_{2}\quad\forall\Psi_{1},\Psi_{2}\in\Phi$

  * •

$d(\Psi_{1},\Psi_{2})=d(\Psi_{2},\Psi_{1})\quad\forall\Psi_{1},\Psi_{2}\in\Phi$

  * •

$d(\Psi_{1},\Psi_{3})\leq
d(\Psi_{1},\Psi_{2})+d(\Psi_{2},\Psi_{3})\quad\forall\Psi_{1},\Psi_{2},\Psi_{3}\in\Phi$

###  2.2 Finite-time Lyapunov Exponent (FTLE) and Lagrangian Coherent
Structure (LCS)

We define the flow map $f:\Omega\rightarrow\Omega$, which collects the
solutions to ODEs (1) for all initial conditions
$\mathbf{x}(t_{0})=\mathbf{x}_{0}\in\Omega$ at all times $t\in\mathbb{R}$. The
flow map represents the arrival location $\mathbf{x}(t_{0}+T)$ at $t=t_{0}+T$
of the particle trajectory with the initial condition
$\mathbf{x}(a)=\mathbf{x}_{0}$ at the initial time $t=t_{0}$. In other words,
this mapping takes a point from $\mathbf{x}(t_{0})$ at $t=t_{0}$ to another
point $\mathbf{x}(t_{0}+T)$ at $t=t_{0}+T$.

An important quantity in studying the separation rate between adjacent
particles with infinitesimal perturbations in their initial locations over a
finite-time interval is the finite-time Lyapunov exponent (FTLE). Let us
consider the initial time $t=t_{0}=0$ and the final time $t=T$. The change in
the initial infinitesimal perturbation is given by

$\delta\mathbf{x}(T)=f(\mathbf{x}+\delta\mathbf{x}(0))-f(\mathbf{x})=\nabla
f(\mathbf{x})\delta\mathbf{x}(0)+O(\|\delta\mathbf{x}(0)\|^{2})\,.$

The magnitude of the leading-order term of this perturbation is given by

$\|\delta\mathbf{x}(T)\|=\sqrt{\langle\delta\mathbf{x}(0),[\nabla
f(\mathbf{x})]^{\prime}\nabla f(\mathbf{x})\delta\mathbf{x}(0)\rangle}$

where $(\cdot)^{\prime}$ denotes the transpose of a matrix. We denote
$\varDelta^{T}_{0}(\mathbf{x})=[\nabla f(\mathbf{x})]^{\prime}\nabla
f(\mathbf{x})$ as the Cauchy-Green deformation tensor.

To find the maximum value of $|\delta\mathbf{x}(T)|$, we align
$\delta\mathbf{x}(0)$ with the eigenvector $\textbf{e}(0)$ associated with the
maximum eigenvalue of $\varDelta^{T}_{0}(\mathbf{x})$. Thus,
$\max_{\delta\mathbf{x}(0)}|\delta\mathbf{x}(T)|=\sqrt{\lambda_{\max}(\varDelta^{T}_{0}(\mathbf{x}))}|\textbf{e}(0)|=e^{\sigma^{T}_{0}(\mathbf{x})|T|}|\textbf{e}(0)|$.
Therefore, the FTLE $\sigma^{T}_{0}(\mathbf{x})$ is defined as

$\sigma^{T}_{0}(\mathbf{x})=\frac{1}{|T|}\ln\sqrt{\lambda_{\max}(\varDelta^{T}_{0}(\mathbf{x}))}\,.$
(2)

Numerical methods for computing the flow map $f(\mathbf{x})$ typically involve
ray tracing, where we solve the ODE (1) or the corresponding level set
equation [20] on a Cartesian mesh at the initial time $t=t_{0}=0$. Then, for
each grid point in the computational domain, we construct the Cauchy-Green
deformation tensor by finite differencing the forward flow map
$f(\mathbf{x}^{i}_{g})$ on the Cartesian mesh at $t=0$. Using equation (2), we
can compute the forward FTLE. Whether we use the Lagrangian (ODE) or the
Eulerian (PDE) approach, the FTLE only requires the flow map from the initial
location at the initial time to the final arrival location at the final time.
It disregards the intermediate trajectory and focuses solely on the particle’s
initial and final positions. As a result, if an initial patch of particles
remains close after a finite period, the FTLE disregards any significant
dispersion during the intermediate time and returns a small value. This
limitation is unsatisfactory because even if the trajectories of these
particles differ significantly, the FTLE fails to capture such unique flow
structures.

It is worth mentioning that while the formulation of FTLE drops the
information of the intermediate time steps, it is indeed possible to address
this issue by extending the flowmap with an extra time dimension. This
approach, as proposed in [47], allows for the preservation of the intermediate
time information. However, it is important to note that solving the resulting
partial differential equation (PDE) in the extended space can be
computationally expensive in terms of both time complexity and space
complexity. The high time complexity arises from solving the PDE, while the
space complexity increases due to the need to store and manipulate the series
of flowmaps. In contrast, our method provides an alternative approach that
avoids directly solving the flowmap while still offering a formulation to
summarize the information obtained from the intermediate time steps. This
alternative approach can be more computationally efficient compared to
directly solving the extended flowmap PDE.

##  3 Local Trajectory Variation Exponent (LTVE)

In this section, we formulate our proposed quantity. Since we consider the
variations among various local trajectories, we name the quantity the Local
Trajectory Variation (LTV) and LTVE for the corresponding exponent.

###  3.1 Theoretical Framework

We define the Local Trajectory Variation (LTV) and the LTV exponent (LTVE)
using a similar approach as the FTLE. Consider the metric space $(\Phi,d)$,
where $\Phi$ is the space of all trajectories and $d$ is the trajectory
metric. Let $\Gamma:\Omega\mapsto\Phi$ be the mapping that assigns each
particle to its trajectory. Therefore,
$\Gamma(\mathbf{x}_{0})=\Psi_{\mathbf{x}_{0}}$. To capture only the intrinsic
shape difference in trajectories, disregarding the effect of physical distance
due to variations in initial positions, we require that the mapping $\Gamma$
instead maps to the displacement trajectory. Following a similar approach as
the FTLE, we perturb the input to $\Gamma$ and measure the perturbation effect
using the trajectory metric
$d\left(\Gamma(\mathbf{x}_{0}+\delta),\Gamma(\mathbf{x}_{0})\right)$.
Formally:

######  Definition 3.1 (Local Trajectory Variation).

Given a trajectory field $\Gamma:\Omega\mapsto\Phi$ and a small
$\delta\in\mathbb{R}^{+}$, denote
$\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})=\\{\mathbf{u}\in\mathbb{R}^{n}|\left|\left|\mathbf{u}\right|\right|=1\land\mathbf{x}+\delta\mathbf{u}\in\Omega\\}$,
we define the LTV at any $\mathbf{x}\in\Omega$ as

$\text{LTV}(\mathbf{x};\delta,T)=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\;d(\Gamma(\mathbf{x}+\delta\mathbf{u}),\Gamma(\mathbf{x}))\,.$

######  Definition 3.2 (Local Trajectory Variation Exponent).

Given $\delta\in\mathbb{R}^{+}$ and the LTV function
LTV$(\mathbf{x};\delta,T)$, the LTVE$:\Omega\mapsto\mathbb{R}$ is defined as

$\text{LTVE}(\mathbf{x};\delta,T)=\frac{1}{\left|T\right|}\ln\left(\frac{1}{\delta}\text{LTV}(\mathbf{x};\delta,T)\right)$

###  3.2 Relation with FTLE

Intuitively, the flow map introduced in the FTLE can be viewed as a
$2$-trajectory, a trajectory of length $2$, and thus the FTLE can be viewed as
a special case of LTVE. Below, we provide a rigorous linkage between FTLE and
LTVE. Let $d$ be the $l^{2}$-norm distance, $\Gamma$ be the discrete
$2$-trajectory mapping, $\sigma(\mathbf{x})$ be the FTLE field, and $f$ be the
flow map induced in the construction of FTLE. Let
$\Lambda_{f}:\Omega\mapsto\mathbb{R}$ be the following function:

$\Lambda_{f}(\mathbf{x}):=\max_{\lambda}\\{\lambda\in\mathbb{R}|\exists\mathbf{u}\in\mathbb{R}^{n}\setminus\\{\mathbf{0}\\},\Delta^{T}_{0}(\mathbf{x})\mathbf{u}=\lambda\mathbf{u}\\}\,.$

Then we start by showing the following lemma.

######  Lemma 1.

If $f\in C^{1}$ on $\Omega$, then $\Lambda_{f}$ is continuous on $\Omega$.

###### Proof.

To begin with, we note that the function $\Lambda_{f}$ is well-defined since
the deformation tensor $\Delta^{T}_{0}$ is symmetric. Now, we have

$\displaystyle\Lambda_{f}(\mathbf{x}+\mathbf{u})-\Lambda_{f}(\mathbf{x})$
$\displaystyle=\max_{\mathbf{v}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|\nabla
f(\mathbf{x}+\mathbf{u})\cdot\mathbf{v}\right|\right|^{2}-\max_{\mathbf{v}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|\nabla
f(\mathbf{x})\cdot\mathbf{v}\right|\right|^{2}$
$\displaystyle\leq\max_{\mathbf{v}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left\\{\left|\left|\nabla
f(\mathbf{x}+\mathbf{u})\cdot\mathbf{v}\right|\right|^{2}-\left|\left|\nabla
f(\mathbf{x})\cdot\mathbf{v}\right|\right|^{2}\right\\}$
$\displaystyle=\max_{\mathbf{v}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left\\{(\left|\left|\nabla
f(\mathbf{x}+\mathbf{u})\cdot\mathbf{v}\right|\right|+\left|\left|\nabla
f(\mathbf{x})\cdot\mathbf{v}\right|\right|)(\left|\left|\nabla
f(\mathbf{x}+\mathbf{u})\cdot\mathbf{v}\right|\right|-\left|\left|\nabla
f(\mathbf{x})\cdot\mathbf{v}\right|\right|)\right\\}\,.$

Since $\Omega$ is closed and bounded, and $f\in C^{1}$ implies
$\left|\left|\nabla f(\mathbf{x})\right|\right|<L$ for some $L\in\mathbb{R}$,
we have

$\Lambda_{f}(\mathbf{x}+\mathbf{u})-\Lambda_{f}(\mathbf{x})\leq
2L\max_{\mathbf{v}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|\nabla
f(\mathbf{x}+\mathbf{u})\cdot\mathbf{v}-\nabla
f(\mathbf{x})\cdot\mathbf{v}\right|\right|\leq 2L\left|\left|\nabla
f(\mathbf{x}+\mathbf{u})-\nabla f(\mathbf{x})\right|\right|\,.$

Similarly, we can obtain another inequality
$\Lambda_{f}(\mathbf{x})-\Lambda_{f}(\mathbf{x}+\mathbf{u})\leq
2L\left|\left|\nabla f(\mathbf{x}+\mathbf{u})-\nabla
f(\mathbf{x})\right|\right|$. These two inequalities imply that
$\left|\Lambda_{f}(\mathbf{x}+\mathbf{u})-\Lambda_{f}(\mathbf{x})\right|\leq
2L\left|\left|\nabla f(\mathbf{x}+\mathbf{u})-\nabla
f(\mathbf{x})\right|\right|$. Then it follows that $f\in C^{1}$ implies the
latter term can be bounded arbitrarily small. ∎

From the definition of FTLE, we can express the FTLE field using the greatest
eigenvalue $\lambda_{\max}(\mathbf{x})$ of $\Delta_{0}^{T}(\mathbf{x})$. In
particular, we have
$\sigma(\mathbf{x})=\frac{1}{2\left|T\right|}\ln\lambda_{\max}$. Hence, for
the FTLE value to be well-defined over $\Omega$, we require $\lambda_{\max}>0$
at all $\mathbf{x}\in\Omega$. Then, it follows from Lemma 1 that
$\lambda_{*}:=\min_{\mathbf{x}\in\Omega}\Lambda_{f}(\mathbf{x})$ is well-
defined and satisfies $0\leq\lambda_{*}\leq\lambda_{\max}(\mathbf{x})$ for all
$\mathbf{x}\in\Omega$. Since we assumed the mapping maps to the displacement
trajectory $\Psi$, any $2$-trajectory evaluated under such a norm will give
the particle displacement at time $t=t_{0}+T$. Namely,

$\text{LTV}(\mathbf{x})=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\;d(\Gamma(\mathbf{x}+\delta\mathbf{u}),\Gamma(\mathbf{x}))=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})-\delta\mathbf{u}\right|\right|_{2}\,.$

For simplicity, we denote
$\mathcal{F}_{\delta}(\mathbf{x}):=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left\\{\left|\left|f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})\right|\right|_{2}\right\\}$.
Recall from the construction of FTLE, we have
$\mathcal{F}_{\delta}(\mathbf{x})\approx\sqrt{\lambda_{\max}}\delta\geq\sqrt{\lambda_{*}}\delta$.
Thus, we can approximate the FTLE as
$\sigma(\mathbf{x})\approx\frac{1}{\left|T\right|}\ln\frac{\mathcal{F}_{\delta}(\mathbf{x})}{\delta}$.
Then consider

$\displaystyle\text{LTV}(\mathbf{x})^{2}$
$\displaystyle=\left[\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})-\delta\mathbf{u}\right|\right|_{2}\right]^{2}=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})-\delta\mathbf{u}\right|\right|_{2}^{2}$
$\displaystyle=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left\\{\left|\left|f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})\right|\right|_{2}^{2}+\delta^{2}\left|\left|\mathbf{u}\right|\right|^{2}_{2}-2\delta\langle
f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x}),\mathbf{u}\rangle\right\\}$
$\displaystyle=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left\\{\left|\left|f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})\right|\right|_{2}\left(\left|\left|f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})\right|\right|_{2}-2\delta\left|\left|\mathbf{u}\right|\right|_{2}\cos\theta\right)\right\\}+\delta^{2}$

where $\theta$ is determined by the angle between
$f(\mathbf{x}+\delta\mathbf{u})-f(\mathbf{x})$ and $\mathbf{u}$. Taking the
trivial bound of $\cos\theta\geq-1$, we obtain
$\text{LTV}(\mathbf{x})^{2}\leq\mathcal{F}_{\delta}(\mathbf{x})\left(\mathcal{F}_{\delta}(\mathbf{x})+2\delta\right)+\delta^{2}$,
and therefore,

$\frac{\text{LTV}(\mathbf{x})^{2}}{\mathcal{F}_{\delta}(\mathbf{x})^{2}}\leq
1+\frac{\delta^{2}}{\mathcal{F}_{\delta}(\mathbf{x})^{2}}+2\frac{\delta}{\mathcal{F}_{\delta}(\mathbf{x})}\implies\frac{\text{LTV}(\mathbf{x})}{\mathcal{F}_{\delta}(\mathbf{x})}\leq
1+\frac{\delta}{\mathcal{F}_{\delta}(\mathbf{x})}\,.$

Since the logarithmic function is increasing, we have

$\ln\frac{\text{LTV}(\mathbf{x})}{\mathcal{F}_{\delta}(\mathbf{x})}\leq\ln\left(1+\frac{\delta}{\mathcal{F}_{\delta}(\mathbf{x})}\right)\leq\ln\left(1+\frac{1}{\sqrt{\lambda_{*}}}\right)\,.$

On the other hand, by considering $g(\mathbf{x}):=f(\mathbf{x})-\mathbf{x}$,
we can rewrite LTV in a different form,
$\text{LTV}(\mathbf{x})=\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|g(\mathbf{x}+\delta\mathbf{u})-g(\mathbf{x})\right|\right|\approx\max_{\mathbf{u}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\left|\left|\nabla
g\cdot\delta\mathbf{u}\right|\right|$. Then since $f\in C^{1}$, we have $g\in
C^{1}$. Hence following the same approach, we can define
$\eta_{*}:=\min_{\mathbf{x}\in\Omega}\Lambda_{g}(\mathbf{x})$ and thus
$\text{LTV}(\mathbf{x})\geq\sqrt{\eta_{*}}\delta$ for all
$\mathbf{x}\in\Omega$. Hence we have
$\ln\left[{\mathcal{F}_{\delta}(\mathbf{x})}/\text{LTV}(\mathbf{x})\right]\leq\ln\left(1+{1}/{\sqrt{\eta_{*}}}\right)$.
Thus it follows from

$\left|\text{LTVE}(\mathbf{x})-\sigma(\mathbf{x})\right|=\frac{1}{\left|T\right|}\left|\ln\left(\frac{\text{LTV}(\mathbf{x})}{\delta}\right)-\ln\left(\frac{\mathcal{F}_{\delta}(\mathbf{x})}{\delta}\right)\right|=\frac{1}{\left|T\right|}\left|\ln\frac{\text{LTV}(\mathbf{x})}{\mathcal{F}_{\delta}(\mathbf{x})}\right|\,,$

we arrive our final conclusion that

$\left|\text{LTVE}(\mathbf{x})-\sigma(\mathbf{x})\right|\leq\frac{1}{\left|T\right|}\max\left\\{\ln\left(1+\frac{1}{\sqrt{\eta_{*}}}\right),\ln\left(1+\frac{1}{\sqrt{\lambda_{*}}}\right)\right\\}$
(3)

for all $\mathbf{x}\in\Omega$.

The introduction of $\lambda_{*}$ and $\eta_{*}$ in equation (3) aims to
establish a universal bound with respect to the domain $\Omega$. By
definition, $\lambda_{*}$ and $\eta_{*}$ represent the lower bounds of
$\Lambda_{f}$ and $\Lambda_{g}$, respectively. This implies that for each
point in the domain $\Omega$, there exists a direction of perturbation that
induces a change in the trajectory. Considering this condition is natural, as
the FTLE should be well-defined at a point when the maximal eigenvalue of the
deformation tensor is positive, indicating trajectory deviation when
perturbed. Furthermore, allowing the bounds to depend on $\mathbf{x}$, the
location in the domain, can yield tighter bounds. Intuitively, if we evaluate
the bound near the LCS, the corresponding maximal eigenvalue should be large,
resulting in a smaller bound between the LTVE and FTLE quantities. Ideally,
our quantity can capture all structures that the FTLE captures by applying the
$l^{2}$-norm with two trajectories.

Instead of solely considering the takeoff and arrival locations of a
trajectory, we have expanded our investigation by adding a third point to each
trajectory data. This addition allows us to explore longer trajectories and
their impact on our quantity. In Section 4.3, we will present an example that
demonstrates the superiority of our quantity over the FTLE when considering
these longer trajectories. This example will showcase a significant difference
between the LTVE with three trajectories and the FTLE, emphasizing the
advantages of our approach.

###  3.3 Numerical Scheme and Error Analysis

To perform the numerical computation of the quantity, it is necessary to
discretize and apply relaxation techniques to make the problem computationally
feasible. Below, we provide a brief summary of these processes and analyze
their associated errors.

First, we focus on the discretization in the time domain. We employ a simple
scheme to approximate the real trajectory. We uniformly sample $k$ points
along the trajectory, where $k$ can be adjusted either adaptively or manually
depending on the problem. The trajectory is then approximated by a polygonal
line connecting these sampling points, which essentially provides a linear
approximation. While this approximation may introduce imprecise results when
$k$ is poorly chosen, opting for a higher-order approximation would increase
the complexity of designing and implementing the trajectory comparison metric.

Now, we analyze the error introduced by this approximation. For any continuous
$\Psi$, we denote $\Psi_{\Delta}$ as the discrete trajectory approximated by
the polygonal line. Similarly, for the continuous trajectory metric $d$, we
denote $d_{\Delta}$ as its discrete counterpart. For any (continuous)
trajectories $\Psi_{1}$ and $\Psi_{2}$, using metric properties, we have:

$d(\Psi_{1},\Psi_{2})\leq
d(\Psi_{1},\Psi_{1,\Delta})+d(\Psi_{1,\Delta},\Psi_{2,\Delta})+d(\Psi_{2},\Psi_{2,\Delta})$

and

$d(\Psi_{1,\Delta},\Psi_{2,\Delta})\leq
d(\Psi_{1},\Psi_{1,\Delta})+d(\Psi_{1},\Psi_{2})+d(\Psi_{2},\Psi_{2,\Delta})\,.$

Combining them, we have
$\left|d(\Psi_{1},\Psi_{2})-d(\Psi_{1,\Delta},\Psi_{2,\Delta})\right|\leq
d(\Psi_{1},\Psi_{1,\Delta})+d(\Psi_{2},\Psi_{2,\Delta})$. Hence the error can
be separated into two parts,

$\displaystyle\left|d(\Psi_{1},\Psi_{2})-d_{\Delta}(\Psi_{1,\Delta},\Psi_{2,\Delta})\right|$
$\displaystyle\leq\left|d(\Psi_{1},\Psi_{2})-d(\Psi_{1,\Delta},\Psi_{2,\Delta})\right|+\left|d(\Psi_{1,\Delta},\Psi_{2,\Delta})-d_{\Delta}(\Psi_{1,\Delta},\Psi_{2,\Delta})\right|$
$\displaystyle\leq
d(\Psi_{1},\Psi_{1,\Delta})+d(\Psi_{2},\Psi_{2,\Delta})+\left|d(\Psi_{1,\Delta},\Psi_{2,\Delta})-d_{\Delta}(\Psi_{1,\Delta},\Psi_{2,\Delta})\right|\,.$

Taking the supremum over the space of trajectories, we have the error bound

$2\sup_{\Psi}\\{d(\Psi_{1},\Psi_{1,\Delta})\\}+\sup_{\Psi_{1,\Delta},\Psi_{2,\Delta}}\left|d(\Psi_{1,\Delta},\Psi_{2,\Delta})-d_{\Delta}(\Psi_{1,\Delta},\Psi_{2,\Delta})\right|\,.$

The first part addresses the error resulting from estimating the trajectory
using a polygonal line under the desired metric. As mentioned earlier, this
error arises from the use of linear approximation for trajectories. One can
reduce this term by increasing the number of sampling points, adjusting the
parameter $k$, or employing a higher-order approximation scheme to achieve a
more accurate solution. The second part of the error analysis stems from the
discretization of the metric itself. Approximating the metric using its
discrete variant introduces additional error. For metrics such as the
$L^{p}$-norm, one can reduce this term to arbitrary precision by employing a
higher order of numerical integration. Explicitly considering the
discretization of the metric is important, especially for metrics like the
Fréchet distance and Hausdorff metric. The computational complexity of these
metrics in the continuous setting is prohibitively high for practical
implementation. Therefore, including an approximation of the metric itself
becomes necessary. The actual order of the error depends on the chosen metric
and underlying field, as mentioned previously. Providing a general bound on
the error without specifying the metric is challenging. In our chosen metric
for the experiment, all these errors are well-bounded. The derivation of these
bounds is either standard or well-studied in the literature, and we omit it
here for brevity.

Next, we focus on the relaxation in the space domain. It is important to note
that solving for the optimal $\mathbf{u}$ can be computationally expensive.
However, we can adopt an approach similar to the FTLE by constructing the
Cauchy-Green deformation tensor for $T$. It should be noted that this process
is more complex compared to the FTLE case, given the higher dimensionality
involved. Nevertheless, as we will see in later numerical experiments, our
relaxation scheme is sufficient for extracting the LCS.

Our relaxation scheme is straightforward. Instead of comparing the
trajectories within the $\delta$-unit circle, we only compare them to nearby
trajectories. The specific choice of neighborhood depends on the scheme
specified by the users. Below, we propose several possible schemes for
selecting the neighborhood:

######  Definition 3.3 (First Order Relaxed-LTV).

Given a trajectory field $\Gamma:\Omega\mapsto\Phi$, fix some small
$\delta\in\mathbb{R}^{+}$, denote $e_{j}$ to be the $j$-th standard basis in
$\mathbb{R}^{n}$, then $\forall\mathbf{x}\in\Omega$, we define

$\text{RLTV}_{1}(\mathbf{x})=\max_{e_{j}\in\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}\max\left\\{d(\Gamma(\mathbf{x}+\delta
e_{j}),\Gamma(\mathbf{x})),d(\Gamma(\mathbf{x}-\delta
e_{j}),\Gamma(\mathbf{x}))\right\\}\,.$

######  Definition 3.4 (Second Order Relaxed-LTV).

Given a trajectory field $\Gamma:\Omega\mapsto\Phi$, fix some small
$\delta\in\mathbb{R}^{+}$, denote $e_{j}$ to be the $j$-th standard basis in
$\mathbb{R}^{n}$, then denote
$\mathcal{N}:=\\{\mathbf{u}\in\mathbb{R}^{n}|\mathbf{u}=\pm a_{1}e_{1}\pm
a_{2}e_{2}\pm...\pm a_{n}e_{n},\,a_{i}=\\{0,1\\}\enspace\forall i\\}$. Then
$\forall\mathbf{x}\in\Omega$, we define

$\text{RLTV}_{2}(\mathbf{x})=\max_{\mathbf{u}\in\mathcal{N}\cap\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}d(\Gamma(\mathbf{x}+\delta\mathbf{u}),\Gamma(\mathbf{x}))\,.$

######  Definition 3.5 (Third Order Relaxed-LTV).

Given a trajectory field $\Gamma:\Omega\mapsto\Phi$, fix some small
$\delta\in\mathbb{R}^{+}$, denote $e_{j}$ to be the $j$-th standard basis in
$\mathbb{R}^{n}$, then denote
$\mathcal{N}:=\\{\mathbf{u}\in\mathbb{R}^{n}|\mathbf{u}=\pm a_{1}e_{1}\pm
a_{2}e_{2}\pm...\pm a_{n}e_{n},\,a_{i}=\\{0,1\\}\enspace\forall i\\}$ and
$\mathcal{N^{\prime}}:=\\{\mathbf{v}\in\mathbb{R}^{n}|\mathbf{v}=\mathbf{u}+e_{i},\,\forall\mathbf{u}\in\mathcal{N},\,\forall
i\\}\setminus\\{2e_{i}\in\mathbb{R}^{n}|\forall i\\}\,.$ Then
$\forall\mathbf{x}\in\Omega$, we define

$\text{RLTV}_{3}(\mathbf{x})=\max_{\mathbf{u}\in\mathcal{N^{\prime}}\cap\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}d(\Gamma(\mathbf{x}+\delta\mathbf{u}),\Gamma(\mathbf{x}))\,.$

The corresponding RLTVE are defined similarly as in the LTVE case and we omit
it here for brevity. While it is indeed possible to adopt a random sampling
approach, such as Monte-Carlo sampling, to approximate the optimal
$\mathbf{u}$ without using relaxation, we have chosen to utilize the
relaxation scheme proposed for improving time complexity. The main advantage
of the relaxation scheme is the ability to reuse trajectories. Since all the
comparison points lie on the mesh points, it is sufficient to compute the
trajectories only at the mesh points. This significantly reduces the
computational complexity. If we were to use random sampling, the complexity
would scale linearly with the sampling size. By utilizing the relaxation
scheme and working with the mesh points, we can achieve better computational
efficiency compared to random sampling approaches.

###  3.4 Algorithm

Below we present the main algorithm for computing the RLTVE field and analyze
its time and space complexity. Given a rectangular computation domain
$\Omega=[x_{1,0},x_{1,1}]\times...\times[x_{n,0},x_{n,1}]$ and the Lipchitz
velocity field $\mathbf{v}:\Omega\mapsto\mathbb{R}^{n}$, fix some small
$\delta\in\mathbb{R}^{+}$ for the relaxation condition, then we adapt a
uniform mesh such that each particles are separated with distance $\delta$
along each standard axis direction.

Function  _RLTVE(_$\mathbf{v}:\Omega\mapsto\mathbb{R}^{n}$
,$\delta\in\mathbb{R}^{+}$_)_:

$M\leftarrow$ uniform meshes with $\Delta x_{j}=\delta$ for all $j\in[1,d]$

$\Psi\leftarrow$ Solve the trajectory field from $\mathbf{v}$ using ODE solver
with $M$ as meshes

$\Psi\leftarrow$ transform to displacement trajectory

$RLTVE\leftarrow$ Int array of equal dimension as $M$ with initial value $0$

for  _$\mathbf{x}\in M$_ do

$\mathcal{N}\leftarrow$ neighbourhood of $\mathbf{x}$ defined by the chosen
scheme.

$RLTVE[\mathbf{x}]\leftarrow\frac{1}{\left|T\right|}\ln\frac{1}{\delta}\max_{\mathbf{u}\in\mathcal{N}\cap\mathcal{B}^{\Omega}_{\delta}(\mathbf{x})}d(T(\mathbf{x}+\delta\mathbf{u}),T(\mathbf{x}))$

end for

return _$RLTVE$_

end

Algorithm 1 RLTVE

###  3.5 Implementation Details and Analysis

In this section, we summarize some technical implementation details and
provide some analysis of the algorithm’s performance. We will consider both
the computational time complexity and the algorithm’s memory requirement.

##### Boundary Condition.

Upon solving the trajectories using the given velocity field, it is possible
for particles to exit the computational domain $\Omega$. In cases where no
prior information is available outside $\Omega$, we need to impose artificial
boundary conditions on the numerical integrator. One approach, as suggested in
[34], is to implement a velocity extension algorithm. This method creates an
artificial layer that gradually stops the particles, similar to the absorbing
boundary condition used in solving the wave equation. However, in our work,
instead of exploring the effects of the velocity extension, we directly
enforce a stopping boundary condition. When the numerical solution from the
ODE leaves the computational domain, we halt the numerical integrator and
project the solution onto the boundary of the computational domain
$\partial\Omega$. We will provide more details on this effect in Section 4.7.

##### Time Complexity.

Now let us analyze the time complexity. First, we note that using uniform
meshes $M$ is not necessary for actual implementation. All particle positions
can be computed using looping indices with a given $\delta$. Therefore, the
time complexity mainly comes from the ODE solver and the computation of RLTVE.
However, it is possible to pre-compute and store trajectory data, reading it
from a file during computation. Still, it takes
$O\left({\left|\Omega\right|}{\delta^{-n}}k\right)$ time to write all the
values for each trajectory in the domain. Next, for computing the RLTVE,
although depending on the scheme of relaxation chosen, the comparison done for
each particle can be approximated as $O(n)$. It is important to note that the
complexity of each computation depends on the metric $d$, and it is generally
not negligible as the number of sampling points in each trajectory increases.
Let us denote the time complexity for computing a single query of the metric
as $T_{d}(n,k)$. The total time complexity for computing the RLTVE would then
be $O\left({\left|\Omega\right|}{\delta^{-n}}nT_{d}(n,k)\right)$. It is worth
mentioning that for the trajectory metric applied in this project, it takes at
least $\Omega(k)$ time to compute. Furthermore, since comparing two arbitrary
trajectories without any special constraints requires going through each
sampling point at least once, we can assume that $T_{d}(n,k)\in\Omega(k)$.
Hence the complexity is at least
$O\left({\left|\Omega\right|}{\delta^{-n}}nk\right)$ but could be potentially
higher. In summary, the total time complexity for general metric is given by
$O\left({\left|\Omega\right|}{\delta^{-n}}nT_{d}(n,k)\right)$.

##### Storage.

Next, we consider the memory required for the algorithm. A naive approach of
saving all trajectories and the computed RLTVE values at each point would take
space of $O\left({\left|\Omega\right|}{\delta^{-n}}(k+1)\right)$. However, we
can reduce it to
$O\left({\left|\Omega^{\prime}\right|}{\delta^{-(n-1)}}k\right)$ without
affecting the time complexity, where
$\left|\Omega^{\prime}\right|={\left|\Omega\right|}/{\max_{j\in[1,n]}\left|x_{j,1}-x_{j,0}\right|}$.
The reduction in space usage is based on the observation that our calculations
are based on comparisons with the neighborhood. Therefore, we can discard
trajectory data that is too far from the current computing position. To
achieve this, we first select an axis along which $\Omega$ has its longest
side. We then compute and store all the trajectories on the co-dimension one
faces that are tangential to the chosen longest axis and touch one of the
boundaries of the corresponding side. Next, we move to the next set of faces,
shifting $\delta$ from the boundary of the longest side. Using the previously
calculated values, we can reduce the number of calculations required in these
faces by $1$. After computing all the faces, we delete the stored values from
the previous layer and replace them with the values calculated in the current
layer. With this approach, we only need to store at most
$O\left({\left|\Omega^{\prime}\right|}{\delta^{-(n-1)}}(k+1)\right)$, which is
on the surface. Furthermore, for the RLTVE values, they can be computed and
written to the output file on-the-fly during the trajectory computation. This
further reduces our memory complexity to
$O\left({\left|\Omega^{\prime}\right|}{\delta^{-(n-1)}}k\right)$.

###  3.6 Trajectory Metric

In this section, we specify our choice of metric. In this paper, we
investigate three different trajectory metrics, and we are going to give a
brief definition of each below. In the following discussion, $\Phi$ denotes
the time-continuous trajectory space, and $\Phi_{\Delta}$ represents the space
of all $k$-discrete trajectory.

####  3.6.1 $L^{2}$-norm

As might be expected that the most natural choice for trajectory metric will
be the $L^{2}$-norm, it is a well-studied norm in classic analysis and serve
as a good intuitive for understanding the physical meaning of the value
captured by our quantity. Moreover it is easy to be computed, in particular
its discrete variant, $l^{2}$ norm, has a linear complexity and is easy to
implement.

######  Definition 3.6 (Normalized $L^{2}$-norm).

The normalized $L^{2}$-norm
$\left|\left|\cdot\right|\right|_{\overline{L}^{2}}:\Phi\mapsto\mathbb{R}$ is
defined as
$\left|\left|\Psi\right|\right|_{\overline{L}^{2}}^{2}=\frac{1}{\left|T\right|}\int_{t_{0}}^{t_{0}+T}\left|\Psi(t)\right|^{2}\,dt$.
For the discrete variant, we naturally take the discrete analog of
$L^{2}$-norm, we define the normalized $l^{2}$-norm as
$\left|\left|\cdot\right|\right|_{\overline{l^{2}}}:\Phi_{\Delta}\mapsto\mathbb{R}$,
$\left|\left|\Psi_{\Delta}\right|\right|_{\overline{l^{2}}}^{2}=\frac{1}{k}\left|\left|\Psi_{\Delta}\right|\right|_{l^{2}}^{2}$
where $\left|\left|\cdot\right|\right|_{l^{2}}$ denote the conventional
$l^{2}$-norm on $\mathbb{R}^{kn}$, such definition is safe by recalling the
isomorphism of $\mathbb{R}^{kn}$ and $\Phi_{\Delta}$ as mentioned in
def_discrete.

####  3.6.2 Fréchet Distance

Fréchet distance is a quantity that captures the minimal maximum separation of
trajectories over reparameterization in the time domain. It is a popular
similarity measure of trajectories used in the field of trajectory analysis
and clustering. This distance measure considers the reparameterization of the
trajectory in the time-domain, which makes it interesting to study, as such a
property might be able to capture structures that the classical method fails
to identify.

######  Definition 3.7 (Fréchet Distance).

Denote $\mathcal{P}$ as the space of function
$f:[t_{0},t_{0}+T]\mapsto[t_{0},t_{0}+T]$ which satisfies $\forall
f\in\mathcal{P}$, $f$ is continuous, surjective and is non-decreasing. Then
the Fréchet distance $d_{F}:\Phi\times\Phi\mapsto\mathbb{R}$ is defined as

$d_{f}(\Psi_{1},\Psi_{2})=\inf_{\alpha,\beta\in\mathcal{P}}\max_{\tau\in[t_{0},t_{0}+T]}\left|\left|\Psi_{1}(\alpha(\tau))-\Psi_{2}(\beta(\tau))\right|\right|_{2}$

where $\left|\left|\cdot\right|\right|_{2}$ is the $l^{2}$ norm induced by the
space $\mathbb{R}^{n}$.

Unfortunately, there are certain issues with the Fréchet distance. The first
issue is that it does not meet the requirement of being a metric.
Specifically, it is only a pseudo-metric, which means it fails to meet the
condition that $d(a,b)=0\iff a=b$. However, we can relax this condition and
still use the Fréchet distance by using its discrete variant. The work [7]
defined the discrete variant of the Fréchet distance and showed that it
preserves the metric property. It also serves as a good approximation to the
Fréchet distance with respect to the sampling size, resolving the first issue.
The second issue is the time complexity. The discrete Fréchet distance can be
computed using dynamic programming, which has a runtime of $O(nm)$, where $n$
and $m$ are the lengths of the two trajectories. In [3], it is proven that for
both continuous and discrete Fréchet distance, there is no strong subquadratic
algorithm. This means that there is no algorithm with a complexity of
$O(n^{2-\delta})$ for some $\delta>0$, unless the strong exponential time
hypothesis (SETH) fails. SETH is a long-standing hypothesis in the study of
algorithmic complexity theory. It states that there is no algorithm that can
solve the $k$-SAT problem with length $n$ in $O((2-\delta)^{n})$ time
$\forall\delta>0$. To address this issue, there are several algorithms that
have been proposed. The state-of-the-art solutions include an exact algorithm
proposed by [1] that runs in weak subquadratic time, as well as an almost-
linear time approximation algorithm proposed by [5]. These algorithms provide
efficient ways to compute the Fréchet distance. In this project, we have
implemented the $O(k^{2})$ solution using dynamic programming since the
efficiency is not the main focus of studying this metric, we are more
interested in whether it produces different results from the $L^{2}$-norm.

######  Definition 3.8 (Discrete Fréchet Distance).

For all $\Psi_{\Delta,1},\Psi_{\Delta,2}\in\Phi_{\Delta}$, the coupling of
$\Psi_{\Delta,1},\Psi_{\Delta,2}$ is defined as

$L(\Psi_{\Delta,1},\Psi_{\Delta,2}):=(\Psi_{\Delta,1}(a_{0}),\Psi_{\Delta,2}(a_{0})),(\Psi_{\Delta,1}(a_{1}),\Psi_{\Delta,2}(b_{1})),...,(\Psi_{\Delta,1}(a_{k-1}),\Psi_{\Delta,2}(a_{k-1}))$

and satisfying $a_{k-1}=b_{k-1}=k-1$ with $a_{0}=b_{0}=0$ and
$(a_{j+1}=a_{j}+1\land b_{j}=b_{j+1})\lor(b_{j+1}=b_{j}+1\land
a_{j}=a_{j+1})\quad\forall j\in[0,k-2]$. Then for each coupling, its length
$\left|\left|L\right|\right|$ is defined as
$\left|\left|L\right|\right|=\max_{i\in[0,n]}\left|\left|\Psi_{\Delta,1}(a_{0})-\Psi_{\Delta,2}(a_{0})\right|\right|_{2}$
then the discrete Fréchet distance $d_{F_{\Delta}}$ is defined as
$d_{F_{\Delta}}(\Psi_{\Delta,1},\Psi_{\Delta,2})=\min\\{\left|\left|L\right|\right|\in\mathbb{R}|L\text{
is a coupling of $\Psi_{\Delta,1}$ and $\Psi_{\Delta,2}$}\\}$.

####  3.6.3 Hausdorff Metric

The last metric we have chosen is the Hausdorff metric, which is another
important metric for comparing the similarity of two different trajectories or
sets of points. Unlike other metrics, the Hausdorff metric does not consider
the ordering of points in the set. It is worth noting that we are not the
first to apply the Hausdorff metric for computing the LCS. In previous works
[37, 2], they also incorporated the Hausdorff metric into their frameworks for
computing the LCS. The Hausdorff metric can be naively computed in $O(k^{2})$
time complexity. However, similar to the previous metric, our main focus is
not on the runtime of this metric but rather on its difference from the
$L^{2}$-norm.

######  Definition 3.9 (Hausdorff Metric).

The one sided Hausdorff metric
$d_{H}^{\prime}:\Phi\times\mathbb{R}^{n}\mapsto\mathbb{R}$ is defined as
$d_{H}^{\prime}(\Psi,\psi)=\inf_{\tau\in[t_{0},t_{0}+T]}\left|\left|\Psi(\tau)-\psi\right|\right|_{2}$
where $\left|\left|\cdot\right|\right|_{2}$ is the $l^{2}$ norm induced by the
space $\mathbb{R}^{n}$.
Then the Hausdorff metric is defined as

$d_{H}(\Psi_{1},\Psi_{2})=\max\left\\{\sup_{t\in[t_{0},t_{0}+T]}d_{H}^{\prime}(\Psi_{1},\Psi_{2}(t)),\sup_{t\in[t_{0},t_{0}+T]}d_{H}^{\prime}(\Psi_{2},\Psi_{1}(t))\right\\}\,.$

Noted that by simply changing the domain of $\sup$ and $\inf$ correspondingly
to $[0,k-1]$, the above definition works for $k$-discrete trajectory as well.

##  4 Numerical Examples

This paper specifically focuses on two-dimensional flows. Extending the
algorithm to three-dimensional examples is relatively straightforward, but we
will not cover that in this paper. The trajectories in all the test cases have
been pre-computed using the fourth-order Runge-Kutta method, with the
underlying velocity known analytically. For the implementation of the main
computation algorithm, we have used the programming language C++.
Additionally, for visualization purposes, we have utilized the contourf
function in MATLAB for simplicity. In all of the experiments, unless otherwise
specified, we used 100 sampling points for the trajectories.

###  4.1 Circular Flow

(a) (b) (c) (d)

Figure 1: (Section 4.1) The computed LTVE using (a) the normalized 2-norm, (b)
the Fréchet distance, and (c) the Hausdorff metric. (d) The corresponding FTLE
field.

We begin with a simple toy example that intuitively demonstrates a structure
in the field. The field is designed with a small ring-shaped flow on a 2D
plane, centered at the origin. At the boundary of the ring, a notable
difference in velocity creates a distinct structure. The velocity field is
given by:

$\frac{d\mathbf{x}}{dt}=\begin{cases}(y,-x)^{T}\quad&(1-\varepsilon)^{2}\leq(x+y)^{2}\leq(1+\varepsilon)^{2}\\\
(0,0)^{T}\quad&\text{otherwise}\end{cases}$

where $\varepsilon$ is some small constant, in our experiment it is taken to
be $0.1$, we run the experiment with domain $[-2,2]\times[-2,2]$ and
$\delta={4}/{399}$ and time range $[0,1]$.

As we examine the visualizations in Figure 1 (a-d), we can observe that all
the visualizations successfully highlight the structure we expected, which
corresponds to the boundaries of the ring. It is important to note that in
Figure 1 (a-c), there is a large blank area that may initially appear as a
visualization issue. However, this blank area corresponds to regions where the
field remains static. Consequently, the LTVE at those points is undefined
because the exact LTV is zero. The same theoretical result holds for the FTLE
case, where the LTVE is expected to be undefined in static regions. However,
in practice, the computation of numerical gradients involved in the FTLE
calculation may introduce computational errors, preventing exact zero values
from being obtained. This discrepancy explains the visual differences between
the LTVE and FTLE visualizations. It is worth noting that this issue only
arises when the field is static and does not affect the visualization of LCS.

###  4.2 Standing Wave

(a) (b) (c) (d)

Figure 2: (Section 4.2) The computed LTVE using (a) the normalized 2-norm, (b)
the Fréchet distance, and (c) the Hausdorff metric. (d) The corresponding FTLE
field.

In this example, we consider a simple flow where we can analytically determine
the trajectory of each individual particle. The FTLE only considers the final
location of the particles and disregards any information from the intermediate
steps. As a result, the FTLE field heavily depends on the chosen final
visualization time. To demonstrate the superiority of our proposed LTVE, we
introduce a velocity field given by:

$\frac{d\mathbf{x}}{dt}=\left(\frac{\omega\varepsilon\cos(\omega
t)}{1+\varepsilon\sin(\omega t)}\right)\mathbf{x}$

Here, $\varepsilon=\varepsilon(\mathbf{x}_{0})$, where $\mathbf{x}_{0}$
represents the initial location of the particle at $t=0$. We model
$\varepsilon$ as $\varepsilon=\varepsilon_{\max}\cos(k\theta)$, where
$\theta=\tan^{-1}\left(\frac{y}{x}\right)$ is the polar angle representation
of $\mathbf{x}_{0}$, and $\varepsilon_{\max}=0.8$ controls the magnitude of
the oscillations. The trajectory of each particle is given by
$\mathbf{x}(t)=(1+\varepsilon\sin(\omega t))\mathbf{x}_{0}$. By considering
this specific flow, we can demonstrate that our proposed LTVE provides a
better understanding of the underlying flow compared to the FTLE. Unlike the
FTLE, which only considers the final location of the particles, the LTVE takes
into account the information from all intermediate steps. This allows for a
more comprehensive analysis of the flow dynamics.

This specific flow exhibits two important properties. First, for
$\theta=\frac{(2n-1)\pi}{2k}$, where $n$ and $k$ are integers, we observe that
$\varepsilon=0$, resulting in stationary particles for all time. Second, for
all other particles, we observe oscillatory behavior with a period of
$2\pi/\omega$. These particles return to their original locations after this
period, leading to an FTLE value of 0 throughout the computational domain. In
this numerical example, we choose $\omega=8\pi$ and $k=7$. We compute the
trajectories for a time period of $[0,1]$. The computational domain is set to
$[-3,3]^{2}$, and we use a grid spacing of $\delta={6}/{299}$ in the
definition of the LTVE. Figure 2(d) shows the FTLE field, which evaluates to 0
with some noise on the order of $10^{-7}$. However, in Figures 2(a-c), the
structure of the flow is clearly visualized, with LTVE values ranging from
$10^{-1}$ to $10^{-2}$. This significant difference in magnitudes demonstrates
the distinguishability of our proposed quantity from the classical FTLE
result.

###  4.3 Evaluation with 3-Trajectories

(a) (b) (c) (d)

Figure 3: (Section 4.3) The computed LTVE using (a) the normalized 2-norm, (b)
the Fréchet distance, and (c) the Hausdorff metric. (d) The corresponding FTLE
field.

In this section, we recompute the LTVE for the Standing Wave example
introduced in Section 4.2 with a different setting. Specifically, we use only
3 sample points for the trajectories and choose a final time of $T=0.125$,
while keeping all other parameters unchanged.

As shown in Figure 3, the results obtained from evaluating the LTVE with only
three sampling points in the trajectories align with the findings reported in
Section 4.2. The LTVE values obtained exhibit a substantial difference,
spanning several orders of magnitude, compared to the FTLE values. This stark
contrast demonstrates that even with the addition of just one extra sampling
point to the trajectories, our LTVE quantity is capable of identifying
structures that the FTLE fails to capture.

###  4.4 Double Gyre Flow

(a) (b)
(c) (d)

Figure 4: (Section 4.4) The computed LTVE using (a) the normalized 2-norm, (b)
the Fréchet distance, and (c) the Hausdorff metric. (d) The corresponding FTLE
field.

(a) (b) (c) (d)

Figure 5: (Section 4.4) The computed LTV using (a) the normalized 2-norm and
the Fréchet distance. (b) We also plot the trajectory of several particles.
(c) The zoom-in to the red boxed region on left and (d) the red boxed region
on right in (b). Figure 6: (Section 4.4) The actual difference between the
LTVE and the FTLE at different time. The estimate from equation (3) is shown
using the dashed line.

In this section, we proceed with a common example used for computing the LCS
of dynamic systems, known as the double gyre flow. We use this example to
illustrate how our quantity is able to accurately capture the LCS. As a
convention, we choose the domain to be $[0,2]\times[0,1]$ with a spacing of
$\delta={1}/{249}$. The time range is taken to be $[0,15]$ with 100 sampling
points along the trajectory. The velocity field is given by:

$\frac{d\mathbf{x}}{dt}=\begin{pmatrix}-\pi A\sin(\pi(ax^{2}+bx))\cos(\pi
y)\\\ \pi A(2ax+b)\cos(\pi(ax^{2}+bx))\sin(\pi y)\end{pmatrix}$

where $A=0.1$, $a=0.1\sin\left({\pi t}/{5}\right)$, and $b=1-2a$. Below, we
provide a visualization of our quantity compared with FTLE. We then analyze
the trajectory in the part that demonstrates differences between the two
metrics in the visualization. We also provide the results of an actual runtime
experiment, where we vary the mesh size and trajectory length. Finally, we
compare the effects of using different relaxation schemes.

The results of the computation are visualized in Figures 4(a-c). Comparing
these results to the FTLE field plotted in Figure 4(d), we can see that our
proposed quantity forms a sharp and clean ridge that agrees with the FTLE
ridge. However, there is a slight difference in the tail of the structure
between the results obtained from the three metrics. To provide a better
visualization, we plotted the LTVE values computed under different metrics and
overlapped them in Figure 5(a). Specifically, we selected a point located at
the tail and plotted the trajectories at that point and in its neighborhood
for comparison, as shown in Figure 5(b). Upon closer inspection, we can
observe that the trajectory at the particular point does not exhibit a strong
expelling or attracting feature. However, there are slight variations in the
trajectories, as seen in the region boxed in red and visualized in Figures
5(c-d). We notice that although the result evaluated from the Fréchet distance
appears to be more consistent with the tail part of the FTLE, Figures 5(c-d)
show that the perturbation of the trajectories is not as significant as
expected for being a LCS. This indicates that even along the FTLE ridge, the
effect of perturbing the trajectories can still vary. The structure captured
at a certain point may be much weaker, and our proposed quantity can capture
such features by using a different metric.

To validate the relationship between the LTVE and FTLE, we computed the actual
differences at different times and compared them with the estimations derived
from equation (3). The resulting differences are illustrated in Figure 6. In
the figure, each blue point represents the theoretical bound estimated from
equation (3), while each red point represents the actual difference evaluated
when fixing the time span to $[0,T]$. We observe that the red points,
representing the actual differences, consistently fall below the blue points,
which correspond to the theoretical estimations. This alignment between the
experimental results and our derived equation (3) confirms the validity of our
theoretical derivation.

(a) (b) (c) (d)

Figure 7: (Section 4.5) The computed LTVE using the normalized 2-norm with (a)
first order relaxation scheme (b) second order relaxation scheme and (c) third
order relaxation scheme, (d) is the FTLE field for comparison.

###  4.5 Relaxation Schemes

Here, we present a comparative study on the effect of different relaxation
schemes on the visualization results of the double gyre flow. From Figure 7,
we can observe that there is not a significant difference between (a) and (b).
However, in (c), our quantity highlights a thicker line of structure.
Furthermore, the structure that we focused on in Figure 5(a) is now consistent
with the result obtained using the Fréchet distance. This result aligns with
our claim in Section 4.4 that the structure appearing at the tail is weaker
and more sensitive to the initial condition. Therefore, by adopting a higher-
order relaxation scheme, it is more likely to capture the structure that only
becomes significant along certain directions.

(a)  Metric mesh size $\delta$ Normalized $L^{2}$-norm Fréchet distance
Hausdorff metric ${1}/{49}$ $0.209$ $21$ $47$ ${1}/{99}$ $0.835$ $85$ $193$
${1}/{149}$ $1.878$ $192$ $436$ ${1}/{199}$ $3.334$ $342$ $776$ ${1}/{249}$
$5.225$ $535$ $1231$
(b)  Metric trajectory length $k$ Normalized $L^{2}$-norm Fréchet distance
Hausdorff metric $5$ $0.175$ $0.355$ $0.461$ $10$ $0.196$ $0.912$ $1.392$ $15$
$0.226$ $1.831$ $2.857$ $20$ $0.250$ $3.051$ $4.893$ $25$ $0.277$ $4.565$
$7.415$ $30$ $0.304$ $6.454$ $10.506$ $35$ $0.330$ $8.578$ $14.166$ $40$
$0.370$ $11.459$ $18.520$ $45$ $0.383$ $14.208$ $23.864$ $50$ $0.427$ $17.457$
$28.483$ $55$ $0.463$ $21.980$ $34.331$ $60$ $0.463$ $24.122$ $41.441$

Table 1: (Section 4.6) (a) The CPU run time in seconds for different mesh
sizes over different metrics. (b) The CPU run time in seconds for different
lengths of trajectory over different metrics.

(a) (b)

Figure 8: (Section 4.6) (a) The log time versus $\log\delta^{-2}$, and (b) The
log time versus the log of the trajectory length.

###  4.6 Runtime Performance

In this section, we provide the actual CPU run time results while varying the
trajectory length and mesh size. All the results were obtained using a laptop
computer with a 4-core 11th Gen Intel i5 CPU and 16GB of RAM. It is important
to note that our code is not fully optimized, and the main purpose of
reporting these runtimes is to showcase the complexity under different
metrics. To conduct our measurements, we utilized the double gyre flow from
Section 4.4 as the underlying field. This choice was motivated by its
complexity, which allows us to present a comprehensive example for visualizing
LCS.

##### Mesh Size.

First, we conducted performance measurements using five selected mesh sizes,
denoted as $\delta$, while keeping other parameters the same as in our
previous experiments. The recorded times are presented in Table 1(a), and we
plotted the logarithm of time against $\log{\delta^{-2}}$ in Figure 8(a).
Notably, the figure exhibits a clear linear relationship with a slope of one,
which aligns with our complexity analysis. Furthermore, based on our
complexity analysis, the $y$-intercept of the plot should depend on the domain
size and the complexity of the metric used to compare trajectories. In this
case, since the domain size is fixed, the slope is primarily determined by the
metric. This observation is supported by the data, as the plot of the
normalized $L^{2}$ norm displays a significantly smaller $y$-intercept
compared to the other two metrics, indicating its linear complexity. Regarding
the Fréchet distance and Hausdorff metric, the computation of the Fréchet
distance involves dynamic programming, resulting in a relatively lower actual
runtime compared to the Hausdorff metric, which requires comparing all pairs
of trajectory points to compute the maximum distance. This explains why the
Fréchet distance has a lower runtime compared to the Hausdorff metric.

##### Trajectory Length.

We then measured the performance for five different trajectory lengths,
keeping the other parameters constant. The results are presented in Table
1(b). Once again, the normalized $L^{2}$-norm outperformed the other two
metrics due to its linear complexity. The Fréchet distance metric performed
slightly better than the Hausdorff metric, consistent with our previous
analysis. Our complexity analysis suggests that the total runtime should be
linear with the complexity of the metric, which is confirmed by the plot in
Figure 8(b). The slope of the plot is upper bounded by the actual complexity
of the metric. The slight difference in runtime between the metrics may be due
to hidden compiler optimizations that are hard to control. Nevertheless, the
complexity analysis provides a valid upper bound for the runtime.

###  4.7 Field with Artificial FTLE Ridge

(a) (b) (c) (d)

Figure 9: (Section 4.7) The computed LTVE using (a) the normalized 2-norm, (b)
the Fréchet distance, and (c) the Hausdorff metric. (d) The corresponding FTLE
field.

In [20], an Eulerian approach was proposed for computing the FTLE field, which
highlighted the issue of artificial ridges arising from the boundary condition
upon solving the trajectories. To illustrate this issue, a specific example
was presented, where the velocity field is given by:

$\frac{d\mathbf{x}}{dt}=\begin{pmatrix}x-y^{2}\ -y+x^{2}\end{pmatrix}$

For this example, the computational domain chosen was $[-6,6]\times[-6,6]$
with a spatial resolution of $\delta={12}/{599}$. The time range was set to
$[0,5]$, with 100 sampling points for each trajectory.

As depicted in Figure 9(a-c), our proposed method yields results that closely
resemble those obtained using the FTLE field. However, it is worth noting that
our method generally provides a cleaner visualization due to the sharp values
at the structural points. It is important to acknowledge that this example
contains an artificial LCS that does not exist in its natural definition. To
address this issue, various techniques can be employed, such as the velocity
extension algorithm proposed in [34] or exploring an Eulerian approach for our
quantity as demonstrated in [20], which may help mitigate the presence of
artificial ridges.

###  4.8 Ocean Surface Current Analyses Real-time (OSCAR)

(a) (b)
(c) (d)

Figure 10: (Section 4.8) The computed LTVE using (a) the normalized 2-norm,
(b) the Fréchet distance, and (c) the Hausdorff metric. (d) The corresponding
FTLE field [47, 44].

To assess the efficacy and reliability of our proposed approach, we conducted
experiments on a real-life dataset sourced from the Ocean Surface Current
Analyses Real-time (OSCAR)111http://www.esr.org/oscar_index.html. The OSCAR
dataset captures ocean surface flow in a wide region spanning from
$-80^{\circ}$ to $80^{\circ}$ latitude and $0^{\circ}$ to $360^{\circ}$
longitude. We acquired the OSCAR data from the JPL Physical Oceanography DAAC
developed by ESR. In this particular numerical example, our focus centers on
the vicinity of the Line Islands, encompassing the latitude range of
$-17^{\circ}$ to $8^{\circ}$ and the longitude range of $180^{\circ}$ to
$230^{\circ}$. The temporal resolution of the OSCAR data is approximately 5
days, and the spatial resolution is $1/3^{\circ}$ in each direction.

To enhance the visual representation of our solution, we performed velocity
data interpolation, yielding a finer resolution of $0.25$ days in the temporal
dimension and $1/12^{\circ}$ in each spatial dimension. Subsequently, we
examined the ocean surface current during the first 50 days of the year 2014
to derive the forward FTLE field $\sigma_{0}^{50}(\mathbf{x})$. The
computational results obtained from our proposed algorithm, utilizing various
metrics, are presented in Figures 10(a-c). Notably, these results agree with
the FTLE field computed in [47, 44], as depicted in Figure 10(d), in general,
but show differences in various regions. This provides insight into the flow
structures that are worth further investigation.

##  5 Conclusions

In this study, we propose a novel quantity based on the trajectory metric
space, which bears resemblance to the FTLE field. We evaluate the efficacy of
this new quantity with several examples, including one that reveals
significant differences between our quantity and the FTLE field. This example
serves as a motivation for the use of trajectory analysis in LCS
visualization. We further test our quantity with other examples that
demonstrate its ability to handle complex flows that arise in natural and
certain flow conditions. In all cases, we observe that our quantity captures
the general structure that the FTLE field captures, albeit with slight
differences contingent upon the specific metric employed. We find that the
differences between our quantity and the FTLE field are minor in the field we
examined, but there may be fields where significant differences exist.
However, the direct application of the Fréchet distance or the Hausdorff
metric is prohibitively expensive, and the time cost does not justify the
marginal improvement in results.

In conclusion, our proposed method represents a generalization of the
traditional FTLE approach with trajectories for visualizing Lagrangian
coherent structures (LCS). We recommend utilizing the normalized $L^{2}$-norm
as the primary metric, as it generally provides efficient and satisfactory
results without significant differences compared to other metrics. Unless a
specific application case presents a compelling reason to justify the use of a
different metric, the normalized $L^{2}$-norm is sufficient for capturing the
desired flow structures. To the best of our knowledge, no such application
case has been identified thus far. However, we plan to further investigate
this aspect in future studies to explore potential scenarios where alternative
metrics may offer distinct advantages. By doing so, we aim to enhance our
understanding of the implications and potential benefits of different metrics
for visualizing LCS.

Although we did not present any experiments of computing our quantity on
higher dimensional domains in this paper, we believe that the framework for
our quantity is general and not specifically dependent on the dimension,
making it applicable to fields with higher dimensions.

In future research, we plan to extend our framework to dynamic cases and
explore scenarios where a different metric could lead to a significantly
different and useful result. Such results would be intriguing because they
would also imply that the classical FTLE approach fails in these cases as it
is a special case of applying the normalized $L^{2}$ norm. Additionally, we
plan to investigate algorithms that can accelerate metric calculations since
we found that metric complexity significantly affects practical runtime. To
achieve dynamic cases, we need to identify a way of updating trajectories
efficiently, enabling online updates of the quantity in a practical way.

## Acknowledgments

Leung is partially supported by the Hong Kong RGC (grant no. 16302223).

## References

  * [1] P.K. Agarwal, R.B. Avraham, H. Kaplan, and M. Sharir.  Computing the Discrete Fréchet Distance in Subquadratic Time.  CoRR, abs/1204.5333, 2012.
  * [2] D. Blazevski and G. Haller.  Hyperbolic And Elliptic Transport Barriers In Three-Dimensional Unsteady Flows.  Physica D: Nonlinear Phenomena, 273-274:46–62, 2014.
  * [3] K. Bringmann.  Why Walking The Dog Takes Time: Frechet Distance Has No Strongly Subquadratic Algorithms Unless SETH Fails.  CoRR, abs/1404.1448, 2014.
  * [4] B.M. Cardwell and K. Mohseni.  Vortex shedding over two-dimensional airfoil: Where do the particles come from?  AIAA J., 46:545–547, 2008.
  * [5] T.M. Chan and Z. Rahmati.  An Improved Approximation Algorithm for the Discrete Fréchet Distance.  Information Processing Letters, 138:72–74, 2018.
  * [6] W.M. Chau and S. Leung.  Within-Cluster Variability Exponent for Identifying Coherent Structures in Dynamical Systems.  Commun. Comput. Phys., 33(3):824–848, 2023.
  * [7] T. Eiter and H. Mannila.  Computing Discrete Fréchet Distance.  Technical report, Technical University of Vienna, 1994.
  * [8] N. Ferreira, J.T. Klosowski, C.E. Scheidegger, and C.T. Silva.  Vector Field k-Means: Clustering Trajectories by Fitting Multiple Vector Fields.  Eurographics Conference on Visualization (EuroVis), 32(3):201–210, 2013.
  * [9] S. Gaffney, A.W. Robertson, P. Smyth, S.J. Camargo, and M. Ghil.  Probabilistic clustering of extratropical cyclones using regression mixture models.  Clim. Dyn., 29:423–440, 2007.
  * [10] S. Gaffney and P. Smyth.  Trajectory Clustering with Mixtures of Regression Models.  Proc. 5th ACM SIGKDD Int’l Conf. on Knowledge Discovery and Data Mining, pages 63–72, 1999.
  * [11] M.A. Green, C.W. Rowley, and A.J. Smiths.  Using hyperbolic Lagrangian coherent structures to investigate vortices in biospired fluid flows.  Chaos, 20:017510, 2010.
  * [12] G. Haller.  Distinguished material surfaces and coherent structures in three-dimensional fluid flows.  Physica D, 149:248–277, 2001.
  * [13] G. Haller.  Lagrangian structures and the rate of strain in a partition of two-dimensional turbulence.  Phys. Fluids A, 13:3368–3385, 2001.
  * [14] G. Haller.  A variational theory of hyperbolic Lagrangian Coherent Structure.  Physica D, 240:574–598, 2011.
  * [15] G. Haller.  Lagrangian coherent structures.  Annu. Rev. Fluid Mech., 47:137–162, 2015.
  * [16] G. Haller and G. Yuan.  Lagrangian coherent structures and mixing in two-dimensional turbulence.  Physica D, 147:352–370, 2000.
  * [17] J. Lee, J. Han, and K. Whang.  Trajectory clustering: a partition-and-group framework.  ACM SIGMOD International Conference on Management of Data, pages 593–604, 2007.
  * [18] F. Lekien and N. Leonard.  Dynamically consistent Lagrangian coherent structures.  Experimental Chaos: 8-th Experimental Chaos Conference, pages 132–139, 2004.
  * [19] F. Lekien, S.C. Shadden, and J.E. Marsden.  Lagrangian coherent structures in $n$-dimensional systems.  Journal of Mathematical Physics, 48:065404, 2007.
  * [20] S. Leung.  An Eulerian approach for computing the finite time Lyapunov exponent.  J. Comput. Phys., 230:3500–3524, 2011.
  * [21] S. Leung.  The backward phase flow method for the finite time Lyapunov exponent.  Chaos, 23(043132), 2013.
  * [22] S. Leung, G. You, T. Wong, and Y.K. Ng.  Recent developments in Eulerian approaches for visualizing continuous dynamical systems.  Proceedings of the Seventh International Congress of Chinese Mathematicians, (2):579–622, 2019.
  * [23] D. Lipinski and K. Mohseni.  Flow structures and fluid transport for the hydromedusae Sarsia tubulosa and Aequorea victoria.  J. Exp. Biology, 212:2436–2447, 2009.
  * [24] S. Lukens, X. Yang, and L. Fauci.  Using Lagrangian coherent structures to analyze fluid mixing by cillia.  Chaos, 20:017511, 2010.
  * [25] S. Mowlavi, M. Serra, E. Maiorino, and L. Mahadevan.  Detecting Lagrangian Coherent Structures from Sparse and Noisy Trajectory Data.  Journal of Fluid Mechanics, 948, sep 2022.
  * [26] A. Nath and E. Taylor.  k-Median Clustering under Discrete Fréchet and Hausdorff Distances.  CoRR, abs/2004.00722, 2020.
  * [27] Y.K. Ng and S. Leung.  Estimating the Finite Time Lyapunov Exponent from Sparse Lagrangian Trajectories.  Commun. Comput. Phys., 26(4):1143–1177, 2019.
  * [28] G.A. Rosi, A.M. Walker, and D.E. Rival.  Lagrangian Coherent Structure Identification Using A Voronoi Tessellation-Based Networking Algorithm.  Experiments in Fluids, 56:1–14, 2015.
  * [29] T. Sapsis and G. Haller.  Inertial particle dynamics in a hurricane.  Journal of the Atmospheric Sciences, 66:2481–2492, 2009.
  * [30] K.L. Schlueter-Kuck and J.O. Dabiri.  Coherent Structure Colouring: Identification Of Coherent Structures From Sparse Data Using Graph Theory.  Journal of Fluid Mechanics, 811:468 – 486, 2016.
  * [31] S.C. Shadden, F. Lekien, and J.E. Marsden.  Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows.  Physica D, 212:271–304, 2005.
  * [32] H. Su, S. Liu, B. Zheng, X. Zhou, and K. Zheng.  A Survey of Trajectory Distance Measures and Performance Evaluation.  The VLDB Journal, 29(1):3–32, oct 2019.
  * [33] W. Tang, P.W. Chan, and G. Haller.  Accurate extraction of Lagrangian coherent structures over finite domains with application to flight data analysis over Hong Kong international airport.  Chaos, 20:017502, 2010.
  * [34] W. Tang, P.W. Chan, and G. Haller.  Accurate Extraction Of Lagrangian Coherent Structures Over Finite Domains With Application To Flight Data Analysis Over Hong Kong International Airport.  Chaos, 20(1):017502, 2010.
  * [35] W. Tang and T. Peacock.  Lagrangian coherent structures and internal wave attractors.  Chaos, 20:017508, 2010.
  * [36] K. Toohey and M. Duckham.  Trajectory Similarity Measures.  SIGSPATIAL Special, 7(1):43–50, may 2015.
  * [37] J. Tyler and A. Wittig.  An Improved Numerical Method For Hyperbolic Lagrangian Coherent Structures Using Differential Algebra.  Journal of Computational Science, 65:101883, 2022.
  * [38] M. Vlachos, D. Gunopulos, and G. Kollios.  Discovering Similar Multidimensional Trajectories.  Proc. 18th Int’l Conf. on Data Engineering, pages 673–684, 2002.
  * [39] J. Wei, H. Yu, J. Chen, and K. Ma.  Parallel clustering for visualizing large scientific line data.  IEEE Symposium on Large Data Analysis and Visualization, 1:47–55, 2011.
  * [40] D. Yao, C. Zhang, Z. Zhu, J. Huang, and J. Bi.  Discovering Similar Multidimensional Trajectories.  2017 International Joint Conference on Neural Networks (IJCNN), pages 3880–3887, 2017.
  * [41] G. You and S. Leung.  An Eulerian method for computing the coherent ergodic partition of continuous dynamical systems.  J. Comp. Phys., 264:112–132, 2014.
  * [42] G. You and S. Leung.  Eulerian based interpolation schemes for flow map construction and line integral computation with applications to coherent structures extraction.  J. Sci. Comput., 74(1):70–96, 2018.
  * [43] G. You and S. Leung.  An improved eulerian approach for the finite time lyapunov exponent.  J. Sci. Comput., 76(3):1407–1435, 2018.
  * [44] G. You and S. Leung.  Fast Construction of Forward Flow Maps using Eulerian Based Interpolation Schemes.  J. Sci. Comput., 82(32), 2020.
  * [45] G. You and S. Leung.  Computing the finite time lyapunov exponent for flows with uncertainties.  J. Comput. Phys., 425(109905), 2021.
  * [46] G. You and S. Leung.  Eulerian algorithms for computing some lagrangian flow network quantities.  J. Comput. Phys., 445(11020), 2021.
  * [47] G. You, T. Wong, and S. Leung.  Eulerian methods for visualizating continuous dynamical systems using Lyapunov exponents.  SIAM J. Sci. Comp., 39(2):A415–A437, 2017.