Robust FEM-based extraction of finite-time coherent sets using scattered, sparse, and incomplete trajectories

TL;DR

This paper introduces FEM-based numerical methods to identify finite-time coherent sets in flows using limited, scattered observational data, aiding the analysis of transport and mixing in complex systems.

Contribution

The paper develops three FEM-based algorithms for approximating the dynamic Laplace operator and introduces a new dynamic isoperimetric problem with Dirichlet conditions for better coherent set detection.

Findings

Efficient extraction of finite-time coherent sets from sparse data.

Reliable identification of flow regions with minimal mixing.

Applicable to real-world geophysical flow data.

Abstract

Transport and mixing properties of aperiodic flows are crucial to a dynamical analysis of the flow, and often have to be carried out with limited information. Finite-time coherent sets are regions of the flow that minimally mix with the remainder of the flow domain over the finite period of time considered. In the purely advective setting this is equivalent to identifying sets whose boundary interfaces remain small throughout their finite-time evolution. Finite-time coherent sets thus provide a skeleton of distinct regions around which more turbulent flow occurs. They manifest in geophysical systems in the forms of e.g.\ ocean eddies, ocean gyres, and atmospheric vortices. In real-world settings, often observational data is scattered and sparse, which makes the difficult problem of coherent set identification and tracking even more challenging. We develop three FEM-based numerical methods to efficiently approximate the dynamic Laplace operator, and introduce a new dynamic isoperimetric problem using Dirichlet boundary conditions. Using these FEM-based methods we rapidly and reliably extract finite-time coherent sets from models or scattered, possibly sparse, and possibly incomplete observed data.

Figures

Equations

$$\mathbf{h}(\Gamma ):=\frac{\frac{1}{|\mathcal{T}|}{\sum }_{t\in \mathcal{T}}{\ell }_{d-1}({\Phi }^t\Gamma )}{\min \{\ell ({\mathcal{M}}_1),\ell ({\mathcal{M}}_2)\}},$$

$$\mathbf{h}:=\min\{\mathbf{h}(\Gamma):\Gamma\mbox{ is a $C^{\infty}$ co-dimension 1 manifold disconnecting $\mathcal{M}$}\}.$$

$$\mathbf{s}:=\underset{f\in C^{\infty }(\mathcal{M},\mathbb{R})}{\inf }\frac{\frac{1}{|\mathcal{T}|}{\sum }_{t\in \mathcal{T}}\Vert \nabla ({\Phi }_{\ast }^{t}f){\Vert }_1}{{\inf }_{\alpha \in \mathbb{R}}\Vert f-\alpha {\Vert }_1}$$

$$\left(\frac{1}{|\mathcal{T}|}\sum _{t\in \mathcal{T}}({\Phi }^t{)}^{\ast }\Delta {\Phi }_{\ast }^t\right)v=\lambda v\text{mbox}onint(\mathcal{M}),$$

$${\int }_{\mathcal{M}}\left(\frac{1}{|\mathcal{T}|}\sum _{t\in \mathcal{T}}({\Phi }^t{)}^{\ast }\Delta {\Phi }_{\ast }^t\right)v\cdot \psi d\ell$$

$$-\frac{1}{|\mathcal{T}|}\sum _{t\in \mathcal{T}}{\int }_{{\Phi }^t(\mathcal{M})}\nabla ({\Phi }_{\ast }^{t}v)\bullet \nabla ({\Phi }_{\ast }^t\psi )d\ell =\lambda {\int }_{\mathcal{M}}v\psi d\ell ,\text{for all }\psi \in H^1(\mathcal{M}).$$

$${\mathbf{h}}_d(A):=\frac{\frac{1}{|\mathcal{T}|}{\sum }_{t\in \mathcal{T}}{\ell }_{d-1}({\Phi }^t(\partial A))}{\ell (A)},$$

$${\mathbf{h}}_d:=\min \{\mathbf{h}(A):A\text{ open, submanifold of mboxint(M), compact closure, C∞ boundary}\}.$$

$${\mathbf{s}}_d:=\underset{f\in C_c^{\infty }(\mathcal{M},\mathbb{R})}{\inf }\frac{\frac{1}{|\mathcal{T}|}{\sum }_{t\in \mathcal{T}}\Vert \nabla ({\Phi }_{\ast }^{t}f){\Vert }_1}{\Vert f{\Vert }_1},$$

$$\left(\frac{1}{|\mathcal{T}|}\sum _{t\in \mathcal{T}}({\Phi }^t{)}^{\ast }\Delta {\Phi }_{\ast }^t\right)v$$

$$v$$

$$-\frac{1}{|\mathcal{T}|}\sum _{t\in \mathcal{T}}{\int }_{{\Phi }^t(\mathcal{M})}\nabla ({\Phi }_{\ast }^{t}v)\bullet \nabla ({\Phi }_{\ast }^t\psi )d\ell =\lambda {\int }_{\mathcal{M}}v\psi d\ell ,\text{for all }\psi \in H_0^1(\mathcal{M}).$$

$$-\left(\frac{1}{|\mathcal{T}|}\sum _{t\in \mathcal{T}}{\hat{D}}^t\right)u=\lambda Mu,$$

$${\hat{D}}_{ij}^t={\int }_{{\Phi }^t(\mathcal{M})}\nabla \left({\Phi }_{\ast }^t{\phi }_i\right)\bullet \nabla \left({\Phi }_{\ast }^t{\phi }_j\right)d\ell \text{and}{M}_{ij}={\int }_{\mathcal{M}}{\phi }_i{\phi }_{j}d\ell$$

$${\hat{D}}_{ij}^t={\int }_{{\Phi }^t(\mathcal{M})}\nabla \left({\Phi }_{\ast }^t{\phi }_i\right)\bullet \nabla \left({\Phi }_{\ast }^t{\phi }_j\right)d{\mu }^t\text{and}{M}_{ij}={\int }_{\mathcal{M}}{\phi }_i{\phi }_{j}d{\mu }^0,$$

$${\int }_{{\Phi }^t(\mathcal{M})}\nabla \left({\Phi }_{\ast }^t{\phi }_i\right)\bullet \nabla \left({\Phi }_{\ast }^t{\phi }_j\right)d{\mu }^t$$

$$={\int }_{\mathcal{M}}(D{\Phi }^t{)}^{-\top }\nabla {\phi }_i\bullet (D{\Phi }^t{)}^{-\top }\nabla {\phi }_{j}d{\mu }^0\text{mbox}(chainrule)$$

$$={\int }_{\mathcal{M}}\nabla {\phi }_i\bullet C_t^{-1}\nabla {\phi }_{j}d{\mu }^0$$

$${\Phi }_{\ast }^t{\phi }_i\approx \sum _{k=1}^n{\alpha }_k^i{\phi }_k^t,i=1,\dots ,n,$$

$${\hat{D}}_{ij}^t={\int }_{{\Phi }^t(\mathcal{M})}\nabla \left({\Phi }_{\ast }^t{\phi }_i\right)\bullet \nabla \left({\Phi }_{\ast }^t{\phi }_j\right)d\ell$$

$$D_{k\ell }^t:={\int }_{{\Phi }^t(\mathcal{M})}\nabla {\phi }_k^t\bullet \nabla {\phi }_{\ell }^{t}d{\mu }^t$$

$${\tilde{h}}^t=\sum _{i=1}^n\frac{a_i}{{\ell }_i^t}{\phi }_i^t,$$

$$D_{k\ell }^t={\int }_{{\Phi }^t(\mathcal{M})}(\nabla {\phi }_k^t\bullet \nabla {\phi }_{\ell }^t)d{\mu }^t\approx \sum _{e^t}{\int }_{e^t}(\nabla {\phi }_k^t\bullet \nabla {\phi }_{\ell }^t){\tilde{h}}^{t}d\ell .$$

$${\int }_{e^t}(\nabla {\phi }_k^t\bullet \nabla {\phi }_{\ell }^t){\tilde{h}}^{t}d\ell =\sum _{x_i^t\in e^t}\frac{a_i}{{\ell }_i^t}{\int }_{e^t}(\nabla {\phi }_k^t\bullet \nabla {\phi }_{\ell }^t){\phi }_i^{t}d\ell .$$

$${\int }_{e^t}(\nabla {\phi }_k^t\bullet \nabla {\phi }_{\ell }^t){\tilde{h}}^{t}d\ell$$

$$=(\nabla {\phi }_k^t\bullet \nabla {\phi }_{\ell }^t)\text{area}(e^t)\frac{1}{3}\sum _{x_i^t\in e^t}\frac{a_i}{{\ell }_i^t}.$$

$${\langle\Phi^{t}_{*}\varphi_{i},\delta_{x^{t}_{m}}\rangle}=\sum_{k=1}^{n}\alpha^{i}_{k}{\langle\varphi^{t}_{k},\delta_{x^{t}_{m}}\rangle},\qquad\mbox{for all $1\leq i,m\leq n$},$$

$${\alpha }_m^i={\Phi }_{\ast }^t{\phi }_i(x_m^t)={\phi }_i({\Phi }^{-t}(x_m^t)).$$

$${\Phi^{t}_{*}\varphi_{i}(x^{t}_{m})h_{\mu^{t}}(x^{t}_{m})}=\sum_{k=1}^{n}\alpha^{i}_{k}{\varphi^{t}_{k}(x^{t}_{m})h_{\mu^{t}}(x^{t}_{m})},\qquad\mbox{for all $1\leq i,m\leq n$}.$$

$${\alpha }_m^i={\phi }_i({\Phi }^{-t}{x}_m^t)={\phi }_i(x_m)={\delta }_{i,m},$$

$$D^{t}_{ij}:=\left\{\begin{array}[]{ll}\int_{\Phi^{t}(\mathcal{M})}\nabla\varphi_{i}^{t}\cdot\nabla\varphi_{j}^{t}\ d\ell,&\hbox{$i,j\in I_{t}$},\\ 0,&\hbox{otherwise.}\end{array}\right.$$

$${\int }_{\mathcal{M}}{\phi }_i^0\cdot {\phi }_j^{0}d\ell ={\int }_{{\Phi }^t(\mathcal{M})}{\phi }_i^0\circ {\Phi }^{-t}\cdot {\phi }_j^0\circ {\Phi }^{-t}d\ell .$$

Full text

1 Introduction

The quantification of transport and mixing processes in nonlinear dynamical systems is key to understanding their global dynamics. In the autonomous setting, where the governing dynamical laws do not change over time, stable and unstable manifolds of low-period orbits have been used to describe the global dynamics. For example, co-dimension 1 invariant manifolds represent an impenetrable barrier to trajectories. In the case of periodically-forced flows, the intersections of stable and unstable manifolds form mobile parcels dubbed “lobes”, which give rise to lobe dynamics [29, 33] as a mechanism for transport across periodically varying invariant manifolds.

In the autonomous setting, the concept of almost-invariant sets, introduced in [9] and further developed in [16, 13] uses the transfer (or Perron-Frobenius) operator to find regions in phase space which are either invariant, or in the absence of any open invariant sets, as close to invariant as possible. Such sets decompose the phase space into regions which are as dynamically disconnected as possible, and between which transport is minimized. Related approaches to finding “ergodic partitions” in phase space have also been developed using long orbits [31] or using the Koopman operator [5], the dual of the transfer operator.

In the case of nonautonomous, aperiodic dynamics, the notion of almost-invariant sets was generalized to coherent sets [19] and shortly thereafter to finite-time coherent sets [22, 14]. While almost-invariant sets remain fixed in phase space, when the underlying governing dynamics is time-varying, it is essential to allow sets that are parameterised by time. Nevertheless, the underlying idea, minimising mixing with the surrounding phase space, is common to both almost-invariant and finite-time coherent sets. As in the case of almost-invariant sets, one numerically constructs a transfer operator for the duration of the dynamics to be considered. This transfer operator will push forward densities over the chosen time interval. After some straightforward normalising procedures, one computes the subdominant singular spectrum and singular vectors of the transfer operator, which indicate the presence and location of finite-time coherent sets. This construction relies on a small amount of diffusion being present in the dynamics; this diffusion can be part of the model, or if the underlying dynamics is purely advective it can be easily numerically induced.

An equivalent definition of finite-time coherent sets was introduced in [15], specifically aimed at the case of purely advective dynamics. The underlying goal is to determine co-dimension 1 manifolds whose size remains small relative to enclosed333The boundary of the domain may take part in the enclosing. volume, under evolution by the dynamics. This idea extends the classical (static) isoperimetric problem [4] of determining the co-dimension 1 manifold with least size given a fixed enclosed volume. The idea to define finite-time coherent sets as those sets with boundaries that have dynamically minimal size relative to enclosed volume is consistent with the minimal $L^{2}$-mixing approach of [14] because if small-amplitude diffusion is added, $L^{2}$ mixing can only occur along the boundary interfaces between coherent regions. This fact has been formally established in [15] in the volume-preserving setting and [18] in the non-volume-preserving case. The persistently minimal size of the boundary of finite-time coherent sets in this sense is also broadly consistent with the approach of [26], although [26] is limited to two-dimensions and closed curves. Related work [27] uses a time-dependent metric (with corresponding time-dependent diffusion) induced by the advective flow to define coherent sets as almost-invariant sets in the material manifold. In approximating the time dependent metric by a single one they end up with a similar operator as the dynamic Laplacian in [15, 18].

Numerical implementation of the method of [15] requires approximation of eigenfunctions of a dynamic Laplace operator, which consists of the Laplace operator composed with transfer operators, or equivalently, combinations of Laplace-Beltrami operators with respect to Riemannian metrics pulled back under the nonlinear dynamics. The papers [15, 18] used low-order methods such as finite-differences to approximate the Laplace operator and Ulam’s method [37] to approximate the transfer operators. In [17], the authors investigated the use of radial basis functions (RBFs) as an approximating basis for $L^{2}$, seeking to exploit the smoothness of the underlying dynamics to achieve accurate results with relatively sparse trajectory information. This approach successfully achieved very accurate results with relatively sparse data, but due to known issues with RBF interpolation [39, 12] care must be taken when selecting the RBF centre spacing and RBF radii. In [3], diffusion maps on possibly scattered trajectory data were used for the approximation of the dynamic Laplacian, and in [24], an object is constructed which is related to a discretised version of the dynamic Laplace operator, again from possibly scattered trajectory data.

The dynamic Laplacian is a self-adjoint, elliptic differential operator and arguably the finite element method (FEM) is the standard numerical method to discretize these. In our specific dynamic context, there are a number of major advantages to this method:

• •

Sparsity: We obtain a sparse discrete problem. Thus, the computational cost grows only linearly with the number of mesh elements and enables one to use high resolution meshes (in contrast to the RBF approach [17]).

• •

Structure preservation: The dynamic Laplacian is self-adjoint and thus its weak form is symmetric. This symmetry is inherited in the discrete eigenproblem, ensuring a purely real spectrum (in contrast to, e.g., the RBF method).

• •

*Arbitrary domains: * FEMs can elegantly handle arbitrarily shaped domains, as often occur in applications, unlike spectral methods, which require tensor product domains.

• •

Natural boundary conditions: Homogeneous Neumann boundary conditions are the natural ones in our context and completely disappear in the weak formulation. Homogeneous Dirichlet conditions can also easily be incorporated, as we introduce and describe here.

We would further like to highlight advantages over other existing methods for extracting (finite time) coherent sets:

There is no free parameter to choose (like a cut-off radius, a neighborhood size etc.) experimentally or heuristically when computing the eigenvectors. When extracting coherent sets from these, some heuristics may be needed to choose how many/which eigenvectors yield relevant information. In many cases, this choice is supported by an eigengap heuristic. 2. 2.

Even in the case of scattered and sparse trajectory data, we always obtain coherence information for every point in phase space. 3. 3.

Our method applies to vector field data and discretely sampled trajectory data alike. The data can be given on a scattered grid, can be sparse and incomplete Of course, the predictive power will decrease with decreasing data density, however, we observe robust and reliable results for highly nonlinear flows with relatively little data. 4. 4.

Our FEM approach does not require derivative information on the flow map, but can make use of it when it is available.

An outline of the paper is as follows: In Section 2 and 2.1 we provide a brief overview of the dynamic isoperimetric problem, the associated dynamic Laplacian eigenproblem, and state the weak form of this eigenproblem. In Section 2.2 we introduce and motivate the new Dirichlet form of the dynamic isoperimetric problem, and state fundamental theoretical results (dynamic Federer-Fleming theorem and dynamic Cheeger’s inequality) corresponding to the existing Neumann theory. In Section 3 we first introduce the discrete form of the weak eigenproblem, which involves gradients of push-forwards of elements under the transfer operator. In Section 3.1 we convert the gradient of the push-forward of an element into the push-forward of the gradient of an element which results in a stiffness matrix containing the Cauchy-Green tensor; standard FEM can then be applied. In Section 3.2 we instead estimate the push-forward of elements under the transfer operator. We first describe the case of volume-preserving dynamics, and then in Section 3.2.1 detail the situations where (i) the dynamics is not volume-preserving or (ii) one wishes to extract coherent features with respect to non-uniform densities; cases (i) and (ii) are handled identically. To implement the transfer operator approach we use two forms of collocation described in Sections 3.2.2–3.2.3. In this case, one uses the standard stiffness matrix for the original elements, multiplied by a matrix representing the transfer operator. In the second collocation form (Section 3.2.3), if the input data comes from trajectories, estimation of the transfer operator is completely avoided (cf. Algorithm 1 and the code in the appendix).

The transfer operator approach is extremely efficient and robust, particularly in the situation where the trajectory information is sparse. Section 3.2.4 describes the minor modifications required to handle trajectories where many or most data points are missing. This is particularly useful if trajectories arise from physical observations, which may not be recorded at certain times due to mechanical breakdown, adverse weather, or other obstructions. Section 3.3 recalls how to extract multiple coherent sets from multiple eigenfunctions, and Section 3.4 interprets the adaptive transfer operator collocation method (section 3.2.3) in terms of graph Laplacians.

Section 4 contains numerical experiments that illustrate the robustness and speed of the proposed methods on a variety of models and data availability. The dynamical systems tested include two- and three-dimensional idealized models, including a model of ocean flow, and experiments with scattered, strongly irregularly distributed, and missing data. In the appendices, we provide brief theoretical background for the non-volume preserving case and provide sample code for one the examples. The corresponding MATLAB package FEMDL can be downloaded from https://github.com/gaioguy/FEMDL.

2 Finite-time coherent sets and dynamic Laplacians

Let $\mathcal{M}$ be a compact subset of $\mathbb{R}^{d}$ with empty or piecewise smooth boundary, or more generally a compact, smooth, Riemannian manifold. For $t\in[0,t_{F}]$, we consider diffeomorphisms $\Phi^{t}\colon\mathcal{M}\subset\mathbb{R}^{d}\to\Phi^{t}(\mathcal{M})\subset\mathbb{R}^{d}$, representing e.g. flow maps from time 0 to time $t$. Let $\Phi^{t}_{*}:C^{\infty}(\mathcal{M})\to C^{\infty}(\Phi^{t}(\mathcal{M}))$ denote the pushforward of smooth functions on $\mathcal{M}$ by $\Phi^{t}$, i.e., $\Phi^{t}_{*}\varphi=\varphi\circ\Phi^{-t}$. Similarly define the pullback $(\Phi^{t})^{*}:C^{\infty}(\Phi^{t}(M))\to C^{\infty}(\mathcal{M})$ by $(\Phi^{t})^{*}\varphi=\varphi\circ\Phi^{t}$. We note that for volume-preserving $\Phi^{t}$, the pushforward $\Phi^{t}_{*}$ is the transfer (or Perron-Frobenius) operator and for general $\Phi^{t}$, the pullback $(\Phi^{t})^{*}$ is the Koopman operator.

We will say that regions in phase space remain coherent if their boundaries resist filamentation under the imposed nonlinear dynamics. Lack of filamentation is connected to mixing in an $L^{2}$ sense (see e.g. Corollary 4.4.1 [28] and the finite-time counterpart in Proposition 2 [14]) because long boundaries present large interfaces over which diffusive mixing can occur in the presence of small diffusion. Lack of filamentation is also consistent with fluid mixing norms for advective flows [30, 36], which are necessarily computed with Sobolev-type norms to mathematically introduce an effect akin to diffusion, as well as with recent variational approaches to Lagrangian coherent structures [26], which control evolved length of material curves.

2.1 The dynamic Laplacian with Neumann boundary conditions

The framework of [15] introduces the dynamic isoperimetric problem: The question is how to disconnect the given manifold $\mathcal{M}$ by a co-dimension 1 manifold $\Gamma$ in such a way that the average evolved size of $\Gamma$ is minimal relative to the volume of the two disconnected pieces. Let $\mathcal{T}\subset[0,t_{F}]$, $|\mathcal{T}|>1$, $0\in\mathcal{T}$, be the collection of times at which we wish to check the size of our evolving disconnector. We define the dynamic Cheeger constant for a smooth co-dimension 1 manifold $\Gamma\subset\mathcal{M}$ that disconnects $\mathcal{M}$ into two pieces $\mathcal{M}_{1},\mathcal{M}_{2}\subset\mathcal{M}$, as

[TABLE]

where $\ell_{d-1}$ is co-dimension 1 volume and $\ell$ is $d$-dimensional volume. The numerator of (1) is the average size of the evolved hypersurfaces $\Phi^{t}\Gamma$ over the time instants in $\mathcal{T}$. We wish to solve

[TABLE]

The minimising $\Gamma$ generates a family $\{\Phi^{t}\Gamma\}_{t\in\mathcal{T}}$ of material (evolving with the flow) disconnectors that are of minimal average size relative to the volume of the disconnected pieces of $\mathcal{M}$. Note that the minimisation problem (2) and its solution $\Gamma$ are frame-invariant because $\ell_{d-1}$ and $\ell$ are invariant under time-dependent translations and rotations.

A sharp relationship between $\mathbf{h}$ and functional minimisation is provided by the dynamic Federer-Fleming theorem [15], which states that $\mathbf{h}=\mathbf{s}$, where

[TABLE]

is the dynamic Sobolev constant. The same result for non-volume-preserving transformations on weighted Riemannian manifolds with general Riemannian metrics is developed in [18].

When replacing the $L^{1}$ optimisation with $L^{2}$ optimisation (which may be thought of as a regularisation), the problem can be solved exactly and efficiently as an eigenproblem for a dynamic Laplace operator:

[TABLE]

subject to homogeneous Neumann boundary conditions. Regularising the $L^{1}$ optimisation problem and computing the exact solution of the $L^{2}$ version of the problem is a standard approach in optimisation, and when evaluating the quality of solutions provided by level sets of the $L^{2}$ eigenfunctions, one can make this assessment using the original $L^{1}$ criterion (2). The minimum of this regularised problem is given by the eigenfunction of (4) corresponding to the second eigenvalue $\lambda_{2}$, and the optimal value of the $L^{1}$ minimisation (which equals $\mathbf{h}$) satisfies the bound $\mathbf{h}\leq 2\sqrt{-\lambda_{2}}$ (see Corollary 3.4 [15] and the discussion following Theorem 4.5 [18]).

The weak form of the eigenproblem

is developed by first integrating against a test function $\psi\in H^{1}(\mathcal{M})$, applying the transformation formula for $\Phi^{t}$ (which we firstly assume to be volume-preserving in the calculation below) and then applying Green’s formula, noting that the boundary conditions are the natural ones:

[TABLE]

Thus we obtain the weak form of the eigenproblem (4): Find eigenpairs $(\lambda,v)$ such that

[TABLE]

Note that we have completely eliminated any consideration of the boundary. If we want to have a continuous-time version of the above eigenproblem, we simply replace the sum over $t\in\mathcal{T}$ with an integral. A related discussion for the case of $T$ not volume preserving is given in the appendix.

2.2 The dynamic Laplacian with Dirichlet boundary conditions

In some situations, one may wish to restrict coherent sets to the class of submanifolds that do not intersect the boundary of $\mathcal{M}$ (which is assumed to be nonempty, otherwise one uses the setting of Section 2.1). We will illustrate one such possibility in Section 4.4, where an initial rectangular ocean domain becomes highly distorted. An appropriate isoperimetric problem in this setting is to search for an open submanifold $A\subset\mbox{int}(\mathcal{M})$ with compact closure, and with an evolving boundary $\Phi^{t}(\partial A)$ that remains small relative to the volume of $A$ for $t\in\mathcal{T}$. This can be achieved by modifying the expression (1) to

[TABLE]

and defining the dynamic (Dirichlet) Cheeger constant $\mathbf{h}_{d}$ by

[TABLE]

We will call the minimising $A\subset\mbox{int}(\mathcal{M})$ the dominant finite-time coherent set.

The corresponding version of the (Dirichlet) Sobolev constant is:

[TABLE]

where $C^{\infty}_{c}(\mathcal{M},\mathbb{R})$ denotes the set of non-identically vanishing, compactly supported, $C^{\infty}$ real-valued functions on $\mbox{int}(\mathcal{M})$. One can prove a dynamic (Dirichlet) Federer-Fleming theorem:

Theorem 1**.**

Let $\mathcal{M}$ be a compact subset of $\mathbb{R}^{d}$ with nonempty $C^{\infty}$ boundary, and $\mathbf{h}_{d}$ and $\mathbf{s}_{d}$ defined as in (7) and (8), respectively. Then $\mathbf{h}_{d}=\mathbf{s}_{d}$.

Proof.

In the setting of $T$ volume-preserving, this result follows from the obvious modifications to the proof of Theorem 3.1 [15]. ∎

The solution of the $L^{2}$ version of the $L^{1}$ optimisation problem (8) is given by the eigenfunction of the eigenproblem (9)–(10) corresponding to the largest eigenvalue $\lambda_{1}<0$.

[TABLE]

Note that in comparison to (4), we have simply added the Dirichlet condition (10), and thus the left-hand-side of (9) is the dynamic Laplace operator acting on a function $v$, and this operator is elliptic and self-adjoint as in the Neumann boundary condition case. In this Dirichlet case, we also have a dynamic Cheeger inequality relating $\mathbf{h}_{d}$ to $\lambda_{1}$.

Theorem 2**.**

Let $\mathcal{M}$ be a compact subset of $\mathbb{R}^{d}$ with nonempty $C^{\infty}$ boundary, and $\mathbf{h}_{d}$ be defined as in (7). Then $\mathbf{h}_{d}\leq\sqrt{-2\lambda_{1}}$, where $\lambda_{1}$ is the largest eigenvalue of (9)–(10).

Proof.

In the setting of $T$ volume-preserving, this follows from the obvious modifications to the proof of Theorem 3.2 [15]. ∎

The weak form of the eigenproblem (9)–(10) is identical to the weak form (5), except that we replace the space of test functions $H^{1}(\mathcal{M})$ with $H^{1}_{0}(\mathcal{M})$.

[TABLE]

In terms of heat flow (cf. [27]), one interprets the Dirichlet boundary conditions as “refrigeration” fixed at temperature zero on the boundary of $\mathcal{M}$. The eigenfunction $\phi_{1}$ corresponding to $\lambda_{1}$ is strictly positive (by convention) in the interior of $\mathcal{M}$ and $\phi_{1}$ is designed so that on average (over $t\in\mathcal{T}$) the rate of heat loss through the evolving boundaries $\Phi^{t}(\mathcal{M}),t\in\mathcal{T}$ is minimised. This is in contrast to the previous Neumann (no flux) boundary conditions of Section 2.1, where the second eigenfunction describes a signed heat distribution that approaches the constant equilibrium heat distribution at the slowest average rate (again averaged over $t\in\mathcal{T}$).

3 Finite Element discretization

We are going to use a Ritz-Galerkin discretization of the weak eigenproblem (5) and a finite element basis for the underlying ansatz and test function spaces. As mentioned in the introduction, one particular advantage of this approach is that if we choose the same basis for the ansatz and the test function space then we end up with symmetric stiffness and mass matrices, i.e. we inherit the self-adjointness of the dynamic Laplacian.

Let $V^{0}_{n}\subset H^{1}(\mathcal{M})$ be a finite-dimensional subspace with basis $\varphi_{1}^{0},\ldots,\varphi_{n}^{0}$. For brevity, we write $\varphi_{i}=\varphi^{0}_{i}$ in the following. The discrete eigenproblem corresponding to (5) is to find pairs $(\lambda,v)$, $\lambda\in\mathbb{C}$, $v=\sum_{i=1}^{n}u_{i}\varphi_{i}\in V^{0}_{n}$, such that

[TABLE]

where $u=(u_{1},\ldots,u_{n})^{T}\in\mathbb{C}^{n}$,

[TABLE]

are the stiffness and mass matrices. In the case where $\Phi^{t}$ is not volume-preserving and for initial mass distribution $\mu^{0}$ (which evolves to $\mu^{t}$ by push-forward: $\mu^{t}=\mu^{0}\circ\Phi^{-t}$), one uses

[TABLE]

(see Section 3.2 on how to approximate $\mu^{t}$).

For ansatz and test functions, we are going to use the usual piecewise linear nodal basis on a triangulation of some set of nodes $x_{1}^{0},\ldots,x_{n}^{0}\in\mathcal{M}$, i.e. triangular $P_{1}$ Lagrange elements. We provide a compact MATLAB implementation based on the code by Strang [35] in the Appendix and a faster one based on code from the iFEM package by Chen [7] on https://github.com/gaioguy/FEMDL.

We describe two distinct approaches: the first is based on transforming the weak eigenproblem (5) into an equivalent form which does not require integration on $\Phi^{t}(\mathcal{M})$. This approach, however, uses the pullback of the Euclidean metric which is described by the Cauchy-Green deformation tensor and therefore requires the evaluation of $D\Phi^{t}$ (i.e. the integration of the variational equation in case of an ordinary differential equation). In the second approach, we try to avoid evaluation of $D\Phi^{t}$. Then, however, quadrature on $\Phi^{t}(\mathcal{M})$ and the evaluation of the transfer operator will be required, but the former is relatively cheap and the latter can be even completely eliminated. To this end, we use ansatz/test function spaces $V^{t}_{n}\subset H^{1}(\Phi^{t}(\mathcal{M}))$ with bases $\varphi_{1}^{t},\ldots,\varphi_{n}^{t}$ on $\Phi^{t}(\mathcal{M})$ for each $t\in\mathcal{T}$.

3.1 The Cauchy-Green approach

Note that for general smooth $\Phi^{t}$ and smooth initial mass $\mu^{0}$ (in the case of volume-preserving $\Phi^{t}$, one may take $\mu^{0}=\mu^{t}=\ell$), we have

[TABLE]

where, $C_{t}:=(D\Phi^{t})^{\top}D\Phi^{t}$ is the (right) Cauchy–Green deformation tensor field. We thus have $\hat{D}_{ij}^{t}=\int_{\mathcal{M}}\nabla\varphi_{i}\bullet C_{t}^{-1}\nabla\varphi_{j}\;d\mu^{0}$ and can solve (12).

One advantage of this approach is that we do not need a numerical approximation of $\mu^{t}$. However, depending on the properties of the given transformation $\Phi^{t}$, the numerical evaluation of the integral (3.1) might be challenging as the tensor $C_{t}^{-1}$ might be of large variation locally. In all the following experiments, we employ a Gauss quadrature of varying degree. This approach has been applied to the numerical experiments in [27].

3.2 The transfer operator approach

Instead of approximating the pullbacks of the Euclidean metric under the flow map, we can instead approximate pushforwards of the basis functions $\varphi_{i}$ by $\Phi_{*}^{t}$. We do this by collocation in two ways:

Each manifold $\Phi^{t}(\mathcal{M}),t\in\mathcal{T}$, is meshed independently (using nodes that do not necessarily arise from trajectories) and thus the bases $\varphi^{t}_{i}$, $i=1,\ldots,n$ are not related for different $t\in\mathcal{T}$. 2. 2.

The mesh for each $\Phi^{t}(\mathcal{M})$, $t\in\mathcal{T}$, is a triangulation of $\{\Phi^{t}x_{i}\}_{i=1}^{n}$, namely, a triangulation of the images of the nodes $x_{i},i=1,\ldots,n$, of the mesh used for $\mathcal{M}$. As we will see, this construction completely removes any need to estimate the action of $\Phi^{t}_{*}$.

Suppose now that we have found coefficients $\alpha^{i}_{k}\in\mathbb{R}$ such that

[TABLE]

in a sense to be made precise in the next two subsections. We then have

[TABLE]

where $\alpha^{j}=(\alpha^{j}_{1},\ldots,\alpha^{j}_{n})^{\top}$ and $D^{t}=(D_{k\ell}^{t})_{k\ell}$.

Note that the stiffness matrix ${D}^{t}$ is symmetric and sparse as the $\varphi^{t}_{i}$ are locally supported; in particular, the time to compute $D^{t}$ scales linearly with $n$. In fact, the vectors $\alpha^{i}$ should also be sparse due to the continuity of $\Phi^{t}$ and the local support of the bases $\varphi^{t}_{1},\ldots,\varphi^{t}_{n}$. If so, then the matrices $\hat{D}^{t}$ will also be sparse and symmetric. Carrying out the above computations for each $t\in\mathcal{T}$, we then directly solve (12).

3.2.1 The non-volume-preserving case

In the case of non-volume-preserving $\Phi^{t}$ and general initial (probability) measure $\mu^{0}$ we obtain

[TABLE]

in (3.2). If we approximate $\mu^{0}$ by point masses on the nodes $x_{i}^{0}$, e.g. $\mu^{0}\approx\tilde{\mu}^{0}=\sum_{i=1}^{n}a_{i}\delta_{x_{i}^{0}}$, with $a_{i}=\int_{\mathcal{M}}\varphi_{i}^{0}\;d\mu^{0}$, then the pushforward of $\tilde{\mu}^{0}$ under $\Phi^{t}$ is exactly given by $\tilde{\mu}^{t}=\sum_{i=1}^{n}a_{i}\delta_{x_{i}^{t}}$. From $\tilde{\mu}^{t}$, we obtain a corresponding approximate density of $\mu^{t}$ as

[TABLE]

where $\ell_{i}^{t}=\int_{\Phi^{t}\mathcal{M}}\varphi_{i}^{t}\;d\ell$, which we can now use in order to approximate the entries of $D^{t}$:

[TABLE]

The outer sum in (19) over all elements $e^{t}$ from the triangulation of $\Phi^{t}(\mathcal{M})$ which are in the intersection of the supports of $\varphi_{k}$ and $\varphi_{\ell}$. Further,

[TABLE]

For the case of piecewise linear triangular elements, the gradients of the $\varphi_{j}$ are constant on each element and so in this case we get

[TABLE]

We now briefly compare (3.2.1) to the integral in the right hand side of (3.2): Note that in the volume-preserving case $a_{i}=\ell_{i}$ and the summand in the RHS of (3.2.1) becomes $\ell_{i}/\ell^{t}_{i}\approx 1$ for sufficiently fine meshes. If $\Phi^{t}$ is affine on the support of $\varphi_{i}^{0}$, then $\ell_{i}/\ell^{t}_{i}=1$; with increasing nonlinearity of the map, a good estimation requires finer mesh elements. The factor $\text{area}(e^{t})\;\frac{1}{3}$ also appears in (3.2) and thus (3.2.1) approximately reduces to (3.2) in the volume-preserving case.

Naturally the symmetry and sparseness properties of the resulting $\hat{D}^{t}$ are retained. We will see in sections 3.2.2 and 3.2.3 that this is the only change required; in particular, the expressions for $\alpha$ do not change at all for non-volume-preserving $\Phi^{t}$ and general initial measure $\mu^{0}$.

The next two subsections describe two different approaches to computing the $\alpha^{i}_{k}$.

3.2.2 Collocation on non-adapted meshes

In this variant, the nodes $x^{t}_{i}$ at time $t\in\mathcal{T}$ are in general unrelated to $x^{s}_{i}$ for some other time $s\in\mathcal{T}$, in particular $\{x^{t}_{i}\}_{t\in\mathcal{T}}$ need not be on a trajectory of $\Phi^{t}$. We wish to obtain the matrix $\alpha$ by collocation. To this end, in the case where $\Phi^{t}$ is volume-preserving, using (16) we should solve the following problem:

[TABLE]

where the $\delta_{x^{t}_{m}}$ are Dirac deltas. Since we are using a nodal basis $\varphi^{t}_{1},\ldots,\varphi^{t}_{n}$ at each $t\in\mathcal{T}$, we have $\langle\varphi^{t}_{k},\delta_{x^{t}_{m}}\rangle=\varphi^{t}_{k}(x^{t}_{m})=\delta_{k,m}$ (Kronecker delta) and so the right hand side of (21) is simply $\alpha_{m}^{i}$. Thus we obtain the explicit expression

[TABLE]

It is clear that $\alpha^{i}_{k}$ is nonnegative and summing the right hand side of (22) over $i$ we see that $\sum_{i}\alpha^{i}_{k}=1$ for each $k=1,\ldots,n$; thus the matrix $\alpha=(\alpha^{i}_{k})_{ik}$ is column-stochastic.

In the case of initial probability measure $\mu^{0}$ and non-volume-preserving $\Phi^{t}$, let $h_{\mu^{t}}=d\mu^{t}/d\ell$ be the density of $\mu^{t}$ with respect to Lebesgue. The inner products in (21), now weighted according to $\mu^{t}$ become

[TABLE]

Cancelling $h_{\mu^{t}}(x^{t}_{m})$ on both sides, we again obtain the simple formula (22) for $\alpha^{i}_{m}$.

3.2.3 Collocation on adapted meshes

We again compute the matrix $\alpha$ by collocation, but instead of choosing a mesh at $t\in\mathcal{T}$ independently of the initial mesh we instead define the nodes of the mesh at time $t$ by $x_{i}^{t}=\Phi^{t}x_{i}^{0}$, namely the images of the nodes at the initial time. The expression (22) now becomes

[TABLE]

thus $\alpha$ is the $n\times n$ identity matrix.

This approach has several major advantages, particularly when one does not have a full model of the dynamical system, but instead only has access to $n$ trajectories $\{x_{i}^{t}\}_{t\in\mathcal{T}},i=1,\ldots,n$, although the approach is also very effective even with a full model which can generate trajectories. The algorithm is as follows.

Algorithm 1: (Input – trajectories $\{x_{i}^{t}\}_{t\in\mathcal{T}},i=1,\ldots,n$; Output – approximate eigenfunctions of dynamic Laplacian)

for each $t\in\mathcal{T}$,

1. (a)

triangulate the set of nodes $\{x_{1}^{t},\ldots,x_{n}^{t}\}$, e.g. by a Delaunay triangulation or, in case of highly nonlinear $\Phi^{t}$, using some alpha complex of the nodes 444A Delaunay triangulation always yields a triangulation of the convex hull of a given set of points. In the case of highly nonlinear $\Phi^{t}$, this may lead to large triangles connecting nodes which a far apart. In contrast, alpha complexes [11] yield triangulations with triangles being constructed only locally. In MATLAB, alpha complexes can be computed using alphaShape/alphaTriangulation.. 2. (b)

compute $D^{t}$ using (3.2) or (18) as appropriate, 2. 2.

compute $M$ using (13) or (14) as appropriate, 3. 3.

form $\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}D^{t}$ and solve the eigenproblem (12).

The formation of the triangulations and the computation of each $D^{t}$ and the single $M$ are very fast. In fact, in our experiments, solving the sparse, symmetric eigenproblem (12) for the eigenvalues near to zero, while also fast, is the most time-consuming step. We emphasise that no computation of the transfer operator is required as it is implicit in the trajectory information.

3.2.4 Computing with incomplete trajectory data

We now consider the situation where we have trajectory data $\{x_{i}^{t}\},i=1,\ldots,n,t=1,\ldots,T$, but some of these data points are “empty” (we don’t have the position data at certain times). We discuss how to modify the method of section 3.2.3 to handle this situation. Let $\mathcal{T}_{i}=\{t\in\mathcal{T}:x_{i}^{t}\mbox{ exists}\}$ and $I_{t}=\{i\in\{1,\ldots,n\}:x_{i}^{t}\mbox{ exists}\}$. Respectively, these sets are the times at which we have positions for particle $i$, and indices of those particles for which we have positions at time $t$. We will need to compare distances between trajectories, so we assume that $\mathcal{T}_{i}\cap\mathcal{T}_{j}\neq\emptyset$ for all pairs $1\leq i,j\leq n$; that is, for each point pair $(i,j)$ there exists at least one time $t$ (depending on the pair $(i,j)$) for which $i,j\in I_{t}$. Note that this in particular includes situations where there is no single time at which all data pairs $i,j$ are available.

At each $t\in\mathcal{T}$ we triangulate the points $\{x_{i}^{t}:i\in I_{t}\}$. Using the corresponding basis $\{\varphi_{i}^{t}:i\in I_{t}\}$, we compute ${D}^{t}$ as

[TABLE]

For the mass matrix $M$, first note that

[TABLE]

Motivated by this, we set

[TABLE]

and use

[TABLE]

The rest of the calculations proceed as in the previous sections.

In the non-volume-preserving case, replace $d\ell$ in (25) and (26) at index $t$ with $d\mu^{t}$. We described how to compute the entries in $D^{t}$ in this case in Section 3.2.1. For $i,j\in I_{t}$ the entries of the time-dependent mass matrices can be approximated by

[TABLE]

and the integrals $\int_{\Phi^{t}(\mathcal{M})}\varphi_{i}\;\varphi_{j}\;\varphi_{k}\;d\ell$ can be computed exactly, using the fact that the $\varphi_{i}$ are affine on each triangle.

Algorithm 2: (Input – incomplete trajectories $\{x_{i}^{t}\}_{t\in\mathcal{T}},i=1,\ldots,n$; Output – approximate eigenfunctions of dynamic Laplacian)

for each $t\in\mathcal{T}$,

1. (a)

triangulate the set of nodes $\{x_{i}^{t}:i\in I_{t}\}$ (Delaunay or alpha complex), 2. (b)

compute ${D}^{t}$ as per (25) or the equivalent version of (18), as appropriate, 3. (c)

compute $M^{t}$ as per (26) or (28), as appropriate, 2. 2.

form $\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}D^{t}$ and use $M$ from (27) to solve the eigenproblem (12).

3.2.5 Discussion of the three approaches

The three variants of the FEM approach to the eigenproblem (4) in the preceeding sections come with different advantages/disadvantages which we now briefly discuss.

The Cauchy-Green approach from Section 3.1 only requires a single triangulation of the underlying manifold $\mathcal{M}$, avoiding potential complications with a complicated geometry of the image manifold(s) $\Phi^{t}(\mathcal{M})$. Additionally, this leads to fewer nonzero entries in the stiffness and mass matrices in comparison to the other two methods, so that the eigenproblem is solved faster (cf. the experiments in the next section). It exploits higher order information on $\Phi^{t}$ via the Cauchy-Green tensor, although this comes at a higher computational cost since $D\Phi^{t}$ has to be computed at every quadrature node (of which typically several are required within each mesh element). Also, for highly nonlinear flows, $D\Phi^{t}$ tends to vary very rapidly in space, which sometimes renders the quadrature on the elements challenging (cf. the experiment on the rotating double gyre in Section 4.1).

In the transfer operator approach (Section 3.2.2), we may use a single triangulation only in the case that $\mathcal{M}=\Phi^{t}(\mathcal{M})$, since then we can choose the same basis of $V_{n}^{t}$ for each $t\in\mathcal{T}$. On the other hand, a major advantage is the removal of the need to compute the derivative of the flow map. Instead, we need to evaluate the inverse flow map, which is possible if $\{x^{t}_{m}\}_{t\in\mathcal{T}}$ arises as a trajectory, or in other settings where the preimages are known or can be easily computed. The matrices $D_{kl}^{t}$ and $(\alpha_{k}^{i})_{ik}$ are sparse with the number of nonzero entries per row/column determined by how many elements support a basis function (six in the case of triangular elements in 2D), although their product $\hat{D}^{t}$ may have considerably more nonzero entries than $D_{kl}^{t}$ and $(\alpha_{k}^{i})_{ik}$, which would make solving the eigenproblem more expensive.

In the adaptive transfer operator approach (Section 3.2.3), we need to construct a separate triangulation for every $t\in\mathcal{T}$ (even in the case $\mathcal{M}=\Phi^{t}(\mathcal{M})$). For each of our examples we find $|\mathcal{T}|=2$ already yields good results. There is no estimation of the transfer operator through the $(\alpha_{k}^{i})_{ik}$, so the stiffness matrix calculation is completely standard for FEM. Further, the only input data are the positions of the nodes at each time instance, which, e.g., is readily available in most experimental or real-world settings. The number of nonzeros of the overall stiffness matrix $\sum_{t}\hat{D}^{t}$ grows linearly with $|\mathcal{T}|$, so that the solution of the eigenproblem becomes more expensive.

3.3 Extraction of (multiple) coherent sets from eigenfunctions

By construction of the dynamic Laplacian, and the theory in Section 2 and the references therein, an optimal 2-partitioning of the domain into 2 coherent sets is given in terms of level sets of the second eigenfunction (often corresponding to values near zero). Extending this idea heuristically to multiple coherent sets, one searches for multiple regions defined by several level sets. A straightforward way to extract these is by computing an optimal level set for each eigenfunction, cf. [15, 17]. In general, this does not guarantee that the resulting coherent set is connected, e.g. it might be the union of two vortices in a fluid flow. This is similar to the situation in the context of time-independent systems where one is interested in almost-invariant sets, methods based on the zero level set [9] (the sign structure) or other thresholds [13] of the relevant eigenvectors.

Various approaches for the extraction of “individual” coherent sets have been proposed: For almost-invariant sets, a least-squares assignment strategy [10] have been proposed. Later, techniques drawn from spectral graph partitioning [6, 2] based on clustering embedded eigenvectors via fuzzy c-means clustering [16, 13] have been proposed. In our examples, we use a variant of this latter approach by employing MATLAB’s kmeans command.

In many cases, gaps in the spectrum of the dynamic Laplacian can be used in order to obtain a heuristic on the number of coherent sets: If there is a large gap after the $k$-th eigenvalue, then one should search for $k-1$ connected coherent sets555for Neumann boundary conditions, since the constant eigenfunction at $\lambda=0$ does not contribute.. This heuristic is motivated by perturbation arguments for the spectra of linear operators, cf. [9, 13, 27].

3.4 Link with graph-based methods

The transfer operator approach of §3.2.3 can be interpreted in terms of graph Laplacians. For any $t$ the matrix $\hat{D}^{t}$ from (13) or (14) is clearly symmetric, and also has zero column and row sums. For the latter property, note that summing both sides in (13) over $j$ yields

[TABLE]

as $\sum_{j}\varphi_{j}\equiv\mathbf{1}$ by construction (and $\Phi^{t}_{*}\mathbf{1}=\mathbf{1}$ and $\nabla\mathbf{1}=\mathbf{0}$).

Moreover, the diagonal entries of $\hat{D}^{t}$ are clearly positive. It is known (see Proposition 3 [38], also p5 [1]) that the off-diagonal entries of $\hat{D}^{t}$ are non-positive for two-dimensional Delaunay meshes and for three-dimensional Delaunay meshes when all dihedral angles of the tetrahedra are less than $\pi/2$. Indeed, in all of our numerical experiments, including in three dimensions, we observe all interior triangles/tetrahedra yield non-positive entries. In view of the above, we proceed under the assumption that all off-diagonal entries of $\hat{D}^{t}$ are all non-positive for the purposes of relating our numerical eigenproblem to a graph Laplacian eigenproblem.

Consider the graph $G^{t}$ which is our triangulation of images of data points at time $t$. Because of the symmetry, zero column sum, and sign structure of $D^{t}$ discussed above, we can write $\hat{D}^{t}=\Pi^{t}-W^{t}$, where

[TABLE]

are the weights assigned to edges $(i,j)$ in our graph $G^{t}$, and $\Pi^{t}_{ii}=\hat{D}^{t}_{ii}$ is a diagonal matrix with the $i^{\rm th}$ diagonal element equal to the degree of node $i$ in $G^{t}$, namely $\Pi^{t}_{ii}=-\sum_{j}\hat{D}^{t}_{ij}$. Thus, we can immediately interpret $D^{t}$ as an (unnormalized) graph Laplacian [8].

We now interpret these off diagonal arc weights in terms of the trajectory dynamics. Suppose that at time $t$, the positions of two nearby, distinct, evolved data points $\Phi^{t}(x_{i}),\Phi^{t}(x_{j})$ are sufficiently close that $(i,j)$ is an arc in $G^{t}$ (this means $\Phi^{t}(x_{i}),\Phi^{t}(x_{j})$ are sufficiently close that the supports of $\varphi^{t}_{i}$ and $\varphi^{t}_{j}$ intersect). Then the corresponding entry of $\hat{D}^{t}_{ij}$ will be negative and will decrease in magnitude as the distance between these points increases. This is because the integrand in (13) or (14) is proportional to a product of two approximately linearly decreasing gradients, which outweighs the approximately linear increase in the integration domain given by the mesh volume. If we continue to increase the distance between $\Phi^{t}(x_{i}),\Phi^{t}(x_{j})$, eventually the arc $(i,j)$ will no longer exist in $G^{t}$ because the Delaunay mesh will find a more efficient way to mesh the points. In summary, we see that the edge weights given by $\hat{D}^{t}_{ij}$ reflect distance between nearby points $\Phi^{t}(x_{i}),\Phi^{t}(x_{j})$.

When we solve the problem $-\frac{1}{|\mathcal{T}|}\sum_{t\in\mathcal{T}}\hat{D}^{t}y=\lambda My$, we are averaging the graph Laplacians, normalized by the symmetric, nonnegative mass matrix $M$. There are two differences in this normalization to the standard graph Laplacian normalization, which is based simply on node degree. Firstly, $M$ is not diagonal, although in any given row or column, the largest entry will lie on the diagonal. Moreover, the small number of off-diagonal entries of $M$ coincide with the off-diagonal entries of $\hat{D}^{0}$ and correspond to arcs $(i,j)$ that are in the graph $G^{0}$. Returning to the interpretation of off-diagonal entries as distances $\hat{D}^{t}_{ij}$, we see that normalizing by the mass matrix automatically handles nonuniformly distributed data, because if initial points $x^{0}_{i},x^{0}_{j}$ are far apart the value of $M_{ij}$ will be commensurately larger; see Figure 8 for a dramatic illustration of this automatic handling. This is in contrast to e.g. the approaches of [21, 24], where only the linear distance between points is taken into account, and not additionally the volume in phase space “represented” by those points. This is important in the case where the data is not uniformly distributed in space, such as e.g. trajectory data from ocean drifters.

4 Numerical experiments

In all subsequent experiments, we employ standard linear triangular/tetrahedral elements (i.e. $P_{1}$ Lagrange elements) and an explicit Runge-Kutta scheme of order (4,5) with adaptive step size control as implemented in MATLAB’s ode45 integrator, using the relative and absolute error tolerance $10^{-3}$. The derivative of the flow map is approximated using central finite differences.

4.1 Experiment: The rotating double gyre

We start with the rotating double gyre flow introduced in [32], a Hamiltonian system with Hamiltonian $H=-\psi$, where $\psi$ is the stream function

[TABLE]

and transition function

[TABLE]

On the square $\mathcal{M}=[0,1]^{2}$, the vector field exhibits two gyres with centers at $(\frac{1}{2},\frac{1}{2})$ and $(\frac{3}{2},\frac{1}{2})$ initially (at $t=0$) which rotate by $\pi/2$ during the flow time $T=1$.

Cauchy-Green approach.

We use $\mathcal{T}=\{0,1\}$ and employ Gauss quadrature of degree 5 for the numerical integration of (3.1). A lower degree yields degraded results. Indeed, as shown in Fig. 1, the (inverse) Cauchy-Green tensor is of quite large variation (cf. our earlier remarks).

Figure 2 shows the spectra on a triangulation of $25\times 25$ equidistant points (1152 triangles) as well as on the Delaunay triangulation of a set of 625 scattered nodes (1236 triangles). There appears to be a gap after the fourth eigenvalue. The method of Section 3.2, however, has a gap after the third eigenvalue for both the regular and scattered data, and for this reason and to enable comparison with prior work [20], we used the leading 3 eigenvectors in all cases. Figures 3 and 4 show the 2nd and 3rd eigenvectors as well as the resulting decomposition into coherent sets obtained by kmeans-clustering of these eigenvectors as described in Section 3.3. Already at this comparatively low resolution, the results are broadly similar to those in [20]. Note, however, that we need to evaluate $C_{t}^{-1}$ (for which we evaluate $D\Phi^{t}$ by finite differencing) at roughly $8000$ points, which takes around 2 seconds. The assembly of the matrices takes 0.02 s and the solution of the eigenproblem 0.05 seconds666All timings measured on a dual core Intel i5 with 2.4 GHz and MATLAB R2016a..

Transfer operator approach.

We repeat the same experiment with $\mathcal{T}=\{0,1\}$ using the approach from Section 3.2.2. Figure 5 (left) shows the spectrum (with a small gap after the third eigenvalue), Figure 6 the second and third eigenvectors as well as the resulting coherent 3-partition. Here, the inverse flow map (without variational equation) had to be evaluated $25\cdot 25=625$ times only, which takes 0.1 s, the assembly of the matrices takes 0.02 s and the solution of the eigenproblem again 0.1 seconds. Note that the only data that we input to our method is the initial and final positions of the 625 nodes.

Adaptive transfer operator approach.

Using the approach from Section 3.2.3 on this experiment, we obtain the spectrum in Figure 5 (right), with a gap after the third eigenvalue, and the corresponding eigenvectors and clustered coherent sets in Figure 7. As for the (non-adaptive) transfer operator approach, the flow map (without variational equation) had to be evaluated $25\cdot 25=625$ times only, leading to the same computation times.

In order to show that this approach is robust with respect to non-uniformly sampled data, we repeat this experiment with a highly non-uniform triangulation: We randomly choose 20000 points in $[0,1]\times[0,0.5]$ and 200 points in $[0,1]\times[0.5,1]$ (both according to a uniform distribution) and use the Delaunay triangulation (cf. Figure 8) of the union of these two point sets. The resulting eigenvectors and coherent 3-partition are almost identical to the one obtained by a uniform point sampling.

Missing data.

We finally analyze how well the transfer operator approach from Section 3.2 performs when some of the trajectory data is missing: To this end, we again use 625 randomly scattered points, but now with a denser sampling in time, namely $\mathcal{T}=\{0,0.2,0.4,0.6,0.8,1\}$. Figure 9 (top row) shows the result for the adaptive transfer operator method, which is essentially the same as in Figs. 6-8.

We now randomly delete 60% of the data, reducing the data from 625 points to 250 data points per time step (and the total number of data points is the same as in Figs. 6 and 7). Repeating the same experiment using the approach for missing data described in Section 3.2.4 we obtain the results shown in the bottom row of Figure 9 which still clearly show the relevant structures. The corresponding spectrum is shown in Fig. 10, the gap after the third eigenvalue is still visible.

4.2 Experiment: A non-volume preserving standard map and nonuniform initial density

For an experiment with a non-volume preserving map and a non-uniform initial measure $\mu^{0}$ on $\mathcal{M}$, we use a test case from the literature [18]. We consider the map $\Phi=\Phi_{2}\circ\Phi_{1}$ on the flat 2-torus $\mathcal{M}=2\pi S^{1}\times 2\pi S^{1}$ with

[TABLE]

The map $\Phi_{1}$ compresses horizontally at $x=\pi/2,3\pi/2$ and expands horizontally at $x=0,\pi$. The map $\Phi_{2}$ is the volume-preserving standard map with parameter 8. The initial measure $\mu^{0}$ has density $h^{0}(x,y)=\frac{1}{8\pi^{2}}\left(\sin(y-\frac{\pi}{2})+1\right)$, shown in Figure 11 (upper left); the ratio of its maximum to minimum values is infinite as the density dips to zero at $y=0$. This highly varying density and strongly nonlinear $\Phi$ make this example particularly challenging.

We used the adaptive transfer operator approach from Section 3.2.3 and the formulae for the non-volume preserving case from Section 3.2.1 with $\mathcal{T}=\{0,1\}$ (one iteration of $\Phi$) on a regular triangulation with $40\times 40$ nodes. We provide an implementation of scattered node triangulations on a periodic domain (which is not standard in MATLAB) at https://github.com/gaioguy/FEMDL.

The results are shown in Fig. 11, which should be compared to Fig. 4 in [18]. The boundaries of the optimal coherent sets in these figures are the cyan level sets, which approximately follow lines of constant $x+y$ value and constant $x$ value, respectively (as was the case in Fig. 7 of [15]). However, the distortion from $\Phi_{1}$ and the non-uniform weight $\mu_{0}$ causes deviations of these level sets from straight lines (see Section 6.4 [18] for further discussion). Note that the eigenvectors are captured reasonably well in Figure 11 (upper right, lower right), even though the approximation to the image density is apparently rather crude (Figure 11 (lower left)).

4.3 Experiment: The Bickley Jet

As a second example, we consider the Bickley jet flow introduced in [34], modeling a meandering jet with additional rotating vortices. The defining stream function for this Hamiltonian system is given by

[TABLE]

with parameter values $U_{0}=62.66\text{ ms}^{-1}$, $L_{0}=1770$ km and $A_{1}=0.0075,A_{2}=0.15,A_{3}=0.3,c_{1}=0.1446U_{0},c_{2}=0.205U_{0},c_{3}=0.461U_{0},k_{1}=2/r_{e},k_{2}=4/r_{e},k_{3}=6/r_{e},r_{e}=6371$ km as in [24]. We consider the associated flow on the domain $\mathcal{M}=[0,20]\times[-3,3]$ with periodic boundary conditions in $x$-direction on the time interval $t\in[0,40]$ days.

Cauchy-Green approach.

Employing $\mathcal{T}=\{0,40\}$ and the Delaunay triangulation on a regular grid of $100\times 30$ points as well as Gauss quadrature of degree 1, we obtain the results shown in Figure 12. The evaluation of $C_{t}^{-1}$ on roughly 6000 quadrature nodes takes about 8 s, the assembly of the matrices 0.07 s and the solution of the eigenproblem 0.4 seconds.

Adaptive transfer operator approach.

Using the same triangulation as for the Cauchy-Green approach, we obtain the results shown in Figure 13. The spectrum in Figure 13 (left) shows a gap after the 2nd eigenvalue (corresponding to the upper/lower separation shown in Figure 13 (top right)), followed by a further six eigenvalues before the next gap (corresponding to the 6 vortices present in the flow). The evaluation of the flow map on the $3000$ nodes of the triangulation takes 1 s, the assembly of the matrices 0.1 s and the solution of the eigenproblem 0.5 seconds.

Missing data.

In order to analyze the performance when some of the data is missing we use the same $100\times 30$ grid, but now with ten intermediate time samplings of the trajectories, i.e. $\mathcal{T}=\{0,4,8,\ldots,40\}$. The result of the adaptive transfer operator approach is essentially the same as in Fig. 13. We then randomly delete 80% of the data, reducing it from 3000 points to 600 points per time step, which in total yields the same number of data points as with $\mathcal{T}=\{0,40\}$ and no missing data. The approach from Section 3.2.4 leads to the results shown in Fig. 14, still clearly showing similar structures to Fig. 13. Furthermore, the spectral gap structure shown in Fig. 14 (left) remains broadly consistent with the results using 3000 points shown in Fig. 13 (left), with an additional separation between eigenvalues 6 and 7, possibly due to a loss of coherence information from the data destruction.

4.4 Experiment: Ocean flow from satellite data

We now consider an unsteady velocity field derived from AVISO satellite altimetry measurements in the domain of the Agulhas leakage in the South Atlantic Ocean. Under the assumption of geostrophicity, the sea surface height $h$ yields a stream function for the velocity of water at the surface and the corresponding equations of motion for particle trajectories are

[TABLE]

where $\varphi$ is the longitude and $\theta$ the latitude of a particle, $A(\theta)=g/{R^{2}2\Omega\sin\theta\cos\theta}$, $g$ is the gravitational constant, $R$ is the mean radius of the earth and $\Omega$ the earth’s mean angular velocity. We choose the same spatial and temporal domain as in [23], namely $[-4,6]\times[-34,-28]$ and a period of 90 days, from $t_{0}=$ Nov 11, 2006 on. Figure 15 shows the forward finite time Lyapunov exponent field of the data.

After a flow duration of 90 days the initial rectangular domain shown in Figure 15 becomes strongly distorted and filamented. In this situation, we are not interested in coherent sets that intersect this heavily filamented boundary and so we restrict to coherent sets in the interior of the domain. This can be achieved by replacing the (natural) homogeneous Neumann boundary conditions with homogeneous Dirichlet boundary conditions. In terms of numerical eigenvector computations, to apply homogeneous Dirichlet boundary conditions, one simply sets all rows and columns of $\hat{D}^{t}$ and $M$ in (13) corresponding to nodes on the boundary of the domain to rows and columns of zeros.

Cauchy-Green approach.

As in [23], we use a uniform grid of $250\times 150$ points, the associated Delaunay triangulation, yielding roughly 73000 triangles and Gauss quadrature of degree 1 (i.e. one quadrature point per element). The evaluation of $C_{t}^{-1}$ takes ca. 46 s, the assembly of the matrices 0.3 s and the solution of the eigenproblem 3 seconds.

Figure 16 (left) shows the largest 10 eigenvalues of the dynamic Laplacian with $\mathcal{T}=\{t_{0},t_{0}+90\}$, the results do not change significantly when using more intermediate time steps. Note that there appears to be a gap after the third eigenvalue. Figure 17 (left column) shows the associated first three eigenvectors. The corresponding coherent sets (as identified by k-means clustering the first three eigenvectors, cf. Fig. 18, left) nicely agree with those identified in [23], Fig. 9 (the rightmost vortex from Fig. 9 in that work appears in a lower eigenfunction which we do not show here).

We deliberately chose the same resolution as in [23] here, the results for the leading three eigenvectors do not change significantly, however, when decreasing the grid resolution down to $100\times 60$ points only.

Adaptive transfer operator approach.

The approach from Section 3.2.3 yields comparable results as shown in the right columns of Fig. 16. The computation times are: time integration: 4 s, assembly: 1.4 s, solution of the eigenproblem: 16 seconds. We can decrease the resolution down to a $150\times 90$ grid before the results start to deteriorate significantly.

Missing data.

We finally analyze how well the coherent sets can be recovered in the ocean flow experiment when data is missing. Again, we use the same set of nodes and $|\mathcal{T}|=10$ equidistant intermediate time steps. We then randomly delete 70% of the data, yielding roughly 12000 data points per time step. For comparison, in the experiment with full data and two time steps above we used around 37500 nodes per time step. The result shown in Figure 20 is qualitatively the same as the one with full data in Figure 18. We note, however, that when we delete even more data, these coherent sets will not be recovered as clearly any more.

4.5 Experiment: The unsteady ABC flow (3D)

For a 3D experiment, we consider the unsteady ABC flow, cf. [25], given by

[TABLE]

on the 3-torus, with parameter values $A=\sqrt{3},B=\sqrt{2},C=1$ on the time interval $t\in[0,1]$.

Cauchy-Green approach.

We employ a Delaunay triangulation on a regular grid of $25\times 25\times 25=15625$ points, yielding about 83000 tetrahedra, Gauss quadrature of degree 1, i.e. one quadrature point per tetrahedron and two time steps, i.e. $\mathcal{T}=\{0,1\}$. The integration of the variational equation takes 11 s, the assembly of the matrices 0.9 s and the solution of the eigenproblem 9 seconds. Figure 21 (left) shows the spectrum of the discrete dynamic Laplacian, Figure 22 (top) the 2nd (left) and 3rd (right) eigenvector.

Transfer operator approaches.

Using the same parameters as for the Cauchy-Green approach, we obtain the spectra shown in Figure 21 (center and right) and the eigenvectors in Figure 22 (center and bottom). This yields qualitatively the same results as the Cauchy-Green approach, the eigenfunctions from the non-adaptive transfer operator approach, however, appear to be less smooth. This seems to be due to the relatively coarse approximation of the transfer operator. Computation times are: non-adaptive: time integration 0.4 s, assembly 0.6 s, computation of the $\alpha$-matrix 1.4 s, eigenproblem 160 s; adaptive: integration 0.4 s, assembly 1.4 s, eigenproblem 40 s. Note that the increased computation time for solving the eigenproblem is due to the fact that the stiffness matrix has more nonzero entries ($10\times$ as many for the non-adaptive and twice for the adaptive TO approach in comparison to the CG approach, cf. the remark in Section 3.2.2 on the sparseness of the stiffness matrix in this case).

Acknowledgements

We thank Daniel Karrasch for helpful comments and contributions to the code, and Michael Feischl for helpful comments. GF thanks the Department of Mathematics at the Technical University Munich for hospitality during a visit in May 2016, and the John von Neumann Professorship scheme for financial support for this visit. GF was supported by an ARC Future Fellowship. OJ was supported by the Priority Programme SPP 1881 Turbulent Superstructures of the Deutsche Forschungsgemeinschaft.

Appendix A The non volume-preserving case

We now briefly sketch the theory behind the case of the underlying flow $\Phi^{t}$ not being volume-preserving, where we wish to track coherent masses according to some smooth initial mass distribution on $\mathcal{M}\subset\mathbb{R}^{d}$ given by a probability measure $\mu^{0}$ ($\mu^{0}=m$ would correspond to volume), and where the domain is possibly curved with the curvature described by a Riemannian metric. The measure $\mu^{0}$ is evolved by $\Phi^{t}$ according to $\mu^{t}:=\mu\circ\Phi^{-t}$. All computational aspects for data embedded in Euclidean space are described in Section 3.2.1.

The expressions (1) and (3) can be naturally extended to cover this situation. Firstly (1), where the size of the disconnector $\Gamma$ is computed according to the evolved measure $\mu^{t}$:

[TABLE]

where $\mu_{d-1}^{t}$ is the induced measure on $d-1$-dimensional surfaces at time $t$. Secondly, (3),

[TABLE]

where the subscripts $m^{t}$ denote that the computations are taken with respect to the Riemannian metric $m^{t}$ on $\Phi^{t}(\mathcal{M})$, which in many cases will be either the Euclidean metric, or if $\mathcal{M}$ is of dimension lower than $d$, the induced metric arising from the Euclidean metric. See [18] for details.

Finally, we replace (4) with

[TABLE]

where $\Delta_{\mu^{t}}$ is the $\mu^{t}$-weighted Laplace operator on $\Phi^{t}(\mathcal{M})$ (see [18] for details).

The weak form of the eigenproblem (34) can be written as [18]

[TABLE]

Again, if we want to have a continuous-time version of the above eigenproblem, we simply replace the average over $t\in\mathcal{T}$ with an integral.

Appendix B Appendix: Code example

We here provide sample code which performs the computations reported on in Section 4.1. Note that this code has been stripped down for readability. A more efficient version can be downloaded from https://github.com/gaioguy/FEMDL.