Computation and optimal perturbation of finite-time coherent sets for aperiodic flows without trajectory integration

Computation and optimal perturbation of finite-time coherent sets for aperiodic flows without trajectory integration

Gary Froyland    Péter Koltai    and Martin Plonka

Understanding the macroscopic behavior of dynamical systems is an important tool to unravel transport mechanisms in complex flows. A decomposition of the state space into coherent sets is a popular way to reveal this essential macroscopic evolution. To compute coherent sets from an aperiodic time-dependent dynamical system we consider the relevant transfer operators and their infinitesimal generators on an augmented space-time manifold. This space-time approach avoids trajectory integration, and creates a convenient linearization of the aperiodic evolution. This linearization can be further exploited to create a simple spectral optimization methodology for diminishing or enhancing coherence. We obtain explicit solutions for these optimization problems using Lagrange multipliers and illustrate this technique by increasing and decreasing mixing of spatial regions through small velocity field perturbations.

1 Introduction

Analyzing complicated flows through their transport and mixing behavior has been and still is attracting a great amount of attention [MMP84, RKLW90, Wig92, HP98, Are02, JW02, Wig05, SLM05, FP09, Thi12, FPG14, KK16, KR18, HKK18], both from geometric and probabilistic points of view. Non-autonomous time-aperiodic dynamics poses additional difficulties, especially the case of finite time, where asymptotic notions cannot be applied.

The current work has two main contributions.

  1. Extending the work from [FK17], that deals with the time-periodic case, to the situation of general aperiodic finite-time dynamics. We detail a method to compute finite-time coherent sets [Fro13] of aperiodic flows that does not require any time-integration of trajectories.

  2. A technique to find a small perturbation of the underlying aperiodic vector field in a prescribed ball in a space or subspace of vector fields, which optimally enhances or destroys the existing finite-time coherent sets. This extends optimization results in [FS17], which considered the time-periodic setup of [FK17], to aperiodic dynamics.

The key construction in [FK17] is the representation of a -periodically forced flow on phase space  as an autonomous flow on time-augmented phase space , where denotes the unit circle. On this time-expanded phase space, the time coordinate is simply advanced at a constant rate:

Finite-time coherent sets on the time interval , as introduced in [FSM10, Fro13], were extracted from singular functions of the transfer operator , where the transfer operator is the linear operator describing the evolution of distributions under the dynamics subject to a small random perturbation. The crucial observation is that singular modes of are eigenmodes of , where denotes the dual of , and that for area-preserving dynamics (corresponding to divergence-free velocity fields ) the dual is the transfer operator of the time-reversed dynamics, i.e., the (again, slightly stochastically perturbed) flow governed by the time-reflected velocity field . The dynamic interpretation of this operator-based characterization is that (finite-time) coherent are those sets that are to a large extent mapped back to themselves by a noisy forward-backward evolution of the dynamics. This operator-based framework not only gives a qualitative framework for coherence; the singular values of provide quantitative bounds for coherence [Fro13, FPG14]—namely the closer the singular value is to one, the less mixing occurs between the coherent set and its exterior under the noisy dynamics.

By concatenating the “forward-time” and “backward-time” velocity fields on time intervals and —see (9) below—we will construct a system on the augmented space  that mimics the forward-backward evolution of the dynamics. In this way, the eigenmodes of the (Fokker–Planck-) generator of this augmented system yield the eigenmodes of ; i.e., the singular modes of , and from these singular modes the desired finite-time coherent sets. Again, the corresponding eigenvalues of the generator can be used to give quantitative bounds on coherence/mixing, see Theorem 7. We note that a numerical approach to extract coherent sets from  by solving the Fokker–Planck equation has been described in [DJM16]. In contrast to [DJM16] we do not require time-integration over , which is especially advantageous once we consider the optimization of coherence and mixing.

Connecting the spectral properties of the generator of the augmented-space system with the (finite-time) dynamical properties of the original system is a generalization of the results from [FK17], where it has been done for time-periodic velocity fields on infinite time intervals. The interplay between the spectra of the dynamics in augmented space and the non-autonomous dynamics in the original space has strong connections to the correspondence between evolution semigroups [How74] and two-parameter semigroups, as elaborated in, e.g., [CL99]; see also [EN00, Section VI.9] for a general introduction. We mention that by a similar construction, spatio-temporal dynamical patterns were extracted in [GD17] by considering the generator of the Koopman operator (the adjoint of the transfer operators considered here) associated with the augmented-space dynamics.

There are several different ways to measure mixing and mixedness under (stochastically perturbed) dynamics, such as considering dispersion statistics or the change of variation in a concentration field; see e.g. [Pro99, LH04, TDG04, Thi08]. Multiscale norms of mixing measure how “oscillatory” a concentration field is [MMP05]; see [Thi12] for a review. The most widely used approaches to the problem of mixing optimization search for switching protocols between some fixed velocity fields in order to optimize some topological [BAS00] or other mixing measure [MMG07, CAG08, OBPG15]. Other strategies include the optimal distribution of concentration sources [TP08] and geometric dynamical systems techniques [Bal15]. An interesting theoretical result is that arbitrary mixedness under advection-diffusion can be achieved in finite time solely by sufficiently increasing the strength of the (otherwise fixed) advective flow [CKRZ08]. If there are no restrictions to the choice of the velocity field, one can choose the one that is optimally mixing the actual concentration at every time instance [LTD11]. We also note that a related problem to mixing enhancement arises in statistical mechanics [LNP13] where the convergence toward the stationary distribution should be accelerated, e.g., to increase the efficiency of sampling.

Instead of focusing on one fixed concentration field, we will bound the mixing characteristics of a flow in terms of the objects that most inhibit mixing: coherent sets. As we mentioned above, finite-time coherent sets are characterized by the singular vectors of the transfer operator, and, equivalently, by the eigenvectors of the generator of the augmented-space process, while the corresponding eigenvalue delivers an upper bound on transport between the coherent set and its exterior. Thus, we can quantitatively access the mixing behavior of a flow on finite time through the spectrum of the augmented generator , and can target these eigenvalues if we want to enhance or diminish mixing.

Given a default velocity field and non-autonomous perturbations from an admissible space  of divergence-free velocity fields, our approach considers “small” that change a targeted eigenvalue of (thus also the singular value of ) locally optimally. This procedure can be iterated to obtain a larger perturbation in a gradient-method fashion. Optimizing singular values of the transfer operator directly is difficult as it would necessarily involve the variation of the nonlinear dynamics under the velocity field . Instead, a linearized optimization of the eigenvalue of the generator  leads to a very simple optimization problem (42), which can be solved via a linear system of the same dimension as . Moreover, the theory holds for infinite-dimensional perturbation spaces as well, since Fréchet differentiability of the transfer operator and its spectrum with respect to perturbing velocity fields has been established in [KLP18]. It should be noted that additive perturbations in the sense might be non-physical or technically not realizable, however that is the topic of flow control [BMT01, KB07], and certainly beyond the scope of the current work.

This work is structured as follows. In section 2 we introduce the -function based formalism to study advection-diffusion systems by the Fokker–Planck equation and its evolution operator, the transfer operator, in forward and backward time. In section 3 we consider purely advective transport between a family of sets and its exterior in terms of fluxes through the boundary of the family; this is a geometric analogue of the operator-based considerations that follow. Section 4 introduces the new reflected augmented generator needed to handle aperiodic, finite-time driving. We state formal connections between the spectrum of the reflected augmented generator and the reflected transfer operator, and provide spectral-based bounds on the maximal possible coherence of sets in phase space under the aperiodic dynamics. Section 5 contains a numerical demonstration of the efficacy of our trajectory-free approach. Section 6 describes the formal setup of the optimization problem designed to manipulate the position of the dominant spectral values of the reflected augmented generator, culminating in an explicit expression for the optimal time-dependent local perturbation of the velocity field. Section 7 specializes the infinite-dimensional results of section 6 to the numerical setting via discretization and includes a variety of examples of coherence reduction and enhancement. We conclude in section 8.

2 Advective-diffusive dynamics

Let be compact with smooth boundary (at least and piecewise ). We consider the time interval and the dynamics


with , for some , and being a standard Wiener process in . The initial point is distributed according to some initial density . The evolution of the density of the governing equation (1) is given by the Fokker–Planck equation or Kolmogorov forward equation [LM13, Section 11.6].


where is the normal derivative on the boundary. Associated to (2) is an evolution operator that transports a density at time 0 to the solution density of (2) at time . The evolution operator is an integral operator with stochastic111Doubly stochastic if the flow is volume-preserving. kernel that satisfies [Fro13, Assumptions 1 and 2].

2.1 Construction of a forward-backward process

For simplicity of presentation we assume that the velocity field is divergence free for all . We note that the remaining arguments in this section may be carried through for general velocity fields. Denote by the canonical scalar product of a Hilbert space . Following [Fro13] in the volume-preserving setting222The operators and are replaced by normalised versions for nonzero divergence velocity fields; these are denoted by and in [Fro13]., coherent sets over the time interval are extracted from the eigenfunctions of corresponding to large eigenvalues, where is the -adjoint of , defined to be the unique linear operator satisfying

for all . The eigenvalues of (the singular values of ) are known to lie in the interval (cf. [Fro13, p.3]). The rationale behind the operator is that describes evolution in forward time, describes evolution under the time-reversed dynamics, and coherent sets are characterized exactly by the property that they are “stable” under a noisy forward-backward evolution of the dynamics.

The adjoint operator is the solution operator to the Kolmogorov backward equation [PS08]:


The operator maps a density at time backward in time according to (4) to produce a density at time . We may simplify (4) using volume preservation:


thus we may write (4) as


Reversing time in (4) to obtain an initial value problem we get


using , and the velocity field . Comparing (2) and (7) we see that the natural evolution of the adjoint problem (the Kolmogorov backward equation) corresponds to the forward problem (the Kolmogorov forward equation) of the time reversed dynamics.

We wish to construct the process over the time interval that corresponds to the operator . We view as evolution on the time interval and we therefore shift (7) by time units, defining to obtain a forward problem on :


We denote the solution operator of this problem as ().

Finally, we concatenate the two forward problems (2) and (8) to make a single process over . We mark objects that live on this extended interval with a hat . Define the velocity field


using the reflection map


The resulting velocity field exhibits discontinuities in and whenever it does not vanish there, but one-sided derivatives exist. In what follows, we will solve the Fokker–Planck equation


over the interval ; more precisely on with -continuous concatenation at . Let us summarize the above construction with the following proposition.

Proposition 1.

The concatenation with comprises initializing (2) at time 0, solving forward using the vector field until time , then continuing to evolve (8) for another time units, but now using the reflected and shifted vector field for corresponding to the reversed dynamics.

3 Cumulative flux from a reflected family of sets

Before proceeding with the operator-based description of finite-time coherence, in this section we analyze the reflected dynamics by its flux through the boundary of a moving (possibly coherent) set. Our intention behind connecting this to a flux in augmented space (i.e., space-time) in Proposition 3 is partially to set the stage for the augmented-space operator-based description in section 4.

We consider a family of -dimensional sets , satisfying the following assumptions (the boundaries are piecewise smooth in space and differentiable in time):

Assumption 2.
  1. There exists a co-dimension 1 parameterisation set such that for each there is a bijective function with being piece-wise smooth, ( in and in ) piecewise.

  2. is well defined for all and all , where .

The cumulative outflow flux under the vector field from the family of sets is given by


where denotes the positive part, is the dimensional surface measure and is the outer normal unit vector.

In augmented state space we consider the augmented set


and define the instantaneous outflow flux from the set by


with denoting the -dimensional surface measure and denoting the augmented velocity field defined by

The result [FK17, Theorem 2] shows that (12)=(14) holds.

We extend this result to the reflected velocity field generated by a general aperiodic velocity field. We define a reflected family of sets in synchrony with the reflected vector field :


In particular, recalling from (10), for every the boundary has a parametrization and


At the right- and the left-sided partial derivatives of with respect to exists but they may not be equal. The family of normal vectors is mirrored in time (see Figure 1): for .

Figure 1: Illustration of the augmented reflected set and normal vectors in the case where .

We define the augmented reflected set

see Figure 1.

Proposition 3.

The cumulative outflow flux from the family of sets under the vector field , is equal to the cumulative absolute flux in and out of the family of sets , , under the vector field ; that is,


Furthermore, the time-integrated flux (17) is equal to the instantaneous absolute flux in augmented space:


Let us first prove the first equality. Therefore we do not yet need objects of the augmented setting. We split the integral


and see that the first integral looks promising. So we only need to treat the second integral (19). We use Fubini’s theorem and substitute using with being the Gram determinant. (Transformation theorem would not work, we need the sign we get from substitution.)

Combining the two calculations above we get the desired result. The second equality involving the objects of the augmented setting follows analogously to [FK17, Theorem 2]. ∎

4 Coherent families of sets and the generator on augmented phase space

In this section we create a spectral mapping theorem for our reflected augmented process (Proposition 4) and derive a bound for the finite-time coherence of a family of sets in terms of the second eigenvalue of a generator (the infinitesimal operator) of our augmented reflected advection-diffusion process (Theorem 7). To do this we build on discrete-time theory from [Fro13] with the periodic continuous-time theory from [FK17].

4.1 The evolution operator for the reflected process

It is well known that is a compact, integral preserving, real and positive operator333See [FK17, Lemma 30, at the end of the proof] for compactness, and [LM13] or [KKS16] for the other properties.. Furthermore is a self-adjoint operator with simple largest eigenvalue . Following [Fro13] one has the second eigenvalue satisfies


where in the volume-preserving444In the nonzero divergence case, is a reference measure describing the initial mass distribution of the (possibly compressible) fluid being evolved, and is the forward evolution of under a normalised version of denoted by in [Fro13]. setting and are both simply the Lebesgue measure. We now consider the problem (11) introduced in section 2 as a time-periodic problem on (we extend periodically). Following the considerations of section 2 the evolution operator from time to time is given by


The situation when is exactly is of particular importance:


Note that is self-adjoint for , .

4.2 The time-augmented generator

We now turn to the augmented reflected system


on , and note that is a standard Wiener process on ; in particular it is not constructed by time reflection. We define augmented versions of , , and , denoting them with bold symbols:

and is a dimensional standard Wiener process. The augmented system is

Considering the time-periodic problem with on we formulate a Fokker–Planck equation in augmented space:


Note there is no diffusion in the direction, as per the definition of . Let denote the augmented generator associated with the augmented Fokker–Planck equation (24) in augmented space:


with continuous concatenation conditions for the other parts of the boundary, . We may also write (24) and (25) in terms of the non-autonomous (“unaugmented”) dynamics:


where is the right-hand side operator of (11) at time , i.e., the generator of the reflected Fokker–Planck equation. Analogously to [FK17, Lemma 22] the following result holds.

Proposition 4.

Let be an eigenfunction of corresponding to the eigenvalue . One has then


for all and .


We will prove (27) using some results and ideas from [FK17]. First we apply [FK17, Lemma 30] piecewise on and concerning well-posedness and regularity. Therefore we consider the original problem (2) and the reflected, shifted, time-reversed problem (11). Now [FK17, Lemma 30] guaranteees the unique existence of a function with

and the regularity

Now we can proceed as in [FK17, Lemma 22]. Let and with . According to the construction above (26) we know

and this implies for all


Now is the evolution operator to the evolution equation

Therefore the function solves (28) uniquely. This gives the claim. ∎

Let be an eigenvalue of with an eigenfunction . Inserting and into Proposition 4 yields


Recalling that is a compact self-adjoint positive operator, it must be that for some . This implies


from which it follows that for some .

4.3 Coherent families of sets

In the specific case where the velocity field is periodic in time, [FK17] shows that the families of sets

have an escape rate (see [FK17, Definition 8]) of at most , where is the first nontrivial eigenvalue of corresponding to the eigenfunction . Because we consider the dynamics on a finite time interval, this notion of escape rate is replaced by the concept of a coherence ratio [Fro13].

In the general setting of aperiodic we will quantify the coherence of families and provide a construction of highly coherent families with associated rigorous coherence bound. To ensure the existence of our coherence quantifier, we place an assumption on the family of sets555Such a family of sets was called “sufficiently nice” in [FK17]. Here we simply call it “nice”..

Definition 5 (Nice family of sets).

A nice family of sets  is such that

  1. for every  the set is closed; and

  2. the function —where denotes the Hausdorff distance—is either left- or right-continuous at for every .

Definition 6 (Coherence ratio).

Let  be a family of measurable sets. Denote by  the law of the process  generated by the SDE (1) initialised with , where denotes a possibly signed Borel measure. For we define the coherence ratio of the family as


It was shown in [FK17, Appendix A.6] that for a nice family of sets the quantity (31) is well defined. If no initial distribution  is specified,  denotes the coherence ratio with respect to the uniform distribution (normalised Lebesgue measure on ).

The following theorem shows that the probability of a trajectory remaining in a family of sets constructed from the positive and negative parts of eigenfunctions of decays no faster than the rate given by the corresponding eigenvalues.

Theorem 7.

Let  with , and assume that the family  of sets with  is nice. If  is scaled such that , then


By Proposition 4 we have that is an eigenfunction of the integral preserving operator  at the eigenvalue . Thus we have , and with it again by the integral-preserving property. With the given scaling of  we have that , and let . Note that by Proposition 4 we have that , thus [FK17, Equation (A.10)] gives

Noting that , the claim follows. ∎

Theorem 7 can be considered a generalisation from the periodic velocity field setting of [FK17, Theorem 19] to the case of general aperiodic velocity fields.

5 Computational aspects

5.1 Numerical discretisation

We use the “Ulam’s discretisation for the generator” approach developed in [FJK13] for autonomous flows and extended in [FK17, Sections 7.2 and 7.3] for nonautonomous flows. In brief, referring to the above papers for further details, the Ulam discretisation for the generator yields a matrix that may be interpreted as a rate matrix of a finite-state, continuous-time Markov chain, with the states corresponding to a partition of into hypercubes (hyperrectangles) in . The entries which correspond to rates between boxes adjacent in the temporal coordinate direction (in which the time evolution is a rigid rotation of constant velocity 1) are given by , where  is discretised into intervals of length . The entries in the remaining space directions are computed from the rate of flux out of the hypercube faces by numerical integration of the component of the velocity field normal to (and pointing out of) the face; see the expression for in [FK17, Section 7.2]. The entries of the rate matrix corresponding to diffusive dynamics (in the space coordinates only, there is no diffusion in the time coordinate) are computed from a finite-difference approximation of the Laplace operator; see the expression for in [FK17, Section 7.2]. We then set . The matrix can also be interpreted as a rate matrix for a finite-state Markov chain; it has an eigenvalue and its spectrum is confined to the left half of the complex plane.

The reflected velocity field from (9) may be substituted for the velocity field used in [FK17] and the methodology of [FK17] employed; this is the approach taken in the numerical experiments below.

Remark 8.

The computation of uses only the outward-pointing velocity field values on the faces of the partition elements—similarly to how the outward flux is defined through the positive part of the inner product in (14)—, discarding the inward-pointing parts. Because of the reflected structure of , a slightly more efficient implementation would be to store the evaluations of the velocity field normal to hypercube faces in both directions (not only in the outward-pointing direction). The outward-pointing components would be used on the time interval , while the inward-pointing components would be used on the time interval , where they are outward-pointing because of the sign flip in (9)—similarly as it happens in the proof of Proposition 3. This would reduce by half the computational effort in evaluating the velocity field components normal to the hypercube faces. As we report below, the assembly of the generator matrices is relatively fast anyway, and we have not tried to optimise our implementation of Ulam’s method for the generator in this reflected setting.

5.2 Example: Bickley jet

We now apply the reflected augmented generator approach to a perturbed Bickley Jet [RBBV07]. The following model describes an idealized zonal jet in a band around a fixed latitude, assuming incompressibility, on which two traveling Rossby waves are superimposed. The velocity field is induced by the stream function

The constants are chosen according to [RBBV07]. The length unit is () and the time unit is days. For the amplitudes and the speed of the Rossby waves we choose

with being the Earth’s radius. Further we choose

The state space is periodic in direction and is given by . For good numerical tractability, we resolve our reflected space-time manifold with a spatially somewhat coarse grid, that is uniform in space () leads to square boxes needed for isotropic diffusion) and sufficiently finely resolved in time (). We choose .

The system described above is equipped with outflow boundary conditions instead of homogeneous Neumann conditions on . This leads to a slightly different spectral structure of the generator, which now generates a semigroup of substochastic operators. Its leading eigenvalue is strictly less than zero.

This Bickley Jet differs from the one investigated in [FK17, Section 7.6], as it is not periodic in time. Further, we expect the eigenvectors computed in [FK17, Section 7.6] to induce sharper and smoother separations, because they utilize the periodicity and describe asymptotic properties. The singular vectors computed here give rougher results, because they only have finite-time information available and only describe finite-time dynamics.

Analogously to [FK17, Section 7] our time-augmentation produces companion eigenvalues. Companion eigenmodes denote eigenmodes that are “higher order harmonics” of existing eigenmodes differing only in temporal modulation, and encoding the same coherence information; see below. For more details on the companion eigenvalues for the Ulam-discretisation we refer to [FK17, Section 7.3]. We will use and verify the relations derived there. Therefore we calculate the eigenvalues and vectors of with the smallest magnitude instead of largest real part using eigs(G,10,’SM’).

Table 1: Eigenvalues () of ordered in ascending magnitude and corresponding approximate singular values () of according to (30). The shaded eigenvalues correspond to companion modes. They do not yield purely real singular values through the exponentiation (30), because the numerically computed companions (33) contain a bias induced by discretization; see [FK17, Section 7.3] for further details.

Table 1 shows a gap after the first and sixth eigenvalue. Let us first discuss the leading eigenvectors. Figure 2 shows the eigenvectors corresponding to the dominant (i.e., smallest real part) eigenvalues. The first eigenvector, the quasistationary (or conditionally invariant) distribution highlights (in blue, see Figure 2(a)) parts of the domain that get pushed out of the region of consideration. The red regions are those parts of phase space that remain longest in under the diffusive dynamics (1). This effect is due to the outflow conditions. Note that this example is in this sense explorative, as our theory in the previous sections was only considering reflecting and not outflow boundary conditions.

The second eigenvector indicates an upper/lower separation. The other four eigenvectors show combinations of coherent vortices.

\thesubsubfigure First eigenvector (initial time).
\thesubsubfigure Second eigenvector (initial time).
\thesubsubfigure Third eigenvector (initial time).
\thesubsubfigure Fourth eigenvector (initial time).
\thesubsubfigure Fifth eigenvector (initial time).
\thesubsubfigure Sixth eigenvector (initial time).
Figure 2: Approximate leading eigenvectors of —that are according to (29) singular vectors of —computed from the Ulam discretization of the reflected augmented generator with a time-space resolution.

To investigate which elements of the spectrum of contain genuinely new dynamical information, we checked that the complex eigenvalues of  from to are all companion eigenvalues equal to (recall from section 5.1 that is the temporal grid spacing)


for an eigenvalue of