Higher-Order Finite Element Approximation of the Dynamic Laplacian

TL;DR

This paper develops and analyzes higher-order finite element methods for approximating the dynamic Laplacian, which helps identify coherent structures in fluid flows, with proven convergence and an accompanying Julia implementation.

Contribution

It introduces higher-order discretization schemes for the dynamic Laplacian and proves their convergence, enhancing numerical accuracy for identifying coherent structures.

Findings

Higher-order schemes improve approximation accuracy.

Eigenvalues and eigenvectors converge under certain conditions.

Efficient Julia implementation is provided.

Abstract

The dynamic Laplace operator arises from extending problems of isoperimetry from fixed manifolds to manifolds evolved by general nonlinear dynamics. Eigenfunctions of this operator are used to identify and track finite-time coherent sets, which physically manifest in fluid flows as jets, vortices, and more complicated structures. Two robust and efficient finite-element discretisation schemes for numerically computing the dynamic Laplacian were proposed in Froyland & Junge (2018). In this work we consider higher-order versions of these two numerical schemes and analyse them experimentally. We also prove the numerically computed eigenvalues and eigenvectors converge to the true objects for both schemes under certain assumptions. We provide an efficient implementation of the higher-order element schemes in an accompanying Julia package.

Equations

$${\Delta }^{dyn}:=\frac{1}{|\mathcal{I}|}\sum _{t\in \mathcal{I}}{T}_t^{\ast }{\Delta }_t{T}_{t,\ast }$$

$$a(u,v)$$

$${\Delta }^{dyn}:=\frac{1}{|\mathcal{I}|}{\int }_{\mathcal{I}}{T}_t^{\ast }{\Delta }_t{T}_{t,\ast }dt,$$

$$a(u,v)=\lambda ⟨u,v{⟩}_{L^2}\text{mbox}forallv\in S.$$

$$D\text{u}=\lambda M\text{u}.$$

$$\frac{{\lambda }_{i,h}-{\lambda }_i}{{\lambda }_i}\le C{h}^{2k}(h\to 0).$$

$${\Vert v_i-v_{i,h}\Vert }_{H^1}\le C{h}^k(h\to 0).$$

$${\Vert v_i-v_{i,h}\Vert }_{L^2}\le C{h}^{k+1}(h\to 0).$$

$$e_{i,h}:=\sqrt{1-|⟨v_i,v_{i,h}{⟩}_{L^2}{|}^2}.$$

$$T(x,y)=(x+y+a\sin (x),y+a\sin (x))$$

$$\dot{x}(t)$$

$$\dot{y}(t)$$

$$\psi (x,y,t)$$

$$a(u,v)=\frac{1}{2}\left(a^0(u,v)+a^1(Tu,Tv)\right)$$

$${\tilde{a}}_h(u,v):=\frac{1}{2}\left(a^0(u,v)+a^1(I_{h}Tu,I_{h}Tv)\right)$$

$$({\tilde{D}}_h{)}_{ij}$$

$$(M_h{)}_{ij}$$

$${\Vert I_{h}v-v\Vert }_{H^1}\to 0\text{as }h\to 0.$$

$$\underset{h\to 0}{\mathrm{limsup}}{\tilde{\lambda }}_{h,1}\le {\lambda }_1$$

$$\underset{n\to \infty }{\mathrm{limsup}}{\tilde{\lambda }}_{h_n,1}$$

$$\le \underset{n\to \infty }{\mathrm{limsup}}{\tilde{r}}_{h_n}(w)\text{mbox}(Ritz-valuesminimizetheRayleigh-quotient)$$

$$\underset{n\to \infty }{\mathrm{limsup}}{\tilde{\lambda }}_{h_n,1}$$

$$=\underset{n\to \infty }{\mathrm{limsup}}{\tilde{a}}_{h_n}(w,w)\text{mbox}(definitionofRayleighQuotient)$$

$$=\underset{n\to \infty }{\mathrm{limsup}}(a(w,w)+{\tilde{a}}_{h_n}(w,w)-a(w,w))$$

$$\le \underset{n\to \infty }{\mathrm{limsup}}\left({\lambda }_1+ϵ+a^1(I_{h}Tw,I_{h}Tw)-a^1(Tw,Tw)\right)$$

$$\underset{h\to 0}{\mathrm{liminf}}{\tilde{\lambda }}_{h,1}\ge {\lambda }_1.$$

$${\Vert {\tilde{v}}_{h_n}\Vert }_{H^1}\le C$$

$$a^1(I_{h_n}T{\tilde{v}}_{h_n},I_{h_n}T{\tilde{v}}_{h_n})=a^1(I_{h_n}T{\tilde{v}}_{h_n}-J_n,I_{h_n}T{\tilde{v}}_{h_n}-J_n)$$

$${\Vert I_{h_n}T{\tilde{v}}_{h_n}-T\tilde{v}\Vert }_{L^2}$$

$$\le C_1{\Vert {\tilde{v}}_{h_n}-{\tilde{v}}_{h_m}\Vert }_{L^2}+{\Vert I_{h_n}T{\tilde{v}}_{h_m}-T{\tilde{v}}_{h_m}\Vert }_{L^2}+{\Vert T{\tilde{v}}_{h_m}-T\tilde{v}\Vert }_{L^2}$$

$$\le 2\epsilon /3+{\Vert I_{h_n}T{\tilde{v}}_{h_m}-T{\tilde{v}}_{h_m}\Vert }_{L^2}.$$

$${\Vert I_{h_n}T{\tilde{v}}_{h_m}-T{\tilde{v}}_{h_m}\Vert }_{L^2}\le \epsilon /3.$$

$${\Vert I_{h_n}T{\tilde{v}}_{h_n}-T\tilde{v}\Vert }_{L^2}\le ϵ.$$

$$a^0(\tilde{v},\tilde{v})$$

Full text

1 Introduction

The dynamic Laplacian is a second-order partial differential operator underlying a range of methods for computing finite-time coherent sets in finite-time non-autonomous dynamical systems. It was introduced in [12] in the context of defining sets that remain coherent in a Lagrangian sense via dynamic isoperimetry. Coherent sets are time-dependent families of sets whose boundaries remain small relative to the volume of the set as the family evolves according to the nonlinear dynamics; extensions to weighted, curved manifolds and non-volume-preserving dynamics were made in [15]. Coherent sets are captured by the eigenvectors of the dynamic Laplacian corresponding to the leading eigenvalues (i.e. those closest to [math]). In particular, level sets of the eigenvector corresponding to the first nontrivial eigenvalue can be used to partition the domain into two coherent sets. A dynamic Cheeger inequality [12, 15] links this eigenvalue to the ratio of boundary size to volume. Moreover, if $n$ eigenvalues are close to zero followed by a spectral gap, this forces the eigenvectors to be close to linear combinations of indicator functions on an $n$-partition [6]. We note that the dynamic Laplacian is the Laplace-Beltrami operator of a weighted manifold [19] and therefore can be used as the time-independent generator of a diffusion process approximating the given (time-dependent) advection-diffusion process. Eigenvectors corresponding to small eigenvalues decay slowest under this diffusion process and its almost-invariant sets in the sense of [7]. Algorithmically, coherent sets can be extracted from the eigenvectors of the dynamic Laplacian via, e.g., clustering techniques [11], optimising eigenbasis separation [8], optimising sublevel sets [17], or sparse eigenbasis approximation [18].

In most cases, the dynamic Laplacian eigenproblem must be solved numerically. To this end, a scheme based on radial basis functions had been proposed [13], which showed high order of convergence, but suffered from a number of drawbacks like high sensitivity with respect to the radius parameter, a non-real spectrum and non-sparseness of the discretized operator. In [14], two finite element schemes were proposed (the “Cauchy-Green” (CG) and the “Transfer Operator” (TO) approach), which eliminated each of these drawbacks.

Experimentally, only piecewise linear elements were considered in [14]. In this paper, we consider higher-order (and in particular quadratic) elements and analyse convergence properties both theoretically and experimentally. We provide an efficient implementation in the Julia package CoherentStructures.jl. We find that using $P^{2}$ elements can give a higher asymptotic order of convergence compared to $P^{1}$ elements in the “CG” approach. For the “CG” approach, classical theory concerning eigenproblems in FEM applies. We also provide some test cases where using $P^{2}$ elements can greatly reduce the amount of information needed to calculate partitions of the domain that show important dynamical features. The question of convergence in the “TO” approaches is more subtle. We prove convergence of eigenvalues and eigenvectors for a family of TO approaches with $P^{1}$ elements, but the proof does not give any insight into the convergence rates that should be expected. For $P^{2}$ elements we do not observe asymptotically higher orders of convergence even for the simple example of a one-dimensional shift-map on the torus. This suggests that using $P^{2}$ elements does not have substantial benefits when using the “TO” approach and that one should use the simpler and well-performing linear $P^{1}$ elements in the “TO” schemes.

2 The Dynamic Laplacian

Let $\mathcal{I}\subset\mathbb{R}$ denote a finite subset of time and consider a finite family $(\Omega_{t})_{t\in\mathcal{I}}$ of open bounded subsets of $\mathbb{R}^{d}$ with Lipschitz boundary. For each $t\in\mathcal{I}$ let $T_{t}:\Omega_{0}\rightarrow\Omega_{t}$ be a volume-preserving diffeomorphism. We assume that $T_{t}$ is sufficiently regular so that $T_{t}$ and $T_{t}^{-1}$ can be smoothly extended to the boundary, and that $0\in\mathcal{I}$ with $T_{0}$ being the identity. A typical setting in which these conditions apply are $T_{t}$ taken to be time-$t$ flow maps of a divergence-free vector-field.

Denote the Laplace operator on $\Omega_{t}$ by $\Delta_{t}$ for each $t\in\mathcal{I}$. Then the dynamic Laplacian (an operator on $\Omega_{0}$) is given by

[TABLE]

where $T_{t}^{*}:L^{2}(\Omega_{t})\to L^{2}(\Omega_{0})$ denotes the pullback by $T_{t}$ defined by $T_{t}^{*}f=f\circ T_{t}$, and $T_{t,*}:L^{2}(\Omega_{0})\to L^{2}(\Omega_{t})$ is the pushforward defined by $T_{t,*}f=f\circ T_{t}^{-1}$.

Standard PDE-theoretic arguments can be used to show that $\Delta^{dyn}$ is a uniformly elliptic second-order partial-differential operator [10, 12], with weak form [14]

[TABLE]

where $\ell^{d}$ is the $d$-dimensional Lebesgue measure. The time set $\mathcal{I}$ can also be a compact interval, and the dynamic Laplacian can be defined by means of an integral over $\mathcal{I}$ [12, 15]

[TABLE]

where $|\mathcal{I}|$ now denotes the length of the interval $\mathcal{I}$, we do not consider this generalization further here.

The uniform ellipticity of the dynamic Laplacian ensures that the bilinear form in (1) is coercive under suitable boundary conditions. These are determined by the choice of the underlying space $S$ on which (1) acts. For natural (Neumann) boundary conditions, we look at $a$ on $\hat{H^{1}}\times\hat{H^{1}}$, where $\hat{H}^{1}=\hat{H}^{1}(\Omega_{0}):=\{v\in H^{1}(\Omega_{0}):\int f\;d\ell^{d}=0\}$ is the space of mean-free $H^{1}$ functions. For homogeneous Dirichlet boundary conditions, $a$ must be taken to act on $H^{1}_{0}\times H^{1}_{0}$. With either choice of boundary conditions, the variational form of the of the eigenproblem $\Delta^{dyn}u=\lambda u$ becomes

[TABLE]

In both cases, ellipticity ensures that there exists a countable sequence of pairwise orthogonal eigenvectors $u_{0},u_{1},\ldots$ corresponding to real eigenvalues $0\geq\lambda_{0}\geq\lambda_{1}\geq\cdots$. Furthermore, the span of the eigenvectors is dense in $L^{2}$ (Dirichlet) and mean free $L^{2}$ functions (Neumann) [10, 12].

2.1 Discretisation with finite elements

A natural discretisation of the eigenproblem (2) is by using a finite element method (FEM) [14]. In the standard Galerkin discretization of (2), $S^{\prime}\subset S$ is taken to be a finite dimensional approximation space spanned by some basis $(\varphi_{1},\dots,\varphi_{n})$. We now find approximate solutions $u\in S^{\prime}$, $\lambda\in\mathbb{R}$, that satisfy (2) with $S^{\prime}$ taking the place of $S$. In matrix form, the coefficients $\textbf{u}=(\textbf{u}_{1},\ldots,\textbf{u}_{n})$ of $u$ with respect to the basis $(\varphi_{1},\dots,\varphi_{n})$ are found by solving the generalized eigenvalue problem

[TABLE]

Here, $D=(D_{i,j})=a(\varphi_{i},\varphi_{j})$ is referred to as the stiffness matrix and $M=(M_{i,j})=\langle\varphi_{i},\varphi_{j}\rangle_{L^{2}}$ as the mass matrix. We use $P^{k}$ Lagrange nodal basis functions on some (triangular or quadrilateral) mesh for $\varphi_{1},\dots,\varphi_{n}$, but other choices are of also possible in general (see [9]).

Using a finite element method for approximating $\Delta^{dyn}$ has a number of advantages: The matrix formulation (3) inherits self-adjointness from the continuous problem and thus the computed eigenpairs are always real. With a suitably localized basis, $D$ and $M$ are both sparse. Sparse Hermitian generalized eigenproblems are well known in the literature, and a number of algorithms exist for efficiently solving them [1]. Also, finite element methods have been widely studied, and there are a range of theoretical results regarding convergence (cf [9]) that are applicable to some solution approaches.

While the entries of the mass matrix can be computed exactly, the integrand in (1) in the entries of the stiffness matrix can be of extremely high variation locally, making an accurate computation possibly expensive. Several approaches for approximating $D$ were suggested in [14], which we investigate further in the following sections.

3 Convergence for the CG approach

In the “CG approach” of [14], a quadrature formula is used to approximate (1). This is the standard way of solving elliptic eigenproblems with finite elements.

We work on $P^{k}$-Lagrange finite element spaces for a family of quasi-uniform meshes 111A family $\{\mathcal{T}_{h}\}_{h>0}$ of meshes is quasi-uniform, if it is shape-regular and there exists a constant $c>0$ so that for any $h>0$ and any element $K\in\mathcal{T}_{h}$ one has that $h_{K}\geq ch$, where $h_{K}$ is the diameter of the element $K$. Shape regularity means that there is a constant $\sigma_{0}$ so that uniformly for all $K$ it holds that $\sigma=\frac{h_{K}}{\rho_{K}}\leq\sigma_{0}$. Here $\rho_{K}$ is the diameter of the largest ball that can be inscribed in $K$[9].

triangulations (in $2$ or $3$ dimensions) $\mathcal{T}_{h}^{0}$ of $\Omega_{0}$, where the family parameter $h$ is related to the mesh size. Let $\varphi^{1}_{h},\ldots,\varphi^{N(h)}_{h}$ be the associated nodal basis, $S_{h}^{0}:=\mbox{span}\{\varphi^{1}_{h},\ldots,\varphi^{N(h)}_{h}\}$ the finite element approximation space, and $\hat{S}_{h}^{0}:=\hat{H}^{1}\cap S_{h}^{0}$. Moreover, let $\lambda_{i,h}$ be the Ritz-values with corresponding Ritz-vectors $v_{i,h}$ in $\hat{S}_{h}^{0}$, i.e. $\lambda_{i,h}$ is the minimum value of $a(v,v)$ over the set $\{v\in\hat{S}_{h}^{0};{\|v\|}_{L^{2}}=1\langle v,v_{j}\rangle=0\mbox{ for }j=1,\cdots,(i-1)\}$ and $v_{i,h}$ is the corresponding minimizer.

Under the assumption that the eigenvectors of $a:\hat{H}^{1}\times\hat{H}^{1}\to\mathbb{R}$ are in $H^{k+1}$, classical results from FEM theory [24, 9] give that (where $C$ stands for a generic constant):

[TABLE]

Moreoever, it holds for simple (and $L^{2}$-normalised) eigenvectors that

[TABLE]

In the case that the problem (1) satisfies the technical assumption of being regularizing [9] (i.e ${\|Pu\|}_{H^{2}}\leq C{\|u\|}_{L^{2}}$), the “Aubin-Nitsche trick” gives

[TABLE]

Note that not all elliptic problems are regularizing. We thus cannot expect (6) to hold in general.

If the entries of the stiffness and the mass matrix are approximated with a quadrature rule of order $2k-1$, the convergence orders are unaffected, provided that the quadrature points contain a $\mathcal{P}^{k}$ unisolvent set [2]. Moreoever, the results generalize to non-simple eigenspaces and Kato’s subspace distance used for the error[2].

3.1 Numerical Experiments

We now aim at reproducing the predicted convergence rates in numerical experiments. As a reference, we compute ($L^{2}$ normalized) eigenpairs $(\lambda_{i},v_{i})$ on a very fine mesh and estimate the error in the eigenvector $v_{i,h}$ by computing the $L^{2}$-distance of $v_{i,h}$ to the closest reference eigenvector, given by the expression

[TABLE]

The $L^{2}$ inner-product is approximated by first interpolating to the fine grid, and then calculating the inner product there using quadrature.

To generalize this to higher-dimensional eigenspaces, let $\tilde{V}$ and $V$ be two $m$-dimensional subspaces of $L^{2}$ with orthonormal bases given by $\{\tilde{v}_{1},\dots,\tilde{v}_{m}\}$ and $\{v_{1},\dots,v_{m}\}$. Let $\tilde{P},P$ be orthogonal projections onto $\tilde{V}$ and $V$ respectively. As a measure of the subspace error [20], we maximize ${\|\tilde{v}-P\tilde{v}\|}_{L^{2}}$ over $\tilde{v}\in\tilde{V}$ with ${\|\tilde{v}\|}_{L^{2}}=1$. This is equivalent to maximizing $\sqrt{1-{\|P\tilde{v}\|}^{2}}$, which has maximum given by $\sqrt{1-\rho_{min}(P)^{2}}$ where $\rho_{min}$ is the smallest singular value of $P$ on $\tilde{V}$ which is also that of the matrix $(\langle\tilde{v}_{i},v_{j}\rangle_{L^{2}})_{i,j=1}^{m}$.

In general, one is not directly interested in the eigenvectors themselves, but in a partition of the domain obtained by a suitable post-processing of the eigenvectors (e.g., as mentioned, thresholding, clustering or sparse eigenbases). This motivates us making a qualitative comparison of such partitions based on eigenvectors for different approximation spaces. Here, we focus in particular on the question of how well Lagrangian coherent sets can be found with as little information about the flow as possible. Evidently, in this low data case some features may not be accurately resolved even if the eigenvectors clearly give some indication of their existence. We therefore also look at how coarse the grid can be made without affecting the topology of such a partition in some test cases.

In the sequel, mesh sizes of the form $n\times m$ refer to a regular triangular mesh with $n\times m$ (non-inner) nodes in each direction.

For all time-integration done in this work we used a relative and absolute tolerance of $10^{-8}$ and the DifferentialEquations.jl Julia package [22] with the BS5() solver [3]. Derivatives of flow maps were approximated with second order central finite differences. The stiffness matrix was calculated by approximating the integral in (1) on each element with nodal basis functions for $u$ and $v$ using quadrature. The results from this were then summed over all elements in the support of a basis function to give the corresponding entry in the stiffness matrix. Similarly, the integral in the mass matrix was also calculated element-wise and then additively combined to give the full mass matrix. The CoherentStructures.jl package internally uses the JuAFEM.jl package [5].

3.1.1 Standard Map

As a simple first example, we consider two iterations of the standard map

[TABLE]

on the $2$-torus $\mathbb{S}^{1}\times\mathbb{S}^{1}$ with parameter $a=0.971635$. This is the first example considered in [13] and is a weakly nonlinear map. Figure 1 shows the eigenvalue and eigenspace errors in dependence of the mesh width for triangular $P^{1}$ and $P^{2}$ Lagrange elements.

The reference solution was computed with $P^{2}$ elements on a regular $1025\times 1025$ grid with quadrature order $5$. In this case, the experimentally observed rates almost perfectly agree with the predictions (4) and (6) even for quadrature order $2$ in the computation of the entries of stiffness- and mass-matrix. Improvements in errors ranged between one and two orders of magnitude when moving from $P^{1}$ to $P^{2}$ elements. To be able to directly compare these results with the corresponding plots for TO methods in Figure 5, we used a two-dimensional eigenspace in Figure 1.

3.1.2 Cylinder Flow

As a second example, we consider the cylinder flow map used in [13] based on [16]. This is a time-dependent flow on the cylinder222More precisely, we looked at the flow on $\mathbb{S}^{1}\times[10^{-2},\pi-10^{-2}]$ as our numerical approximation of the flow-map was not well-behaved at the boundary $\mathbb{S}^{1}\times[0,\pi]$ defined by the non-autonomous ordinary differential equation

[TABLE]

where $A(t)=1+0.125\sin(2\sqrt{5}t)$, $G(\psi)=1/(\psi^{2}+1)^{2}$ and $g(x,y,t)=\sin(x-\nu t)\sin(y)+y/2-\pi/4$. Here the parameters $c=0.5,\nu=0.5,\varepsilon=0.25$ were used, the time-interval was $[0,40]$.

The reference solution was computed on a regular mesh on $1025\times 1025$ nodes with triangular $P^{2}$-Lagrange elements and quadrature order $8$. We used quadrature order $5$ for the numerical experiments (since we did not observe the same rates for smaller orders).

The orders of convergence observed for the eigenvalue and eigenvector errors is surprisingly low; the experimental values are almost exactly half of the orders predicted in (4), (5) or (6). The observed order of the eigenvector errors for $P^{1}$ elements is particularly poor. We also note that the slopes shown in the figure do not remain consistent when varying the quadrature order.

For this flow (and other highly stretching and highly nonlinear systems), the mean diffusion tensor $A=[DT_{t}]^{-1}([DT_{t}]^{-1})^{T}$ used in the weak formulation of the dynamic Laplacian (1) has extremely high variation locally (cf. Figure 3). We suspect that even the reference grid of $1025\times 1025$ nodes is not sufficiently resolved to be in a regime where the convergence orders predicted by the theory can be observed.

3.2 Bickley jet

The Bickley jet flow is a well-known test case first introduced in [23]. The flow is defined by a stream-function

[TABLE]

with constants $U_{0},L_{0},A_{i},k_{i},c_{i}$ on a cylindrical domain (see [23]). We considered the flow for a timespan of $40$ days. $8$-partitions were computed by k-means clustering on $200\times 60$ values of the leading eigenvectors of the FEM approximation. The mesh widths in Figure 4 below are the lowest (at the aspect ratio $10:3$) for which the topology of the clustering result does not change. It is possible to obtain the same result (topologically) with $P^{2}$ elements at significantly reduced cost (in terms of the number of quadrature point used) to $P^{1}$ elements. In many cases, the computation requiring the most runtime is the computation of diffusion tensors by time-integration at the quadrature points. As shown in Figure 4, in this case this is reduced by a factor of 25 by using $P^{2}$ elements.

4 Convergence for the TO approach

The CG approach has the disadvantage of requiring the numerical approximation of the derivative $DT_{t}$ and a subsequent quadrature which might be challenging if $DT_{t}$ has high variation locally.

An alternative approach [14] is to rewrite the weak form (1) of the dynamic Laplacian operator using the transfer operator $T_{t}$. In the sequel we discuss the case $t\in\mathcal{I}:=\{0,1\}$, although the results hold for general finite $\mathcal{I}$. Letting $a^{t}$ be the weak-form of the Laplacian on $\Omega_{t}$, we have

[TABLE]

where we write $T:=T_{1}$, and for brevity overload this notation so that $T$ acts on elements of $\Omega_{0}$ in the usual way and acts on functions in $L^{2}(\Omega_{0})$ as $Tu:=T_{1,*}u$. In the transfer operator (TO) approach, we approximate $T$ on a suitable subspace of $H^{1}$ by an operator of the form $I_{h}T$, where $I_{h}$ is some projection operator. Let $\mathcal{T}^{1}_{h}$ be a quasi-uniform family of triangulations on $\Omega_{1}$, $S^{1}_{h}$ the approximation space with Lagrange $P^{k^{\prime}}$ finite elements and $\hat{S}_{h}^{1}:=S_{h}^{1}\cap\hat{H}^{1}(\Omega_{1})$. Note that in general, the collection $\mathcal{T}_{h}^{0}$ (introduced at the beginning of Section 3) and $\mathcal{T}^{1}_{h}$ are unrelated. Additionally, we assume that the meshes used are nested in the sense that $S_{h^{\prime}}\subset S_{h}$ for $h<h^{\prime}$.

We consider the following options for $I_{h}$ in the sequel:

the $L^{2}$-orthogonal projection onto $S_{h}^{1}$ ($L^{2}$-Galerkin), 2. 2.

the $H^{1}$-orthogonal projection onto $S_{h}^{1}$ ($H^{1}$-Galerkin), 3. 3.

the (Lagrange $P^{k}$) nodal interpolation on $S_{h}^{1}$ (collocation).

Case 3. corresponds to the “TO approach” from [14], with the difference between the adaptive and non-adaptive approaches being the choice of the mesh $\mathcal{T}^{1}_{h}$.

We now define the ($h$-dependent) symmetric bilinear form

[TABLE]

As done in [14], eigenpairs are numerically approximated by calculating generalized eigenpairs of $(\tilde{D}_{h},M_{h})$ given by

[TABLE]

As $a^{0}$ is elliptic on $\hat{H}^{1}(\Omega_{0})$, it follows immediately from (9) that $\tilde{a}_{h}$ is elliptic uniformly in $h$ on $\hat{S}_{h}^{0}$, with lower bound of the ellipticity constant being half of that of $a^{0}$.

Properties of $I_{h}$.

In order to investigate the convergence of the numerically calculated eigenpairs to those of the dynamic Laplacian, we will need to consider the following conditions for $I_{h}T$ acting on $S_{h}^{0}$:

1. (i)

exactness for constant functions: $I_{h}T\chi_{\Omega_{0}}=\chi_{\Omega_{1}}$, 2. (ii)

$H^{1}$-convergence: $\lim\limits_{h\rightarrow 0}{\|I_{h}Tf-Tf\|}_{H^{1}}=0$ for any $f\in\cup_{h}S_{h}^{0}$, 3. (iii)

$L^{2}$-stability: $\sup_{h<h^{\prime}}{\|I_{h}T|S_{h^{\prime}}^{0}\|}_{L^{2}\rightarrow L^{2}}\leq C$.

We now consider whether these three conditions are satisfied by the three variants of $I_{h}$:

1. (i)

It is easy to see that condition (i) is satisfied for all three variants, by the fact that $T:=T_{1}$ is volume-preserving and $I_{h}$ is exact on constant functions. 2. (ii)

Condition (ii) is satisfied by the $H^{1}$-Galerkin approach. Under the technical assumption that $\mathcal{T}^{1}_{h}$ is quasi-uniform [9], we also get (ii) for the $L^{2}$-Galerkin approach. 3. (iii)

Condition (iii) is trivially satisfied by the $L^{2}$-Galerkin approach.

We do not consider condition (iii) for the $H^{1}$-Galerkin approach. It remains to show (ii) and (iii) for the case of $I_{h}$ being the nodal interpolation operator. In the appendix, we give proofs of the following two Lemmas for $P^{1}$-triangular elements under the assumption that both meshes are quasi-uniform. For $h^{\prime}>h>0$, let $\mathcal{T}^{0}_{h^{\prime}}$ and $\mathcal{T}^{1}_{h}$ be families of quasi-uniform triangulations of $\Omega_{0}$ and $\Omega_{1}$ respectively. Let $I_{h}:C(\overline{\Omega_{1}})\rightarrow S_{h}^{1}$ be the $P^{1}$-Lagrange nodal interpolation operator corresponding to $\mathcal{T}^{1}_{h}$.

Lemma 1** ($L^{2}$-stability of $P^{1}$-nodal interpolation).**

There exists a constant $C>0$ so that ${\|I_{h}Tv\|}_{L^{2}}\leq C{\|v\|}_{L^{2}}$ for all $v\in S_{h^{\prime}}^{0}$ and all $h<h^{\prime}$.

Lemma 2** ($H^{1}$-convergence of $P^{1}$-nodal interpolation).**

If $v$ is piecewise $C^{\infty}$ on $\Omega_{1}$ then

[TABLE]

Convergence for homogeneous Neumann Boundary.

Let $r(v)=a(v,v)/\langle v,v\rangle_{L^{2}}$ and $\tilde{r}_{h}(v)=\tilde{a}_{h}(v,v)/\langle v,v\rangle_{L^{2}}$ be the Rayleigh quotients for $a$ and $\tilde{a}_{h}$ respectively. Let $\lambda_{1}$ and $\tilde{\lambda}_{h,1}$ be the corresponding smallest Ritz values on $\hat{H}^{1}$ and $\hat{S}_{h}^{0}$, respectively and $v$, $\tilde{v}_{h}$ the corresponding Ritz vectors, normalized to have unit $L^{2}$-norm.

From the aforementioned properties of $I_{h}T$ and the fact that the Ritz vectors minimize the Rayleigh quotient we now derive a convergence result for $\tilde{v}_{h}$ to $v$ in two parts:

Theorem 3**.**

Assume that $I_{h}$ satisfies the $H^{1}$-convergence property (ii), then

[TABLE]

Proof.

Pick any sequence $h_{n}$ with $h_{n}\rightarrow 0$. Then for arbitrary $w\in\hat{S}_{h}^{0}$

[TABLE]

Now for arbitrary $\epsilon>0$, pick333That we can do this follows from the classical theory of convergence of eigenvalues of the Galerkin approximation such as equation (4) in the case that the eigenvectors are sufficiently smooth at the boundary. some $w\in\hat{S}_{h^{\prime}}^{0}$ where $h^{\prime}=h^{\prime}(\epsilon)$ with ${\|w\|}_{L^{2}}=1$ so that $r(w)<\lambda_{1}+\epsilon$. Then by the assumption that our meshes are nested, $w$ is in the domain of $\tilde{a}_{h_{n}}$ for sufficiently large $n$ and we have:

[TABLE]

For this fixed $v$ we know that $a^{1}(I_{h}Tw,I_{h}Tw)-a^{1}(Tw,Tw)\rightarrow 0$ by the $H^{1}$ convergence assumption, and by $H^{1}$-continuity of $a^{1}$. This gives the result. ∎

Theorem 4**.**

Assume that in addition to exactness on constant functions (i), $I_{h}$ satisfies both the $H^{1}$-convergence condition (ii) and the $L^{2}$-stability condition (iii), then

[TABLE]

Proof.

Assume that $I_{h}$ satisfies both properties, but (for the sake of contradiction) that there exists a sequence $h_{n}\rightarrow 0$ monotonically with $\tilde{\lambda}_{h_{n},1}\rightarrow\beta<\lambda_{1}$. Note that by positivity of $\tilde{a}_{h}$ the sequence cannot diverge to $-\infty$, and by Theorem 3 it cannot diverge to $+\infty$. Hence, using ellipticity of $a^{0}$ on $\hat{H}^{1}$ we get that $|\tilde{v}_{h_{n}}|_{H^{1}}\leq C$ for some $C>0$, otherwise the $a^{0}$ term would go to infinity. Here and below, $C$ represents a generic constant independent of $h_{n}$. By the Poincaré-Friedrichs theorem, it follows that

[TABLE]

Recall that the Banach-Alaoglu theorem states that the closed unit ball in $H^{1}$ is weakly sequentially compact, while the Rellich-Kondrachev theorem asserts that norm-bounded $H^{1}$ sets are $L^{2}$-precompact. Hence we can assume that (going to a subsequence) $\tilde{v}_{h_{n}}\rightharpoonup\tilde{v}$ in $H^{1}$ and $\tilde{v}_{h_{n}}\rightarrow\tilde{v}$ in $L^{2}$ for some $\tilde{v}\in H^{1}$. Note that the two limits are indeed the same since $\cdot\mapsto\langle\cdot,v\rangle_{L^{2}}$ is a continuous linear functional on $H^{1}$ for all $v\in L^{2}$.

As $\tilde{v}_{h_{n}}\in\hat{H}^{1}$, it follows that $\tilde{v}\in\hat{H}^{1}$, the $L^{2}$-convergence gives us further that ${\|v\|}_{L^{2}}=1$.

Without loss of generality, assume that $\ell^{d}(\Omega_{1})=1$ and define $J_{n}:\Omega_{0}\to\mathbb{R}$ as the constant function taking the value $\int_{\Omega_{1}}I_{h_{n}}T\tilde{v}_{h_{n}}d\ell^{d}$, and observe that

[TABLE]

We know from the positivity of $a^{0}$ that $a^{1}(TI_{h_{n}}\tilde{v}_{h_{n}},TI_{h_{n}}\tilde{v}_{h_{n}})\leq\tilde{r}_{h_{n}}(\tilde{v}_{h_{n}})$ and as the right hand side converges, the left hand side is bounded by a constant that does not depend on $n$. Hence, by ellipticity of $a^{1}$ on mean-free $H^{1}$ functions, it follows from (11) and the Poincaré-Friedrichs inequality that ${\|I_{h_{n}}T\tilde{v}_{h_{n}}-J_{n}\|}_{H^{1}}<C$, yielding that $I_{h_{n}}T\tilde{v}_{h_{n}}$ is $H^{1}$-bounded as long as $J_{n}$ is – which is the case by $L^{2}$-stability, the $L^{2}$-$L^{2}$ continuity of $T$, and the fact that ${\|\tilde{v}_{h_{n}}\|}_{L^{2}}=1$. Hence $I_{h_{n}}T\tilde{v}_{h_{n}}$ is (uniformly in $n$) $H^{1}$-bounded. Going to a subsequence, it is therefore weakly $H^{1}$- and strongly $L^{2}$-convergent (as above). As the weak-$H^{1}$ and $L^{2}$-limits coincide, to show that the sequence converges to $T\tilde{v}$ weakly in $H^{1}$, it is enough to show $L^{2}$-convergence.

Let $C_{1}$ be the stability constant from the $L^{2}$-stability condition. We know that as $v_{h_{n}}\rightarrow\tilde{v}$ in $L^{2}$, it holds that $T\tilde{v}_{h_{n}}\rightarrow T\tilde{v}$. Therefore given $\varepsilon>0$, pick $m$ so that ${\|T\tilde{v}_{h_{m}}-T\tilde{v}\|}_{L^{2}}\leq\varepsilon/3$ and and $C_{1}{\|\tilde{v}_{h_{n}}-\tilde{v}_{h_{m}}\|}_{L^{2}}\leq\varepsilon/3$ for all $n\geq m$. Also, by the nesting property of the meshes $\tilde{v}_{h_{n}}-\tilde{v}_{h_{m}}\in\hat{S}_{h_{n}}^{0}$ and hence

[TABLE]

Now using the $H^{1}$-convergence property, pick $n_{0}\geq m$ so that for $n\geq n_{0}$, we have

[TABLE]

It follows that for all $n\geq n_{0}$

[TABLE]

Summarizing, we now have proved that

$\tilde{v}_{h_{n}}\rightharpoonup\tilde{v}$ in $H^{1}$, 2. 2.

$I_{h_{n}}T\tilde{v}_{h_{n}}\rightharpoonup T\tilde{v}$ in $H^{1}$, and 3. 3.

$\tilde{v}\in\hat{H}^{1}$ with ${\|\tilde{v}\|}_{L^{2}}=1$

The functions $u\mapsto\sqrt{a^{0}(u,u)}$ and $u\mapsto\sqrt{a^{1}(u,u)}$ are norms that are equivalent to the $\hat{H}^{1}$-norm by the Poincaré-Friedrichs inequality444The weak topologies also coincide, as for any $f\in\hat{H}^{1}$, the Riesz representation theorem gives a $\tilde{f}\in\hat{H}^{1}$ so that $\langle\cdot,f\rangle_{H^{1}}=a_{0}(\cdot,\tilde{f})$. By the Banach-Steinhaus theorem, norms are weakly sequentially lower-semicontinuous. The $L^{2}$ convergence of $I_{h_{n}}Tv_{h_{n}}$ also gives that $J_{n}\rightarrow 0$ in $L^{2}$, and therefore also in $H^{1}$ as $J_{n}$ is constant). Thus we have:

[TABLE]

yielding

[TABLE]

This is a contradiction to the defintion of $\lambda_{1}$, which requires that $a(\tilde{v},\tilde{v})\geq\lambda_{1}$ and thus the claim is proved. ∎

Theorem 5**.**

The conclusion of the previous theorem also holds for $I_{h}$ coming from the $H^{1}$-Galerkin approach.

Proof.

The previous proof required the $L^{2}$-stability only in two places. The first was in bounding $J_{n}$, this is trivially bounded as ${\|I_{h}\|}_{H^{1}\rightarrow H^{1}}=1$ and $\tilde{v}_{h_{n}}$ is $H^{1}$-bounded. The second was in showing that $I_{h_{n}}T\tilde{v}_{h_{n}}\rightharpoonup T\tilde{v}$. But this follows from the fact that for $f\in H^{1}$:

[TABLE]

As $I_{h}f\rightarrow f$ in $H^{1}$ and $T\tilde{v}_{h_{n}}\rightharpoonup T\tilde{v}$, this gives the claim. ∎

Theorem 6** (Eigenvector convergence).**

Assume that the eigenspace corresponding to $\lambda_{1}$ is one-dimensional. Then $|\langle\tilde{v}_{h},v\rangle_{L^{2}}|\rightarrow 1$ for $h\rightarrow 0$.

Proof.

The proof of Theorem 4 shows that for any monotone sequence $h_{n}\rightarrow 0$ for which $\lambda_{h_{n}}$ converges to $\lambda_{1}$, every in $L^{2}$ convergent subsequence of $\tilde{v}_{h_{n}}$ converges to some $\tilde{v}\in H^{1}$ with $r(\tilde{v})=\lambda_{1}$. Theorems 4 and 3 show that $\lambda_{h_{n}}\rightarrow\lambda_{1}$ for any $h_{n}\rightarrow 0$. As $\tilde{v}_{h_{n}}$ is bounded in $H^{1}$ and therefore pre-compact in $L^{2}$, this means that $\tilde{v}_{h_{n}}\rightarrow\tilde{v}$ in $L^{2}$ with $r(\tilde{v})=\lambda_{1}$, from which it follows that $\tilde{v}=\pm v$, giving the result for a subsequence. We get the final result by noting that if $\tilde{v}_{h_{n}}$ has a subsequence $\tilde{v}_{h_{n_{k}}}$ for which $|\langle\tilde{v}_{h_{n_{k}}},v\rangle|\rightarrow c\neq 1$, by the above it must have a (further) subsequence that converges to a different value, a contradiction. ∎

Theorem 7** (Higher eigenpairs).**

The results of the previous theorems generalize to the $m$-th eigenpair for $m>1$.

Proof.

The generalization of Theorem 3 can be proven by the same method as the case $m=1$: instead of chosing a $w$ with $r(w)<\lambda_{1}+\epsilon$, choose a series of pairwise-orthogonal $(w_{i})_{i=1\dots m}\in\hat{S}^{0}_{h}$ with the property that $r(w_{i})\leq\lambda_{i}+\epsilon$.

The proof of Theorem 4 must only be modified to show that the limit point $\tilde{v}_{m}$ is $L^{2}$-orthogonal to the previous eigenvectors $\operatorname{span}\{v_{1},\dots,v_{m-1}\}$. Writing $\tilde{v}_{m,h}$ for the $m$-th Ritz-vector of $\tilde{a}_{h}$ on $\hat{S}_{h}^{1}$, this readily follows (by induction) from the fact that $\tilde{v}_{m,h_{n}}\perp\operatorname{span}\{\tilde{v}_{1,h_{n}},\dots,\tilde{v}_{m-1,h_{n}}\}$ and thus also $\langle\tilde{v}_{i,h_{n}},Tv_{i}\rangle\rightarrow 0$ in $L^{2}$ for $i=1,\dots,m-1$.

The proof of the generalization of Theorem 6 then follows from the generalization of Theorem 4. ∎

Convergence for homogeneous Dirichlet boundary.

For the case of a homogeneous Dirichlet boundary, we must replace $\hat{H}^{1}$ with $H^{1}_{0}$. Here we require the additional condition:

1. (i’)

Boundary preservation: for $v\in H^{1}_{0}(\Omega_{1})\cap S_{h}^{0}$ it holds that $I_{h}Tv\in H^{1}_{0}(\Omega_{1})$

We note that the $I_{h}$ from collocation and $H^{1}$-Galerkin have this property.

Theorem 8**.**

The same convergence results hold also for the Ritz-values for homogeneous Dirichlet boundary conditions, provided that $I_{h}$ satisfies (i’) in addition to satisfying (ii) and (iii) from the previous section.

Proof.

Same as the proofs for the Neumann case, except that coercivity of the bilinear forms involved on $\hat{H}^{1}$ is replaced with that on $H^{1}_{0}$. The boundary preservation property ensures we do not need $J_{n}$ in the proof of Theorem 4, the proof otherwise goes on exactly the same lines. ∎

4.1 Numerical Experiments

The three TO approaches yield a stiffness-matrix $D$ of the form

[TABLE]

where $A$ is the representation matrix of $I_{h}T$ and $D_{0},D_{1}$ result from the bilinear form of the static Laplacian. In the case of the $L^{2}$- or $H^{1}$-Galerkin approach, the matrix $A$ has the form

[TABLE]

where $\tilde{A}=\langle\varphi_{i},T\varphi_{j}\rangle$ and $G_{i,j}=\langle\varphi_{i},\varphi_{j}\rangle$. Note that sparsity of $G$ and $\tilde{A}$ does not necessarily imply sparsity of $A$. A naive calculation of the full matrix $A$ therefore renders the computation of (some) eigenvalues of $D$ too expensive for larger problems. We therefore do not include the Galerkin-TO approaches in the numerical experiments below, further work is needed to determine if the Galerkin-TO methods can be modified to overcome this issue.

All TO results shown below are therefore the collocation-based ones. For the “non-adaptive TO” experiments, we used identical regular meshes for the initial and final triangulations $\mathcal{T}_{h}^{0}$ and $\mathcal{T}_{h}^{1}$. For the “adaptive TO” experiments, the initial mesh was regular but the final mesh was a Delaunay-triangulation (with the VoronoiDelaunay.jl package [21]) of the images of the initial mesh points under the flow map $T$. The theory outlined above requires (i) nested meshes, and (ii) uniform bounds on shape-regularity and quasi-uniformity constants of the meshes as they are refined; each of these properties is difficult to guarantee in general for the image triangulations. We nevertheless include the adaptive TO in the plots for comparison.

We note that the finest grid used for the TO experiments is $257\times 257$. This is not as fine as the finest grid used in the CG approach ($513\times 513)$. The reason for this is that the stiffness matrix $D$ from the TO approach does not have the same banded structure as that coming from the CG approach, making a solution of the eigenproblem more expensive (for the example of the cylinder flow: 1h 28min for the non-adaptive TO using a $P^{1}$ $513\times 513$ mesh vs. 14 minutes for the CG approach with $P^{2}$ elements on a $1025\times 1025$ mesh). The reference solution here is the same as that used for the CG experiments.

4.1.1 Standard Map

We investigate the convergence of collocation-based TO approaches on the Standard map test case.

We see an improvement by a full order of magnitude in the errors in both the eigenvalues and the eigenspaces for $P^{2}$ elements over $P^{1}$ elements. Concerning the convergence order, this experiment suggests that the order for $P^{2}$ elements is not higher than for $P^{1}$ elements in the TO approach, which match those obtained for eigenvectors based on $P^{1}$ elements in the CG approach.

4.1.2 Cylinder Flow

We observe a smaller error reduction in the eigenvalue and eigenvector errors for $P^{2}$ vs $P^{1}$ elements for the cylinder flow (Fig. 6) compared to the standard map example. The reduction in convergence orders compared to the standard map (i.e. comparing Fig. 6 to Fig. 5) mirrors the reductions observed in the corresponding CG experiments (i.e. comparing Fig. 2 to Fig. 1), where the orders are around half those seen for the standard map. Similarly to the standard map results, the $P^{2}$ elements did not perform asymptotically better than $P^{1}$ elements.

4.1.3 Bickley Jet

We repeat the experiment from Section 3.2 with the (non-adaptive) TO approach, see  Fig. 7. As in Section 3.2, the mesh widths in Fig. 7 below are the lowest for which the topology of the clustering result does not change. Here, the dominant computational cost is given by the evaluation of the flow map, which has to be evaluated once for each basis function of the finite element space. Note that there are three basis functions per element for $P^{1}$ and six per element for $P^{2}$ elements. While for $P^{2}$ elements, a coarser mesh was sufficient, the number of basis functions was comparable (even a little larger) than for the $P^{1}$ case – which is in contrast to the corresponding CG experiment.

4.2 TO Convergence Rates

Determining the theoretical convergence orders of the TO-methods is an outstanding task. In order to investigate whether we can expect higher-order convergence, we consider the collocation-based non-adaptive TO approach on a simple example in one dimension. In this example we can largely isolate the dynamics from the errors to focus on errors arising from translations of the basis functions.

The only volume-preserving diffeomorphisms of the circle $\mathbb{S}^{1}=\mathbb{R}/\mathbb{Z}$ to itself are rigid rotations. Let $T:\mathbb{S}^{1}\rightarrow\mathbb{S}^{1}$ be given by $T(x)=(x+\alpha)\mod 1$. Rigid rotations commute with the Laplace operator; thus $\Delta^{dyn}=(\Delta+T^{*}_{1}\Delta T_{1,*})/2=\Delta$, and the dynamic Laplacian is equal to the static Laplacian. The first nontrivial eigenspace of the dynamic Laplacian is therefore spanned by $v(x)=\sin(2\pi x)$ and $u(x)=\cos(2\pi x)$ which are orthogonal eigenvectors for the eigenvalue $\lambda_{1}=4\pi^{2}$.

Any errors in the corresponding discrete bilinear form

[TABLE]

therefore arise solely from discretisation errors related to the rotation of the $\varphi_{i}$ by $\alpha$. In our experiments, we consider $\alpha=0.15$. For the non-adaptive collocation TO approach, denote the leading nontrivial eigenvector of $\Delta^{dyn}$ by $\tilde{v}_{h,1}$. We now look at the order of convergence for the error $\|\tilde{v}_{h,1}-v\|_{L^{2}}$ as $h\rightarrow 0$. We compute the errors independently to the numerical experiments by using Fourier coefficients. With the Fourier coefficients $A_{h}=\langle v,\tilde{v}_{1,h}\rangle_{L^{2}}$ and $B_{h}=\langle u,\tilde{v}_{1,h}\rangle_{L^{2}}$ of $\tilde{v}_{h,1}$, $Pv_{h,1}=A_{h}v+B_{h}u$ is the $L^{2}$-orthogonal projection of $\tilde{v}_{1,h}$ onto the first nontrivial eigenspace. As a measure of eigenvector error we now choose $e_{h}=\|P\tilde{v}_{h,1}-\tilde{v}_{1,h}\|_{L^{2}}^{2}$. All integrals were approximated by the trapezoid rule with $10^{8}$ quadrature points on $[0,1]$; see Figure 8.

The convergence rates in Figure 8 are identical to those found for the standard map in Figure 1 and Figure 5. The CG approach in this example is simply computing the eigenvectors and eigenvectors of the standard Laplace operator and therefore unsurprisingly one recovers the theoretical orders of convergence. These numerical results suggest that when using the non-adaptive TO method, we cannot expect an asymptotically higher convergence rate for $P^{2}$ elements in comparison to $P^{1}$-Elements even for very simple flow maps.

5 Conclusion

We compared the use of $P^{1}$ and $P^{2}$ elements in collocation-based CG and TO approximations of the dynamic Laplacian. In the CG approach applied to weakly nonlinear dynamics, $P^{2}$ elements can significantly reduce the computational cost by providing an asymptotically higher order of convergence. A benefit of $P^{1}$ elements for the CG approach is that a first-order method of quadrature can be used in the discretization, whereas using first-order quadrature for $P^{2}$ elements results in a singular mass matrix. This does not affect the asymptotic order of convergence, but nevertheless introduces a constant factor which may be relevant in some low-data cases.

In the non-adaptive TO approach, there seems to be little benefit gained by using $P^{2}$ elements as opposed to $P^{1}$ elements for the non-adaptive TO approximation. The adaptive TO is inherently $P^{1}$-based, and thus does not benefit from a $P^{2}$ discretization either. In general, numerical experiments suggest that collocation-based TO approaches have equal (or in some cases better) rates of convergence compared to the $P^{1}$ CG discretisation.

The $P^{2}$ CG discretisation generally had a higher order of convergence, though the numerically observed convergence rates varied significantly when the quadrature order was changed or when applied to more nonlinear and hyperbolic dynamics. It is difficult to compare the CG and TO approaches in general as the former relies on being able to calculate derivatives of the flow map, whereas the latter is purely data based. This makes the TO approach applicable to some cases where the CG approach cannot be used.

A hindrance to using the TO approach for finely resolved meshes is the fact that here the call to eigs takes much longer compared to when one uses the CG discretisation on the same mesh. We suspect that this is due to the fact that unlike in the TO approach, the CG discretisation preserves the banded structure of the stiffness matrix. More work is needed to determine how the eigenproblem can be solved efficiently in this case, or whether it is possible to avoid the eigenproblem completely but still be able to compute coherent sets. More work is also needed to determine the true rates of convergence of eigenvalues and eigenvectors. We proved that they do converge for the non-adaptive TO approach, but were only able to conjecture what the true rates are based on numerical experiments. It also remains to be seen whether Galerkin TO approaches can be modified to be computationally efficient.

6 Acknowledgements

OJ and NS acknowledge support by the Priority Programme SPP 1881 Turbulent Superstructures of the Deutsche Forschungsgemeinschaft and a Universities Australia / DAAD travel grant. GF is partially supported by an ARC Discovery Project and a joint Universities Australia / DAAD travel award. GF thanks the Faculty of Mathematics at the Technical University Munich for hospitality during his research visits. We would like thank Christian Ludwig for the reference to [4] used in the proof of Lemma 1. We also thank Alavaro de Diego and Daniel Karrasch for helpful comments and their contributions to the CoherentStructures.jl package.

Appendix A Proofs

Throughout all proofs, $C$ refers to a constant depending only on the mesh and dimension.

A.1 Proof of Lemma 1

The following proof of Lemma 1 uses ideas from the proof of Lemma 2.1 from [4]. We start with a helper lemma:

Lemma 9**.**

Let $\mathcal{T}^{0}_{h^{\prime}}$ and $\mathcal{T}^{1}_{h}$ be quasi-uniform meshes on $\Omega_{0}$ and $\Omega_{1}$ respectively for some values of $h,h^{\prime}>0$. Let $F:\Omega_{0}\rightarrow\Omega_{1}$ be smooth on $\Omega_{0}$ with smooth extension to the boundary. Then there exists $C>0$ independent of $h$ and $h^{\prime}$ so that for all $\tau^{\prime}\in\mathcal{T}^{0}_{h^{\prime}}$ it holds that:

[TABLE]

Proof.

By quasi-uniformity of the meshes, we can find positive constants $a^{\prime},A^{\prime}$ so that for any triangle $\tau^{\prime}\in\mathcal{T}^{0}_{h^{\prime}}$ there is a point $\tau^{\prime}_{x}$ so that:

[TABLE]

This gives that where $y:=F(\tau^{\prime}_{x})$:

[TABLE]

where the last inclusion follows from the mean-value theorem, and ${\|DF\|}_{\infty}$ exists as $F$ is smoothly extensible to the boundary.

Now write $a,A$ for the shape-regularity/quasi-uniformity constants of $\mathcal{T}_{h}^{1}$, so that a formula like (12) holds for $\mathcal{T}_{h}^{1}$ also.

Then if $\tau\in\mathcal{T}_{h}^{1}$ intersects $B(y,Eh^{\prime})$, it must hold that $B(\tau_{x},ah)\subset B(\tau_{x},Ah)\subset B(y,Eh^{\prime}+2Ah)\subset B(y,c(h^{\prime}+h))$ with suitable $c>0$. As for different $\tau$, the sets $B(\tau_{x},ah)$ are disjoint, we get

[TABLE]

which gives the claim. ∎

Lemma 10**.**

Let $\mathcal{T}_{h}$, $h>0$ be a family of quasi-uniform meshes on an open subset of $d$-dimensional space. Let $S_{h}$ be the space of functions representable by $P^{k}$-Lagrange elements on the mesh. Then there exist $C,C^{\prime}>0$ so that for all $v\in S^{h}$:

[TABLE]

where $p$ ranges over the nodes of the triangulation.

Proof.

This follows directly from well-known results about the spectrum of mass-matrices, see [9, p.386] ∎

Proof of Lemma 1.

Throughout the proof, $C$ refers to a constant that does not depend on $h$ or $h^{\prime}$ and whose exact value can change from line to line. Let $h^{\prime}\geq h>0$ and $v\in S^{0}_{h^{\prime}}$. Then by Lemma 10

[TABLE]

where $p$ ranges over the nodes of $\mathcal{T}^{1}_{h}$. Using Lemma 9, we know that for any triangle $\tau^{\prime}\in\mathcal{T}^{0}_{h}$, at most $C\left(\frac{h^{\prime}+h}{h}\right)^{d}$ triangles from $\mathcal{T}^{1}_{h^{\prime}}$ can intersect with $F(\tau^{\prime})$, up to a constant factor this therefore bounds the number of vertices $p$ in $\mathcal{T}^{1}_{h}$ for which $F^{-1}(p)$ lies in a given triangle. Moreover, as we are using $P^{1}$-Lagrange elements, $|v(F^{-1}(p))|$ is bounded by $|v(p^{\prime})|$ for some vertex $p^{\prime}$ of the triangle that contains $F^{-1}(p)$. This gives, using Lemma 10 and the fact that $h^{\prime}\geq h$ that:

[TABLE]

where $p^{\prime}$ ranges over the nodes of $\mathcal{T}^{0}_{h^{\prime}}$. ∎

A.2 Proof of Lemma 2

Proof.

Without loss of generality, look at triangles in dimension $2$. Since ${\|I_{h}v-v\|}_{L^{2}}\leq C{\|I_{h}v-v\|}_{L^{\infty}}$ and $v$ is uniformly continuous, we have that

[TABLE]

What remains to be shown is ${\|\nabla(I_{h}v-v)\|}_{L^{2}}\rightarrow 0$.

We first show that ${\|\nabla I_{h}v\|}_{L^{\infty}}\leq C{\|\nabla v\|}_{L^{\infty}}$ (where $C$ does not depend on $h$). It is enough to prove this for any triangle $\tau$ in the mesh $\mathcal{T}^{1}_{h}$. Assume first that the triangle $\tau$ has vertices $0,e_{1}$ and $e_{2}$, without loss of generality also $v(0)=0$. Then $\nabla I_{h}v=(v(e_{1}),v(e_{2}))^{T}$.

By a mean value inequality555As $v$ is only piecewise $C^{\infty}$, we cannot use the mean value theorem directly. However, as we only care about a weaker statement of the form $|f(1)-f(0)|\leq{\|f^{\prime}\|}_{L^{\infty}([0,1])}$ this is not a problem. , $\max\{|v(e_{1})|,|v(e_{2})|\}\leq{\|\nabla v\|}_{L^{\infty}(\tau)}$. Thus ${\|\nabla I_{h}v\|}_{L^{\infty}(\tau)}\leq{\|\nabla v\|}_{L^{\infty}(\tau)}$. Shape regularity immediately gives ${\|\nabla I_{h}v\|}_{L^{\infty}(\tau)}\leq C{\|\nabla v\|}_{L^{\infty}(\tau)}$ for general triangles, taking suprema over all triangles gives the claim.

By assumption $v$ is $C^{2}$ except for on a nowhere dense set of measure zero. By standard FEM theory, ${\|\nabla(I_{h}v-v)\|}_{L^{2}(\Omega^{\prime})}\rightarrow 0$ on all sub-meshes $\Omega^{\prime}$ on which $v$ is $C^{\infty}$, and with $h\rightarrow 0$ we can choose $\Omega^{\prime}$ so that $\lambda^{d}(\Omega^{\prime})\rightarrow\lambda^{d}(\Omega)$. Moreover, $|\nabla I_{h}v|<C$ almost everywhere, and hence discontinouities do not cause problems. It follows that ${\|\nabla(I_{h}v-v)\|}_{L^{2}}\rightarrow 0$. ∎