Linear response for the dynamic Laplacian and finite-time coherent sets

TL;DR

This paper develops a linear response theory for the eigenfunctions of the dynamic Laplace operator to understand how finite-time coherent sets in nonlinear dynamical systems change under small parameter variations, with efficient numerical methods demonstrated.

Contribution

It introduces a novel linear response framework for the dynamic Laplace operator's eigenfunctions, enabling analysis of coherent set sensitivity to parameter changes.

Findings

Linear response theory for dynamic Laplace eigenfunctions established.

Numerical methods based on finite-element approach developed.

Numerical examples demonstrate practical applicability.

Abstract

Finite-time coherent sets represent minimally mixing objects in general nonlinear dynamics, and are spatially mobile features that are the most predictable in the medium term. When the dynamical system is subjected to small parameter change, one can ask about the rate of change of (i) the location and shape of the coherent sets, and (ii) the mixing properties (how much more or less mixing), with respect to the parameter. We answer these questions by developing linear response theory for the eigenfunctions of the dynamic Laplace operator, from which one readily obtains the linear response of the corresponding coherent sets. We construct efficient numerical methods based on a recent finite-element approach and provide numerical examples.

Equations

$${\partial }_{t}x=v(t,x,\epsilon ),$$

$${\partial }_t{\partial }_{\epsilon }{\phi }_{\epsilon }^{t_0,t}(x){|}_{\epsilon =0}={\partial }_{x}v(t,T_0(x),0){\partial }_{\epsilon }{\phi }_{\epsilon }^{t_0,t}(x){|}_{\epsilon =0}+{\partial }_{\epsilon }v(t,T_0(x),0).$$

$${\partial }_{t}x=v(x,t)$$

$$T_{\epsilon }(x):={\phi }^{t_0,t_1+\epsilon }(x)={\phi }^{t_0,t_1}(x)+\epsilon {\partial }_t{\phi }^{t_0,t}{|}_{t=t_1}(x)+o(\epsilon )=T_0(x)+\epsilon {\dot{T}}_0(x)+o(\epsilon ),$$

$$\|f\|_{C^{2}(\Omega,\mathbb{R}^{n})}=\sum_{\alpha_{j},|\alpha|\leq 2}\max_{x\in\Omega}\bigg{\|}\frac{\partial^{|\alpha|}f}{\partial x^{\alpha_{1}}_{1}\dots\partial x^{\alpha_{n}}_{n}}(x)\bigg{\|},$$

$$T_{\epsilon }=T_0+\epsilon {\dot{T}}_0+R_{\epsilon }$$

$$\mathbf{h}(\Gamma ):=\frac{\frac{1}{2}({\ell }_{n-1}(\Gamma )+{\ell }_{n-1}(T_{\epsilon }(\Gamma )))}{\min \{\ell ({\Omega }_1),\ell ({\Omega }_2)\}},$$

$$\mathbf{h}:=\min \{\mathbf{h}(\Gamma ):\Gamma \text{ is a }{C}^{\infty }\text{ codimension 1}\text{ manifold disconnecting }\Omega \}.$$

$${\Delta }_{\epsilon }^D:=\frac{1}{2}\left({\Delta }_{\Omega }+T_{\epsilon }^{\ast }{\Delta }_{T_{\epsilon }(\Omega )}{T}_{\epsilon ,\ast }\right).$$

$$A_{\epsilon }:=\frac{1}{2}\left(I+(D{T}_{\epsilon }^{\top }D{T}_{\epsilon }{)}^{-1}\right).$$

$${\Delta }_{\epsilon }^D{u}_{\epsilon }={\lambda }_{\epsilon }{u}_{\epsilon }\text{ on int}(\Omega ),$$

$$\nabla u_{\epsilon }\bullet A_{\epsilon }\nu =0(\text{resp.}{u}_{\epsilon }=0)\text{ on }\partial \Omega$$

$$-\frac{1}{2}\left({\int }_{\Omega }\nabla u_{\epsilon }\bullet \nabla \phi d\ell +{\int }_{T_{\epsilon }(\Omega )}\nabla (T_{\epsilon ,\ast }{u}_{\epsilon })\bullet \nabla (T_{\epsilon ,\ast }\phi )d\ell \right)={\lambda }_{\epsilon }{\int }_{\Omega }{u}_{\epsilon }\phi d\ell \forall \phi \in H,$$

$$-{\int }_{T_{\epsilon }(\Omega )}\nabla (T_{\epsilon ,\ast }{u}_{\epsilon })\bullet \nabla (T_{\epsilon ,\ast }\phi )d\ell$$

$$-a_{\epsilon }(u_{\epsilon },\phi )={\lambda }_{\epsilon }⟨u_{\epsilon },\phi ⟩\forall \phi \in H.$$

$${\dot{u}}_0:=\frac{d}{d\epsilon }{u}_{\epsilon }{|}_{\epsilon =0}.$$

$${\dot{A}}_0=-{\left((D{T}_0{)}^{-1}(D{\dot{T}}_0)(D{T}_0{)}^{-1}(D{T}_0{)}^{-\top }\right)}^{\text{sym}}$$

$$\lim_{\varepsilon\rightarrow 0}\bigg{\|}\frac{A_{\varepsilon}-A_{0}}{\varepsilon}-\dot{A}_{0}\bigg{\|}_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}=0.$$

$$(D{T}_{\epsilon }{)}^{-1}$$

$$\Vert (D{T}_0{)}^{-1}(\epsilon D{\dot{T}}_0+D{R}_{\epsilon }){\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}$$

$$\le \Vert (D{T}_0{)}^{-1}{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}\left(|\epsilon |\Vert D{\dot{T}}_0{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}+\Vert D{R}_{\epsilon }{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}\right)$$

$$\le |\epsilon |\Vert (D{T}_0{)}^{-1}{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}\left(\Vert D{\dot{T}}_0{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}+C\right).$$

$$|\epsilon |<{\left(\Vert (D{T}_0{)}^{-1}{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}\left(\Vert D{\dot{T}}_0{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}+C\right)\right)}^{-1}$$

$$(D{T}_{\epsilon }{)}^{-1}$$

$$\Vert {\hat{R}}_{\epsilon }{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}$$

$$+\sum _{i\ge 2}\Vert (D{T}_0{)}^{-1}{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}^i{\left(\epsilon \Vert D{\dot{T}}_0{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}+\Vert D{R}_{\epsilon }{\Vert }_{C^1(\Omega ,\mathcal{B}({\mathbb{R}}^n))}\right)}^i,$$

$$(D{T}_{\epsilon }{)}^{-1}(D{T}_{\epsilon }{)}^{-\top }$$

$$=(D{T}_0{)}^{-1}(D{T}_0{)}^{-\top }$$

$$-\epsilon \left((D{T}_0{)}^{-1}(D{\dot{T}}_0)(D{T}_0{)}^{-1}(D{T}_0{)}^{-\top }+(D{T}_0{)}^{-1}(D{T}_0{)}^{-\top }(D{\dot{T}}_0{)}^{\top }(D{T}_0{)}^{-\top }\right)+{\tilde{R}}_{\epsilon }$$

$$=(D{T}_0{)}^{-1}(D{T}_0{)}^{-\top }+2\epsilon {\dot{A}}_0+{\tilde{R}}_{\epsilon },$$

$${\dot{a}}_0(\psi ,\phi ):={\int }_{\Omega }{\dot{A}}_0\nabla \psi \bullet \nabla \phi d\ell .$$

$${\dot{a}}_0(\psi ,\phi )=\underset{\epsilon \to 0}{\lim }\frac{a_{\epsilon }(\psi ,\phi )-a_0(\psi ,\phi )}{\epsilon }.$$

$$a_{\epsilon }(\psi ,\phi )-a_0(\psi ,\phi )-\epsilon {\dot{a}}_0(\psi ,\phi )={\int }_{\Omega }(A_{\epsilon }-A_0-\epsilon {\dot{A}}_0)\nabla \psi \bullet \nabla \phi d\ell .$$

Full text

1 Introduction

Finite-time coherent sets [19, 10, 11] are regions in the compact phase space of a nonlinear dynamical system that minimally mix over a finite time duration, and therefore play an important role in the analysis of how material objects are transported in fluids. Spectral methods for identifying finite-time coherent sets were developed in [19, 10] directly from transfer operators, and later in [11, 15] using the dynamic Laplacian, which was derived as a zero-diffusion limit of the transfer operator construction in [10]. Various implementations of these two approaches and related methods may be found in [29, 17, 12, 35, 7, 24, 5, 26, 13, 9, 14]. In the present paper we use the approach of [11, 15], which defines finite-time coherent sets through the notion of dynamic isoperimetry: those sets whose boundary size relative to volume remains small under the finite-time dynamics. These persistently small boundaries represent evolving fluid interfaces across which there is minimal mixing. The key technology for finding these coherent sets is the dynamic Laplace operator, defined in [11]; the leading eigenfunctions of this operator encode the finite-time coherent sets. The dynamic Laplacian and its eigenfunctions may be efficiently approximated using a specialised finite element method [13], and individual coherent sets may be automatically separated using algorithms such as SEBA [18].

Throughout, we will represent the finite-time dynamics by a single application of a transformation $T_{0}:\Omega\to T_{0}(\Omega)$, $\Omega\subset\mathbb{R}^{n}$ compact; $T_{0}$ may arise, for example, as a flow map of some nonautonomous ordinary differential equation. The question we investigate in this work is how coherent sets behave under perturbation of the dynamics. For some $\varepsilon_{0}>0$ we consider a family of maps $T_{\varepsilon}$, $\varepsilon\in(-\varepsilon_{0},\varepsilon_{0})$, where we think of $T_{0}$ as governing the original, unperturbed dynamics. As $\varepsilon$ is varied from zero, the dynamic Laplacian corresponding to $T_{\varepsilon}$, its eigenfunctions, and the corresponding coherent sets, all vary from those objects computed with $T_{0}$. Given sufficient regularity of $\varepsilon\mapsto T_{\varepsilon}$ we may hope for some regular dependence of coherent sets on $\varepsilon$. The notion of coherent sets has found application in fluid flows from the laboratory scale to the planetary scale, and the dynamic Laplacian has proven to be an efficient way of extracting coherent objects (such as the Gulf Stream and ocean eddies [18]). In the context of system perturbations due to climate change, an important step in quantifying potential impacts would be the prediction of responses of coherent geophysical features.

The approach we take is based on the ideas of linear response [32, 2] in dynamical systems. Linear response is classically concerned with the derivatives of physical invariant measures $\mu_{\varepsilon}$ of autonomous maps $T_{\varepsilon}$ with respect to the parameter $\varepsilon$. The physical invariant measure $\mu_{\varepsilon}$ is the leading eigenfunction (or eigendistribution) of the transfer operator $\mathcal{L}_{\varepsilon}$ for $T_{\varepsilon}$, and formulae for $\frac{\partial}{\partial\varepsilon}\mu_{\varepsilon}$ involve $\frac{\partial}{\partial\varepsilon}\mathcal{L}_{\varepsilon}$. In order for physical invariant measures to exist, usually some hyperbolicity of the dynamics is required. For Anosov maps (and more general dynamical systems, like Axiom A), the differentiability (which include the linear response) of the eigendata of the transfer operator associated to these dynamical systems have been obtained [21, 22]; linear response results are also available for uniformly hyperbolic flows [6, 33]. Aside from smooth dynamics, linear response has also been treated for unimodal maps [3] and intermittent maps [4]; there are also results for the existence of linear response for stochastic dynamical systems [25, 20, 1]. Linear response is not guaranteed; see e.g. [2] for details on failure of linear response.

In the present paper, we wish to follow this general notion of linear response, namely computing the derivative of a quantity with respect to a parameter. Specifically, we replace a family of transfer operators $\mathcal{L}_{\varepsilon}$ with a family of dynamic Laplace operators arising from a family of diffeomorphisms $T_{\varepsilon}$. This new linear response problem requires a very different functional analytic setup and has rather well-behaved responses to perturbations. We prove that differentiability of $\varepsilon\mapsto T_{\varepsilon}$ implies differentiability of the dynamic Laplacian. Further, if a particular eigenvalue is algebraically simple, this eigenvalue and the associated eigenfunction have a linear response (are differentiable with respect to $\varepsilon$). We obtain a formula for the derivative of the eigenvalues with respect to the parameter; the derivative of the second eigenvalue quantifies the instantaneous rate of change of global mixing as the parameter is varied. We then derive a formula for the linear response of the eigenfunctions; the derivatives of the dominant eigenfunctions of the dynamic Laplacian immediately yield derivatives of the corresponding finite-time coherent sets. Building on the finite-element method (FEM) based approaches in [13] we develop numerical schemes for numerically computing these linear responses, and illustrate these schemes on the standard map and the Meiss-Mosovsky map. In addition to computing the response of coherent sets, we observe that our first-order approximations of the perturbed eigenvectors, computed using linear response, produce coherent sets that are rather close to the true coherent sets, even for large extrapolation values.

An outline of the paper is as follows: in Section 2 we introduce differentiability hypotheses on the dynamics. In Section 3 we define the dynamic Laplacian, coherent sets, and our linear response problem. Section 4 contains the proof of the weak differentiability of the dynamic Laplacian with respect to the perturbing parameter, and the proof of the existence of linear response of eigenvectors. In Section 5 we derive a linear system whose solution provides the linear response, and in Section 6 we develop two FEM-based approaches to numerically solve this linear system and estimate the linear responses. We conclude in Section 7 with numerical demonstrations of the theory.

2 Perturbations

Let $\Omega\subset\mathbb{R}^{n}$ be a compact, connected domain with smooth boundary. We consider a family of maps $T_{\varepsilon}:\Omega\rightarrow T_{\varepsilon}(\Omega)$, $\varepsilon\in(-\varepsilon_{0},\varepsilon_{0})$, where $T_{0}$ represents the original, unperturbed dynamics. For simplicity, we assume that $T_{\varepsilon}$ is volume-preserving, and consider a single application of $T_{\varepsilon}$. The methods we propose are easily extendable to non-volume-preserving $T_{\varepsilon}$, curved manifolds [15], and multiple applications of $T_{\varepsilon}$ (see [11]).

Special families we have in mind are:

$T_{\varepsilon}$ is given by the flow map $\varphi_{\varepsilon}^{t_{0},t_{1}}$ of some ordinary differential equation

[TABLE]

where the vector field $v$ depends on a parameter $\varepsilon$ and $t_{0},t_{1}\in\mathbb{R}$ are chosen such that the flow map is defined for all $x$. Under appropriate assumptions on $v$ we have $T_{\varepsilon}=T_{0}+\varepsilon\dot{T}_{0}+o(\varepsilon)$, where $\dot{T}_{0}(x):=\partial_{\varepsilon}\varphi_{\varepsilon}^{t_{0},t_{1}}(x)|_{\varepsilon=0}$ and $\partial_{\varepsilon}\varphi_{\varepsilon}^{t_{0},t}(x)|_{\varepsilon=0}$ satisfies the variational equation

[TABLE] 2. 2.

As a further specialisation of 1. we interpret the time $t$ itself as the parameter $\varepsilon$, i.e. we consider

[TABLE]

with the flow map $T_{0}=\varphi^{t_{0},t_{1}}$. In this case we have that

[TABLE]

where $\dot{T}_{0}(\cdot):={\partial_{t}}\varphi^{t_{0},t}|_{t=t_{1}}(\cdot)=v(\cdot,t_{1})$.

The precise setting we consider is the following: Let Diff${}^{2}(\Omega,\mathbb{R}^{n})$ be the space of $C^{2}$-diffeomorphisms from $\Omega$ to $\mathbb{R}^{n}$ which is endowed with the $C^{2}$-norm

[TABLE]

where $\alpha=(\alpha_{1},\dots,\alpha_{n})\in\mathbb{N}_{0}^{n},|\alpha|=\sum_{j=1}^{n}\alpha_{j}$. We assume that the map $\varepsilon\mapsto T_{\varepsilon}$ is $C^{1}$ from $(-\varepsilon_{0},\varepsilon_{0})\subset\mathbb{R}$ to $\text{Diff}^{2}(\Omega,\mathbb{R}^{n})$. From Taylor’s theorem (see [27], XIII §6) for sufficiently small $\varepsilon_{0}>0$, one has

[TABLE]

for $\varepsilon\in(-\varepsilon_{0},\varepsilon_{0})$, where $\dot{T}_{0},R_{\varepsilon}\in\text{Diff}^{2}(\Omega,\mathbb{R}^{n})$, and $\|R_{\varepsilon}\|_{C^{2}(\Omega,\mathbb{R}^{n})}=o(\varepsilon)$. Since all $T_{\varepsilon}$ are $C^{2}$ diffeomorphisms we have that for any $\varepsilon\in(-\varepsilon_{0},\varepsilon_{0})$ the maps $DT_{\varepsilon}$ and $DT^{-1}_{\varepsilon}$ are in $C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))$, where $\mathcal{B}(X)$ is the space of bounded linear maps from some Banach space $X$ into itself.

3 The dynamic Laplacian

We are ultimately interested in analysing the response of coherent sets to perturbations of the dynamics. As coherent sets can be detected via level sets of leading nontrivial eigenfunctions of the dynamic Laplace operator, we need to understand how these eigenfunctions respond to perturbations in the dynamics, i.e. how they change with $\varepsilon$.

Following [11], when dividing a manifold $\Omega$ into two coherent sets, one seeks a dynamically minimal interface $\Gamma$ disconnecting $\Omega$; the interface $\Gamma$ forms the shared boundary of the two coherent sets. More precisely, if $\Gamma\subset\Omega$ is a codimension-1 submanifold disconnecting $\Omega$ into $\Omega_{1}$ and $\Omega_{2}$, we compute the dynamic Cheeger value of $\Gamma$:

[TABLE]

where $\ell_{n-1}$ is the induced $n-1$-dimensional volume and $\ell$ is the $n-$dimensional volume. We seek the minimising $\Gamma$ to obtain the dynamic Cheeger constant [11]:

[TABLE]

In the case where we do not wish the interface $\Gamma$ to intersect the boundary of $\Omega$ we can alternatively consider a Dirichlet dynamic Cheeger constant; see §2.2 [13]. These two options are summarised in [13, Figures 2 and 3], respectively.

A minimizing $\Gamma$ can be linked to level sets of eigenfunctions of a dynamic Laplace operator, see [11, 15, 13]. Denote the pushforward resp. pullback of a function $f:\Omega\to\mathbb{R}$ by $T_{\varepsilon,*}f:=f\circ T_{\varepsilon}^{-1}$ resp. $T_{\varepsilon}^{*}f:=f\circ T_{\varepsilon}$ and let $\Delta_{\Omega}$ be the Laplace operator on $\Omega$. The dynamic Laplace operator [11] is

[TABLE]

Define the matrix-valued function $A_{\varepsilon}\in C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))$ by

[TABLE]

We are interested in the eigenproblem

[TABLE]

with homogeneous Neumann (resp. Dirichlet) boundary conditions

[TABLE]

($\nu$ denotes the outer normal on $\Omega$). The spectral properties of the family $\Delta^{D}_{\varepsilon}$ are given by Theorem 4.1 [11]. A discussion of the interpretation of the (natural) Neumann boundary conditions is given immediately after Theorem 3.2 [11]; the Dirichlet boundary condition case is developed in [13]. Throughout the paper, we will assume that all eigenvalues of $\Delta^{D}_{\varepsilon}$ are algebraically simple.

The weak form of the eigenvalue problem (7)–(8) is given by

[TABLE]

where $H$ denotes $H^{1}(\Omega)$ in the case of Neumann and $H^{1}_{0}(\Omega)$ in the case of homogenous Dirichlet boundary conditions. Note that if we let $\varphi=u_{\varepsilon}$ in (9), all integrals are positive; thus, the eigenvalues $\lambda_{\varepsilon}$ are negative (or [math]). Note further that

[TABLE]

so that (9) can be written as

[TABLE]

with the bilinear form $a_{\varepsilon}(u_{\varepsilon},\varphi)=\int_{\Omega}A_{\varepsilon}\nabla u_{\varepsilon}\bullet\nabla\varphi\;d\ell$, and $\langle\cdot,\cdot\rangle$ the $L^{2}$ scalar product on $\Omega$. By the above considerations we may also write $\Delta^{D}_{\varepsilon}=\sum_{i,j=1}^{n}\partial_{j}A_{\varepsilon,ij}\partial_{i}$.

4 Existence of a Linear Response

Throughout we assume that $u_{\varepsilon}$ is scaled so that $\|u_{\varepsilon}\|=1$, where $\|\cdot\|$ is the $L^{2}(\Omega,\ell)$ norm. In order to answer the question of how coherent sets of $T_{\varepsilon}$ depend on $\varepsilon$, we are going to show that the map $\varepsilon\mapsto u_{\varepsilon}$ is differentiable at [math] as a map from $(-\varepsilon_{0},\varepsilon_{0})$ to $H$ and devise a method for computing the linear response

[TABLE]

We begin with a lemma about the regularity of the coefficient function $A_{\varepsilon}$ of the dynamic Laplace operator. In Proposition 1 we show that we can differentiate (in a weak sense) the map $\varepsilon\mapsto\Delta^{D}_{\varepsilon}$. Finally, we apply a general regularity theorem for the spectral data of elliptic operators to obtain the differentiability of the maps $\varepsilon\mapsto\lambda_{\varepsilon}$ and $\varepsilon\mapsto u_{\varepsilon}$.

Let $Q^{\text{sym}}=\frac{1}{2}(Q+Q^{\top})$ denote the symmetric part of a matrix $Q$.

Lemma 1**.**

The matrix-valued function $\dot{A}_{0}\in C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))$ given by

[TABLE]

satisfies

[TABLE]

Proof.

We recall from (2) that for sufficiently small $\varepsilon_{0}>0$, we have that $T_{\varepsilon}=T_{0}+\varepsilon\dot{T}_{0}+R_{\varepsilon}$ for $\varepsilon\in(-\varepsilon_{0},\varepsilon_{0})$ with $\|R_{\varepsilon}\|_{C^{2}(\Omega,\mathbb{R}^{n})}=o(\varepsilon)$, yielding

[TABLE]

Using the fact that $\|R_{\varepsilon}\|_{C^{2}(\Omega,\mathbb{R}^{n})}=o(\varepsilon)$, we have that $\|DR_{\varepsilon}\|_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}=o(\varepsilon)$ and so there exists $C<\infty$, that is independent of $\varepsilon$, such that

[TABLE]

Choosing $\varepsilon$ small enough to satisfy

[TABLE]

we get $\|(DT_{0})^{-1}(\varepsilon D\dot{T}_{0}+DR_{\varepsilon})\|_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}<1$. We can now use the Neumann series representation for the RHS of (13) to obtain

[TABLE]

where $\widehat{R}_{\varepsilon}=(DT_{0})^{-1}DR_{\varepsilon}+\sum_{i\geq 2}(-1)^{i}\left((DT_{0})^{-1}(\varepsilon D\dot{T}_{0}+DR_{\varepsilon})\right)^{i}$. Noting that

[TABLE]

and using the fact that $\|DR_{\varepsilon}\|_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}=o(\varepsilon)$, we have that $\|\widehat{R}_{\varepsilon}\|_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}=o(\varepsilon)$. Hence, using (14) we get

[TABLE]

where $\|\widetilde{R}_{\varepsilon}\|_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}=o(\varepsilon)$. We conclude that $\|A_{\varepsilon}-A_{0}-\varepsilon\dot{A}_{0}\|_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}=o(\varepsilon)$. ∎

We define $\dot{\Delta}_{0}^{D}:=\mathop{\rm div}(\dot{A}_{0}\nabla)$ and consider the associated bilinear form

[TABLE]

Proposition 1**.**

The bilinear form (15) is a weak derivative of the weak form of $\Delta_{\varepsilon}^{D}$ at $\varepsilon=0$ in the sense that for $\psi,\varphi\in H$

[TABLE]

Proof.

We have

[TABLE]

Lemma 1 yields $A_{\varepsilon}=A_{0}+\varepsilon\dot{A}_{0}+R_{\varepsilon}$ with $\|R_{\varepsilon}\|_{C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))}=o(\varepsilon)$. We therefore immediately get

[TABLE]

∎

We now state a theorem concerning differentiability of the spectral data for the eigenproblem:

[TABLE]

of some general uniformly elliptic second order differential operator $L_{A}=\sum_{i,j=1}^{n}\partial_{j}A_{ij}\partial_{i}$ with coefficients $A=(A_{ij})$. Let $\Lambda(L_{A})\subset\mathbb{R}\times H$ be the set of eigenpairs $(\lambda,u)$ of $L_{A}$.

Theorem 1** ([23]).**

Let $O\subset\mathbb{R}^{n}$ be a bounded domain and $A_{0}\in C^{k}(\overline{O},\mathbb{R})^{n^{2}}$, $k\geq 1$, the coefficients of the uniformly elliptic operator $L_{A_{0}}$. Let $(\lambda_{0},u_{0})\in\Lambda(L_{A_{0}})$ and assume $\lambda_{0}$ is algebraically simple. Then there exists a neighbourhood $\mathcal{U}\subset C^{k}(\overline{O},\mathbb{R})^{n^{2}}$ of $A_{0}$ and $C^{k}$-functions $\bm{\lambda}:\mathcal{U}\rightarrow\mathbb{R}$ and $\mathbf{u}:\mathcal{U}\rightarrow H_{0}^{1}(O)$ such that:

$\bm{\lambda}(A_{0})=\lambda_{0}$* and $\mathbf{u}(A_{0})=u_{0}$;* 2. 2.

$(\bm{\lambda}(A),\mathbf{u}(A))\in\Lambda(L_{A})$* for every $A\in\mathcal{U}$.*

Let $(A_{\varepsilon,ij})$ be the entries of $A_{\varepsilon}$. Note that since $A_{\varepsilon}$ is in $C^{1}(\Omega,\mathcal{B}(\mathbb{R}^{n}))$, we have that $(A_{\varepsilon,ij})\in C^{1}(\Omega,\mathbb{R})^{n^{2}}$. We note further that $L_{A_{\varepsilon}}=\Delta_{\varepsilon}^{D}$ is uniformly elliptic [11, 15] and so this theorem applies to the eigenproblem (7)–(8) setting $\overline{O}=\Omega$. We note that the proof in [23] does not make use of the assumption of zero Dirichlet boundary data and in fact also applies to the Neumann boundary case.

In the subsequent results, $H$ denotes $H^{1}(\Omega)$ in the case of homogeneous Neumann boundary conditions, and $H^{1}_{0}(\Omega)$ in the case of homogeneous Dirichlet boundary conditions. The following theorem establishes the existence of derivatives of the maps $\varepsilon\mapsto u_{\varepsilon}$ and $\varepsilon\mapsto\lambda_{\varepsilon}$ from $(-\varepsilon_{0},\varepsilon_{0})$ to $H$. Let $V_{0}:=$ span$\{u_{0}\}^{\perp}\subset H$.

Theorem 2**.**

Let $\lambda_{0}$ be algebraically simple and $(\lambda_{\varepsilon},u_{\varepsilon})\in\Lambda(\Delta_{\varepsilon}^{D})$ for $\varepsilon\in(-\varepsilon_{0},\varepsilon_{0})$. Then there exists a function $\dot{u}_{0}\in H$ and $\dot{\lambda}_{0}\in\mathbb{R}$ such that

[TABLE]

Furthermore, $\dot{u}_{0}\in V_{0}$.

Proof.

Let $\mathcal{U}\ni A_{0}$ be the neighborhood and $\mathbf{u}:\mathcal{U}\rightarrow H$ and $\bm{\lambda}:\mathcal{U}\rightarrow\mathbb{R}$ the maps according to Theorem 1. Since these maps are $C^{1}$, there exist bounded linear maps $B_{1}:C^{1}(\Omega,\mathbb{R})^{n^{2}}\rightarrow H$ and $B_{2}:C^{1}(\Omega,\mathbb{R})^{n^{2}}\rightarrow\mathbb{R}$ satisfying

[TABLE]

and

[TABLE]

Define $\dot{u}_{0}:=B_{1}(\dot{A}_{0})\in H$ and ${\dot{\lambda}_{0}}:=B_{2}(\dot{A}_{0})$. Using $A_{\varepsilon,ij}=A_{0,ij}+\varepsilon\dot{A}_{0,ij}+r^{\varepsilon}_{ij}$ with $\|r^{\varepsilon}_{ij}\|_{C^{1}(\Omega,\mathbb{R})}=o(\varepsilon)$, the differentiability results follow.

In order to show that $\dot{u}_{0}\in V_{0}$, we note that for small $\varepsilon$, we have $u_{\varepsilon}=u_{0}+\varepsilon\dot{u}_{0}+g^{\varepsilon}$, where $g^{\varepsilon}\in H$ is such that $\|g^{\varepsilon}\|_{H}=o(\varepsilon)$. We therefore have

[TABLE]

thus, considering the leading term of order $\varepsilon$ we see that $\langle u_{0},\dot{u}_{0}\rangle=0$ and therefore $\dot{u}_{0}\in V_{0}$. ∎

5 A formula for the linear response

We will now derive a linear system that yields the linear response $\dot{u}_{0}$ as its solution. To this end, we first show that the (weak) derivative of the products $\varepsilon\mapsto\lambda_{\varepsilon}u_{\varepsilon}$ and $\varepsilon\mapsto\Delta_{\varepsilon}^{D}u_{\varepsilon}$ can be computed by the usual product rule.

Lemma 2**.**

For $\varphi\in H$,

[TABLE]

Proof.

From Theorem 2 we have that $u_{\varepsilon}=u_{0}+\varepsilon\dot{u}_{0}+g^{\varepsilon}$ and $\lambda_{\varepsilon}=\lambda_{0}+\varepsilon\dot{\lambda}_{0}+\mu^{\varepsilon}$, where $\|g^{\varepsilon}\|_{H}=o(\varepsilon)$ and $|\mu^{\varepsilon}|=o(\varepsilon)$. Thus,

[TABLE]

with $\|f^{\varepsilon}\|=o(\varepsilon)$, so that $|\langle f^{\varepsilon},\varphi\rangle|\leq\|f^{\varepsilon}\|\|\varphi\|=o(\varepsilon)$ for each $\varphi\in H$. This yields (19).

From Proposition 1, $a_{\varepsilon}(\psi,\varphi)=a_{0}(\psi,\varphi)+\varepsilon\dot{a}_{o}(\psi,\varphi)+o(\varepsilon)$ for all $\psi,\varphi\in H$. Hence,

[TABLE]

yielding (20).

∎

The following theorem establishes the existence of a unique solution of the linear system (25) in a weak sense.

Theorem 3**.**

Let $\dot{\lambda}_{0}$ and $\dot{u}_{0}$ be as in Theorem 2. These linear responses $(\dot{u}_{0},\dot{\lambda}_{0})\in V_{0}\times\mathbb{R}$ are the unique solution to the equations:

[TABLE]

and

[TABLE]

Proof.

We begin by showing that $\dot{\lambda}_{0}$ and $\dot{u}_{0}$ as in Theorem 2 solve (21) for all $\varphi\in V_{0}$ and (22). Subtract (19) from (20); we obtain 0 on the LHS because $u_{\varepsilon}$ is the eigenfunction associated to the eigenvalue $\lambda_{\varepsilon}$. Rearranging the RHS we immediately obtain that (21) is satisfied for all $\varphi\in H$.

We now write $H=\mathop{\rm span}\{u_{0}\}\oplus V_{0}$ and consider (21) for $\varphi$ according to this decomposition. Substituting $\varphi=u_{0}$ into (21) yields

[TABLE]

The LHS is zero since $u_{0}$ is the eigenfunction with eigenvalue $\lambda_{0}$; rearranging to solve for $\dot{\lambda}$ yields (22). Thus, (21) holding for all $\varphi\in H$ is equivalent to (21), (22), proving the statement, except for uniqueness.

Suppose that there is another pair $(\dot{w},\dot{\nu})\in V_{0}\times\mathbb{R}$ satisfying (21) for all $\varphi\in H$. Subtracting (21) with $(\dot{w},\dot{\nu})$ from (21) with $(\dot{u}_{0},\dot{\lambda}_{0})$ we obtain

[TABLE]

We again use the decomposition $H=\mathop{\rm span}\{u_{0}\}\oplus V_{0}$. Substituting $\varphi=u_{0}$, and arguing as previously, we see that the LHS of (24) is zero and therefore that $\dot{\nu}=\dot{\lambda}_{0}$, i.e.

[TABLE]

which implies that $\dot{u}_{0}-\dot{w}$ is a weak eigenfunction with eigenvalue $\lambda_{0}$. Because $\lambda_{0}$ is simple, we must have $\dot{u}_{0}-\dot{w}\in\mathop{\rm span}\{u_{0}\}$. Recalling that $\dot{u}_{0},\dot{w}\in V_{0}$, this implies that $\dot{u}_{0}-\dot{w}=0$. Thus with $\dot{\lambda}_{0}$ as in (22), there is a unique solution $\dot{u}_{0}$ to (21). ∎

We note that the strong form of (21) is given by the equation

[TABLE]

with boundary conditions

[TABLE]

In order to see this, multiply (25) with a test function $\varphi$ and apply the divergence theorem,this yields

[TABLE]

The boundary integrals either vanish if $\varphi\in H^{1}_{0}(\Omega)$ (the Dirichlet case) or if the (natural) boundary condition (26) is satisfied.

Remark 1**.**

We note that the expression (25) is reminiscent of the classical linear response formula for the invariant density of a deterministic dynamical system. In this setting, one has a family of transfer operators $\{\mathcal{L}_{\varepsilon}\}$ generated by a family of maps $\{T_{\varepsilon}\}$. The (typically assumed unique) fixed point $h_{\varepsilon}$ of $\mathcal{L}_{\varepsilon}$ is the invariant density of $T_{\varepsilon}$. It is easy to verify the identity $(I-\mathcal{L}_{\varepsilon})(h_{\varepsilon}-h_{0})=(\mathcal{L}_{\varepsilon}-\mathcal{L}_{0})h_{0}$. Dividing through by $\varepsilon$ and taking the limit as $\varepsilon\to 0$, one is able to show in certain situations that the limits $\dot{h}_{0}:=\lim_{\varepsilon\to 0}(h_{\varepsilon}-h_{0})/\varepsilon$ and $\dot{\mathcal{L}}_{0}:=\lim_{\varepsilon\to 0}(\mathcal{L}_{\varepsilon}-\mathcal{L}_{0})/\varepsilon$ exist in suitable senses, see e.g. [28]. This leads to $(I-\mathcal{L}_{0})\dot{h}_{0}=\dot{\mathcal{L}}_{0}h_{0}$, which is of the form (25) with $h_{0},\dot{h}_{0},\mathcal{L}_{0},\dot{\mathcal{L}}_{0}$ replaced by $u_{0},\dot{u}_{0},\Delta^{D}_{0},\dot{\Delta}^{D}_{0}$, respectively, noting that $\lambda_{0}=1$ and $\dot{\lambda}_{0}=0$.

6 Computing the linear response numerically

We now describe how to compute the linear response $\dot{u}_{0}$ numerically. To this end, we approximately solve the weak form (21) using the method described in [13]. That is, we consider (21) on a finite-dimensional approximation space $V_{N}\subset H$, denoting the approximations of $\lambda,\dot{\lambda},u_{0},\dot{u_{0}}$ by $\tilde{\lambda},\dot{\tilde{\lambda}},\tilde{u}_{0},\dot{\tilde{u}}_{0}$, respectively. Instead of choosing $V_{N}$ as a subspace of $V_{0}$ (as would be required by (21)), we enforce $\dot{\tilde{u}}_{0}\in\tilde{V}_{0}:=\mathop{\rm span}(\tilde{u}_{0})^{\bot}$ by adding an additional constraint. In practice, the approximation space will be realised as a finite element space, typically using linear triangular Lagrange elements.

In [13], two different variants of a finite-element discretisation of the basic eigenproblem for the dynamic Laplacian have been proposed, one based on the evaluation of the right Cauchy Green deformation tensor (the CG method) and one based on an explicit approximation of the transfer operator associated to $T_{\varepsilon}$ (the TO method). We now describe how to use both variants in order to compute $\dot{u}_{0}$.

6.1 The CG Method

Let $\varphi_{1},\ldots,\varphi_{N}$ be a basis for $V_{N}$. As described in [13] we obtain an approximation $\tilde{\lambda}_{0},\tilde{u}_{0}$ of the eigenpair $\lambda_{0},u_{0}$ by solving the matrix eigenproblem

[TABLE]

where

[TABLE]

are the stiffness and mass matrix, respectively, and ${\tilde{\bm{u}}}_{0}\in\mathbb{R}^{N}$ is the vector of coefficents of $\tilde{u}_{0}$ with respect to the chosen basis.

Similarly, we define the Galerkin approximation $\dot{\tilde{u}}_{0}\in V_{N}$ of $\dot{u}_{0}$ by requiring it to satisfy (21) for $\varphi=\varphi_{j},j=1,\ldots,N$. This yields the linear system

[TABLE]

for the coefficient vector $\dot{\tilde{\bm{u}}}_{0}$ of $\tilde{u}_{0}$ with respect to the basis $\varphi_{1},\ldots,\varphi_{N}$. Here,

[TABLE]

is the “linear response” matrix.

Instead of choosing $V_{N}$ as a subspace of $V_{0}$, we enforce $\dot{\tilde{u}}_{0}\in\tilde{V}_{0}:=\mathop{\rm span}(\tilde{u}_{0})^{\bot}$ by adding an additional constraint on the coefficient vectors ${\tilde{\bm{u}}}_{0}$ and $\dot{\tilde{\bm{u}}}_{0}$, namely

[TABLE]

which we append to (29). We combine (29) and (31) into a single linear system which allows to solve for $\dot{\tilde{u}}_{0}$ and $\dot{\tilde{\lambda}}_{0}$ simultaneously:

[TABLE]

Note that according to our standing assumption, $\lambda_{0}$ is simple and so $\tilde{\lambda}_{0}$ is simple if the elements are fine enough (cf. [8], Lemma 3.65 and [34]) and the kernel of the matrix $K-\tilde{\lambda}_{0}M$ is spanned by $\tilde{\bm{u}}_{0}$. Thus, on $\tilde{V}_{0}$, the matrix $K-\tilde{\lambda}_{0}M$ is nonsingular and ($K,M$ and $L$ are symmetric)

[TABLE]

i.e. the right hand side $(\dot{\tilde{\lambda}}_{0}M-L)\tilde{\bm{u}}_{0}$ of (29) is in $\tilde{V}_{0}$. The system (29) therefore has a unique solution on $\tilde{V}_{0}$ or equivalently:

Proposition 2**.**

The linear system (32) has a unique solution.

6.2 The TO Method

The second variant of the finite-element based computation of the linear response $\dot{\tilde{u}}$ employs an explicit approximation of the transfer operator $T_{0,*}$ associated to $T_{0}$. It yields an alternative way to compute the matrices $K$ and $L$ in (30) – everything else remains unchanged from the previous section. In particular, we again solve the linear system (32) in order to obtain the approximate linear response $\dot{\tilde{u}}_{0}$.

Approximating the transfer operator.

In addition to $V_{N}$, we choose a finite-dimensional subspace $V_{N}^{1}\subset H_{0}^{1}(T_{0}(\Omega))$ in the case of Dirichlet boundary conditions (resp. $V_{N}^{1}\subset H^{1}(T_{0}(\Omega))$ in the case of Neumann boundary conditions, cf. the discussion on the appropriate spaces in the preceeding section). Let $\varphi_{1}^{1},\ldots,\varphi_{N}^{1}$ be a basis of $V_{N}^{1}$. In order to approximate

[TABLE]

we choose a set $\{x_{1}^{1},\ldots,x_{N}^{1}\}$ of sample points in $T_{0}(\Omega)$ and require that

[TABLE]

for $j,m=1,\ldots,N$. If the $\varphi_{j}^{1}$ are a nodal basis with respect to the sample points $x_{1}^{1},\ldots,x_{N}^{1}$, then $\varphi_{k}^{1}(x_{m}^{1})=\delta_{km}$ and thus $\alpha_{jm}=T_{0,*}\varphi_{j}(x_{m}^{1})=\varphi_{j}(T_{0,*}^{-1}(x_{m}^{1})).$ In particular, if the sample points $x_{m}^{1}$ are chosen as the image points $x_{m}^{1}=T_{0}(x_{m})$ of the sample points $x_{1},\ldots,x_{N}$ in $\Omega$ and the basis $\varphi_{1},\ldots,\varphi_{N}$ is a nodal basis with respect to the points $x_{1},\ldots,x_{N}$, then $\alpha_{jm}=\varphi_{j}(T_{0,*}^{-1}(x_{m}^{1}))=\varphi_{j}(x_{m})=\delta_{jm}$, i.e. the representation matrix ${\bm{\alpha}}:=(\alpha_{jm})_{jm}$ of $T_{0,*}$ with respect to the two nodal bases $\varphi_{1},\ldots,\varphi_{N}$ and $\varphi_{1}^{1},\ldots,\varphi_{N}^{1}$ is the identity matrix. This latter case is the “adaptive” TO method from [13], where here we are considering only two discrete time-instances.

Approximating the stiffness matrix $K$.

With $K^{0}:=\left(\int_{\Omega}\nabla\varphi_{j}\bullet\nabla\varphi_{k}\;d\ell\right)_{j,k}$ and $K^{1}:=\bm{\alpha}^{\top}K^{0}\bm{\alpha}$ we obtain

[TABLE]

as an approximation to the stiffness matrix in (28).

Approximating the stiffness response matrix $L$.

We next describe an alternative way to compute the matrix $L$ in (30) based on the explicit approximation of the transfer operater described above. We will use only function evaluations of $T_{0}$ and $\dot{T}_{0}$. First, we manipulate the expression for $L$.

Proposition 3**.**

For $f,g\in H^{1}(\Omega)$,

[TABLE]

Proof.

We recall from (11) that $\dot{A}_{0}=-\left((DT_{0})^{-1}D\dot{T}_{0}(DT_{0})^{-1}(DT_{0})^{-\top}\right)^{\text{sym}}$. Next, we compute

[TABLE]

Similarly, we have

[TABLE]

We thus obtain

[TABLE]

∎

The right-hand-side of the expression in Proposition 3 has two types of terms, we now discuss their approximation. Given a function $g\in H^{1}(\Omega)$ we approximate it in $V_{N}$ by

[TABLE]

where $g_{j}=g(x_{j})$, since we assume $\varphi_{1},\ldots,\varphi_{N}$ to be a nodal basis on the nodes $x_{1},\ldots,x_{N}$. In the case of the adaptive TO method, the approximation of the pushforward $T_{0,*}g$ is therefore given by $\sum_{j,k}g_{j}\alpha_{jk}\varphi_{k}^{1}=\sum_{k}g_{k}\varphi_{k}^{1}$ (because $\alpha_{jk}=\delta_{jk})$. Finally, following [13] we approximate

[TABLE]

in particular, we have $\nabla T_{0,*}\varphi_{j}\approx\nabla\varphi_{j}^{1}$. Now we discuss the approximation of the term $D(T_{0,*}\dot{T}_{0})$ in (3). We denote by $\dot{T}_{0,1},\dots,\dot{T}_{0,n}$ the component functions of $\dot{T}_{0}$. Correspondingly,

[TABLE]

Since each $\dot{T}_{0,i}$ is a scalar-valued function, we approximate $\nabla(T_{0,*}\dot{T}_{0,i})$ in exactly the same way as we approximated $\nabla T_{0,*}g$, namely, we write $\nabla(T_{0,*}\dot{T}_{0,i})\approx\sum_{s}w_{s}^{i}\nabla\varphi^{1}_{s}$, where $w_{s}^{i}=\dot{T}_{0,i}(x_{s})$. Thus

[TABLE]

where $W=(w_{s}^{k})_{ks}$ and $D\Phi=(\partial_{k}\varphi_{s}^{1})_{sk}$. We then obtain the approximation

[TABLE]

which, using Proposition 3, yields

[TABLE]

as approximations for the entries of the matrix $L$ in (30).

7 Experiments

The code for the following experiments is available in the FEMDL333Available at https://github.com/gaioguy/FEMDL package. Since we identify coherent sets as level sets of eigenfunctions, and are interested in the evolution of coherent sets, we will begin this section with a short note about the evolution of level sets.

7.1 Level-Set Evolution

We wish to describe the change of the level sets of $u_{\varepsilon}$ as we perturb the parameter $\varepsilon$. From the level-set method [31], we note the following. For $\varepsilon\in[-\varepsilon_{0},\varepsilon_{0}]$, let $\Gamma_{\varepsilon}=\{x\in\Omega:u_{\varepsilon}(x)=c\}$ be a differentiable family of nontrivial closed curves in $\Omega$. Because $u_{\varepsilon}$ varies as $\varepsilon$ increases from 0, so too do the curves $\Gamma_{\varepsilon}$. Define the function $s:\Gamma_{0}\to\mathbb{R}$ to be the instantaneous speed of the curve $\Gamma_{0}$ with respect to $\varepsilon$ in the direction normal to $\Gamma_{0}$. Then $s$ satisfies the level-set equation

[TABLE]

and

[TABLE]

Extending this formula to all of $\Omega$ we obtain a vector field

[TABLE]

which describes the instantaneous evolution of level sets of $u_{0}$. We will use the vector field $v_{\rm level}$ in the following experiments to visualise the evolution of coherent sets with small changes in $\varepsilon$. We are not directly concerned with the possibility that the level sets occasionally undergo topological bifurcations as $\varepsilon$ is varied, however, it is well known [31] that such bifurcations are seamlessly captured by smooth evolution of the $u_{\varepsilon}$ with $\varepsilon$.

7.2 The Standard Map

We start with the standard map on the flat 2-torus, given by

[TABLE]

The parameter $a+\varepsilon$ controls the nonlinearity of the map and we investigate how varying $\varepsilon$ from [math] affects coherent sets. For the computations, we use a Delaunay triangulation on a regular grid of $100\times 100$ points on the 2-torus and Gauss quadrature of degree 2 in order to approximate the integrals in the CG approach.

In Figure 1 (left), we show the eigenvector $u_{0}$ at the second eigenvalue $\lambda_{0}=-1.08$444Note that the first eigenfunction is constant. The experiment works equally well for higher eigenfunctions, cf. the code at https://github.com/gaioguy/FEMDL. of the dynamic Laplacian $\Delta^{D}_{0}$ for the nominal value $a=0.98$ (corresponding to $\varepsilon=0$), which identifies a coherent set in the center of the domain (red). Figure 1 (second and third from left) displays $\dot{u}_{0}$ and $u_{0}+\varepsilon\dot{u}_{0}$ for $\varepsilon=0.5$, which – even though only a linear extrapolation – is quite similar to the exact eigenvector $u_{\varepsilon}$ for $\varepsilon=0.5$ at the second eigenvalue $\lambda_{\varepsilon}=-1.23$ of $\Delta^{D}_{\varepsilon}$ shown in Figure 1 (right). In fact, using again Gauss quadrature of order 2, we obtain the relative $L^{2}$-error $\|u_{\varepsilon}-(u_{0}+\varepsilon\dot{u}_{0})\|/\|u_{0}\|\approx 0.03$. We also obtain $\dot{\lambda}_{0}=-0.23$, which results in the estimate $\lambda_{0}+\varepsilon\dot{\lambda}_{0}=-1.19$ for $\lambda_{\varepsilon}$, i.e. using $\dot{\lambda}_{0}$ we get an estimate of $\lambda_{\varepsilon}$ with, again, a relative error of $3\%$. Because the more negative $\lambda_{0}$ is, the less coherent the associated coherent sets, $\dot{\lambda}_{0}<0$ indicates a loss of coherence as $\varepsilon$ is increased. These numbers and figures are obtained from the CG approach. The results from the TO approach are similar and in fact visually indistinguishable, so we do not show them here. Note that this is an advantage for the TO approach since its computational effort is considerably lower and it only requires point evaluations of the flow map.

Figure 2 shows the velocity field for the level-set curves at $\varepsilon=0$ which describes how the coherent set boundaries move in the fixed frame at $t_{0}$ as $a$ is varied from its nominal value $0.98$. We also show the level set at the value $c=0.1447$ which was selected from a line search of $c\in[0,\max_{x}u_{0}(x)]$ that minimises the dynamic Cheeger value in (3) with $\Gamma_{c}=\{x\in\Omega:u_{0}(x)=c\}$. As the parameter is increased from $a=0.98$ to a larger value, the boundary of the coherent set moves according to the velocity field $v_{\rm level}$ shown in Figure 2. We further compare the level sets of $u_{0}$, $u_{\varepsilon}$ and $u_{0}+\varepsilon\dot{u}_{0}$, i.e. the prediction of $u_{\varepsilon}$ by the linear Taylor approximation at $\varepsilon=0.5$. The predicted level set is indistinguishable from the true level set. Note that we can obtain predictions for the perturbed level sets very cheaply by computing contours for $u_{0}+\varepsilon\dot{u}_{0}$, $\varepsilon\in[0,0.5]$.

7.3 The transitory double gyre

In our second experiment, we consider the transitory flow from [30]. This is a non-periodic time-variant Hamiltonian system with Hamiltonian $H=-\psi$, where $\psi$ is the stream function

[TABLE]

and $s(t)$ is the transition function

[TABLE]

On the square $\Omega=[0,1]^{2}$, the vector field initially (at $t_{0}=0$) exhibits two gyres (if considered as a steady flow), with centers at $(\frac{1}{4},\frac{1}{2})$ and $(\frac{3}{4},\frac{1}{2})$. At the terminal time $t_{1}=1$, the vector field exhibits these gyres rotated by $\pi/2$ (again if considered as a steady flow). In this experiment we view the flow time $t_{1}$ as the perturbation parameter, and analyse the effect on the coherent sets under the unsteady flow as the flow time is increased.

For the computations, we approximate the flow map and the solution of the variational equation by Matlab’s ode45 command, i.e. an explicit Runge Kutta-scheme with adaptive step size control. We construct a Delaunay triangulation of a regular $100\times 100$ grid of nodes on the square and use Gauss quadrature of degree 5 in order to approximate the integrals in the CG approach.

In Figure 3 (left) we show the eigenvector $u_{0}$ at the second eigenvalue $\lambda_{0}=-50.4$ (all numbers from the TO approach) of the dynamic Laplacian $\Delta_{0}^{D}$ for $t_{1}=0.6$ (corresponding to $\varepsilon=0$). This figure is consistent with earlier experiments using transfer operator methods [16, Figure 7(a)] and Ulam approximation of the dynamic Laplacian [11, Figure 8 (left)]. The eigenvector identifies two coherent sets displayed in red and blue. The two plots in the center of Figure 3 display $\dot{u}_{0}$ (center left) and $u_{0}+\varepsilon\dot{u}_{0}$ (center right) for $\varepsilon=0.2$, corresponding to a flow time $t_{1}+\varepsilon=0.6+0.2=0.8$. Even for this quite large extrapolation value, we obtain a result that is very similar to the exact eigenvector $u_{\varepsilon}$ at the second eigenvalue $\lambda_{\varepsilon}=-61.6$ shown in the very right column in Figure 3. The corresponding relative $L^{2}$-error is $\|u_{\varepsilon}-(u_{0}+\varepsilon\dot{u}_{0})\|/\|u_{0}\|\approx 0.1$. We further obtain $\dot{\lambda}_{0}=-50.4$, which results in the estimate $\lambda_{0}+\varepsilon\dot{\lambda}_{0}=-60.5$ for $\lambda_{\varepsilon}$, i.e. using $\dot{\lambda}_{0}$ we get an estimate of $\lambda_{\varepsilon}$ with a relative error of $2\%$. As expected, lengthening the flow time leads to a loss of coherence, indicated by the negative sign of $\dot{\lambda_{0}}$, and approximately quantified by the magnitude of $\dot{\lambda_{0}}$.

Note that while agreeing qualitatively, there are some visible differences in $\dot{u}_{0}$ as computed by the CG vs. the TO approach (top row vs. bottom row of Figure 3). The TO method, however, is using considerably less information on the dynamics than the CG approach: the only dynamical data used in TO are the images of the $10^{4}$ grid nodes at the final time and their derivatives with respect to $\varepsilon$. In the CG approach, we have to time-integrate the variational equation for each quadrature node in each element of the triangulation which here amounts to $\approx 10^{5}$ time integrations. If one chooses a correspondingly finer grid for TO – so that the numerical effort is comparable to CG – the figures and the prediction on how the coherent set changes become visually indistinguishable.

Using either of the latter approaches, we identify coherent sets in the time frame at $t_{0}$ by selecting the value $c=0.8412$ from a line search of $c\in[0,\max_{x}u_{0}(x)]$ that minimises the dynamic Cheeger value in (3) with $\Gamma_{c}=\{x\in\Omega:u_{0}(x)=c\}$. Figure 4 shows the velocity field $v_{\rm level}$ of the level-set curves at $\varepsilon=0$, which describes how the coherent set boundaries move in the fixed frame at $t_{0}$ as $t_{1}$ is extended beyond $t_{1}=0.6$. Our linear extrapolation again appears to be accurate, even for a macroscopic extension of the flow time, as the change in the level-set contour from $\varepsilon=0$ (solid line) to $\varepsilon=0.2$ is consistent with the prediction by the velocity field.

Acknowledgements

FA is supported by a UNSW PhD scholarship, GF is partially supported by an ARC Discovery Project. OJ acknowledges support by the DFG within the Priority Programme SPP 1881 Turbulent Superstructures. A visit of FA and GF to the TUM Department of Mathematics was supported by a Universities Australia / DAAD Joint Research Co-operation Scheme, and they thank TUM for hospitality. A visit of OJ to the UNSW School of Mathematics and Statistics was also supported by this scheme and he thanks UNSW for hospitality.