Computing Coherent Sets using the Fokker-Planck Equation
CComputing Coherent Setsusing the Fokker-Planck Equation
Andreas Denner Oliver Junge Daniel Matthes ∗ October 17, 2016
Abstract
We perform a numerical approximation of coherent sets in finite-dimensionalsmooth dynamical systems by computing singular vectors of the transfer opera-tor for a stochastically perturbed flow. This operator is obtained by solution ofa discretized Fokker-Planck equation. For numerical implementation, we employspectral collocation methods and an exponential time differentiation scheme. Weexperimentally compare our approach with the more classical method by Ulam thatis based on discretization of the transfer operator of the unperturbed flow.
Many fluid flows at the onset of turbulence exhibit regions which disperse slowly with time(so called coherent sets [10, 13]) while other regions disaggregate comparatively quickly.Often, the boundary of a slowly dispersing region can be associated to a lower-dimensionalobject (a so called
Lagrangian coherent structure [17, 15, 16]) which serves as a transportbarrier for Lagrangian particles within a coherent region [11].Based on the related concept of almost invariant (resp. metastable ) sets in time-invariant dynamical systems [5, 8], recently a framework has been proposed for the com-putation of coherent sets via singular vectors of a certain smoothed transfer operatorwhich describes the evolution of probability densities on phase space under the given dy-namics [7]. Numerically, this operator can be approximated by a Galerkin [22, 5] methodvia evaluating the flow map explicitly by time integration of the vector field. Dependingon the type of approximation space chosen, a diffusion operator has either to be appliedexplicitly or is already implicitly present in the discretization ( via numerical diffusion).In this manuscript, we propose to compute coherent sets by directly solving the Fokker-Planck equation instead. More precisely, instead of computing the evolution of the basisof our approximation space under the deterministic dynamics and then applying diffusion,we directly compute the evolution of this basis under the stochastic push forward operatorgiven by the solution operator of the Fokker-Planck equation. This advection-diffusionequation can efficiently be discretised using spectral collocation (cf. also [9]). In order todeal with aliasing in the case of dominating advection we use a skew symmetric form ofthe advection term and in order to deal with stiffness in time due to the Laplace operatorwe employ an exponential time differentiation (etd) integrator. ∗ Center for Mathematics, Technische Universit¨at M¨unchen, 85747 Garching bei M¨unchen a r X i v : . [ m a t h . D S ] O c t s a key advantage of our method we only need to sample the vector field at eachtime instance on a fixed grid of rather coarse resolution. In particular, we do not need tointegrate trajectories of (Lagrangian) particles and no interpolation of the vector field topoints off the grid. We consider a time-dependent ordinary differential equation ˙ x = b ( t, x ), b : R × X → R d ,on some bounded domain X ⊂ R d . We fix some initial and final time t , t ∈ R andassume that the vector field b is continuous and locally Lipschitz w.r.t. x for all t ∈ [ t , t ],such that the associated flow map Φ = Φ( · , t , t ) : X → X is uniquely defined. In ordernot to obscure the key ideas we restrict to the case of X being a hyperrectangle andfurthermore b being periodic in x and divergence free, i.e. the flow map Φ being volumepreserving.Roughly speaking, we would like to compute a set A ⊂ X which disperses slowlyunder the evolution of Φ, i.e. which roughly retains its shape while being moved aroundby the flow Φ. A such set A will be called coherent [13].Inspired by the associated notion of an almost invariant (resp. metastable ) set in theautonomous setting [5], we start formalizing this request by asking for a pair of sets A , A ⊂ X such that A will approximately be carried to A by Φ in the sense that κ ( A , A ) := m ( A ∩ Φ − A ) m ( A ) ≈ , (1)where m denotes Lebesgue (i.e. volume) measure. Evidently, with A = Φ( A ) we obtain κ ( A , A ) = 1 for any A ⊂ X , so this is not a well defined problem yet. In fact, (1)does not impose any condition on the geometries of the sets A and A . In particular, theimage set A = Φ( A ) might be stretched and folded all over the domain X – but this isnot the type of coherent set we have in mind. Froyland [7] observed that κ ( A , Φ( A ))is not close to 1 for every set A any more as soon as one artificially adds some randomperturbation to the dynamics (cf. [5, 4] for related ideas in the autonomous context). More formally, Froyland proposes the following approach [7], see also [12]. We will employthe transfer operator (resp.
Frobenius-Perron operator or push forward ) P associated toΦ. This operator describes how densities u : X → [0 , ∞ ) are evolved by Φ: If a set ofpoints in X is initially distributed according to some density u then after flowing thesepoints forward with Φ they will be distributed according to P u . In our case of Φ beingvolume preserving, P : L = L ( X, m ) → L is given by P u = u ◦ Φ − . In the sequel, we restrict our attention to L = L ( X, m ) ⊂ L as we want to use thecanonical scalar product on L . We are going to add an ε -small random perturbation tothe flow map Φ via complementing the action of P by a diffusion operator D ε : L → L , D ε f ( y ) = (cid:90) X α ε ( y − x ) f ( x ) dx (2)2 A t D ε A t T PD ε A t D ε D ε P ε A t := D ε PD ε A t ε Figure 1: Graphical illustration of the action of the perturbed transfer operator P ε as thecomposition of diffusion operators at initial and final time and the transfer operator P .where α ε : X → [0 , ∞ ) is a bounded kernel with (cid:82) X α ε ( x ) dx = 1 and α ε → δ in adistributional sense as ε →
0. Here, we use a diffusion with bounded support of radius ε ,namely α ε ( x ) = 1 B ε ( x ) ∩ X /m ( B ε ( x ) ∩ X ). With P and D ε , we finally define the evolutionoperator [7] P ε : L → L , P ε := D ε PD ε , cf. Figure 1. Note that P ε is stochastic , i.e. positive and P ε X = 1 X , where 1 X denotesthe characteristic function on X , since we assume Φ to be volume preserving. Further,with this choice of α ε , P ε is compact and as the vector field b is divergence-free has asimple leading singular value σ = 1 [5, 7].In order to compute coherent sets, note that if κ ( A , A ) ≈
1, then κ ( A c , A c ) ≈ m ( A ∩ Φ − A ) = (cid:104) A , A ◦ Φ (cid:105) = (cid:104)P A , A (cid:105) , where (cid:104)· , ·(cid:105) is thestandard inner product on L . We therefore define ρ ε ( A , A ) := (cid:104)P ε A , A (cid:105) m ( A ) + (cid:104)P ε A c , A c (cid:105) m ( A c ) (3)and look for some pair ( A , A ) of sets which maximizes this quantity. We getmax A ,A ρ ε ( A , A ) − A ,A (cid:104)P ε ( c A − c − A c ) , c A − c − A c (cid:105) =: ( (cid:63) ) , where c = (cid:112) m ( A c ) /m ( A ) and c = (cid:112) m ( A c ) /m ( A ) and since the functions f := c A − c − A c and g := c A − c − A c statisfy (cid:107) f (cid:107) = (cid:107) g (cid:107) = 1 and (cid:104) f, X (cid:105) = (cid:104) g, X (cid:105) = 0,we obtain, by relaxing to arbitrary functions f, g with zero mean,( (cid:63) ) ≤ max (cid:107) f (cid:107) = (cid:107) g (cid:107) =1 {(cid:104)P ε f, g (cid:105) : (cid:104) f, X (cid:105) = (cid:104) g, X (cid:105) = 0 } . (4)This problem is much easier to solve than ( (cid:63) ), where we need to maximize over charac-teristic functions. As for fixed f ∈ L max (cid:107) g (cid:107) =1 (cid:104)P ε f, g (cid:105) = (cid:28) P ε f, P ε f (cid:107)P ε f (cid:107) (cid:29) = (cid:107)P ε f (cid:107) and (cid:104)P ε f, X (cid:105) = 0, we obtain( (cid:63) ) ≤ max (cid:107) f (cid:107) =1 {(cid:107)P ε f (cid:107) : (cid:104) f, X (cid:105) = 0 } = (cid:107)P ε (cid:107) L ( V,m ) , where V = { f ∈ L : (cid:104) f, X (cid:105) = 0 } is the orthogonal complement of span 1 X .3 A t P ε A t P ε Figure 2: In contrast to the approach using diffusion at initial and final time cf. 1, in ourapproach diffusion is built into the model of the dynamical system.This operator norm is given by σ ( P ε ) [7], the second largest singular value of P ε and the maximizing function is f = v , the associated right singular function. As v isan approximation to c A − c − A c the common heuristics is to identify an associatedcoherent set by A := { x ∈ X : v ( x ) > θ } where θ ∈ R is some appropriately chosen threshold [5], [8]. Correspondingly, u = P ε v is the associated left singular function and A = { x ∈ X : u ( x ) > } .To sum up, we can compute a pair A , A of coherent sets via computing the secondsingular value and its corresponding singular vectors of a slightly perturbed Frobenius-Perron operator. Often a low order spatial discretisation such as Ulam’s method is used[26], which is a Galerkin projection on indicator functions of a partition of X into boxes X i . Usually P is then computed numerically via sampling K test points x i in every box X i , compute Φ( x i ) and count how many fall in B j : P ji = x i ) ∈ X j K .
The method automatically adds sufficient numerical diffusion so that we can actuallydirectly employ the unperturbed operator P . However, with increasing resolution thisnumerical diffusion decreases which results in all singular values approaching the valueone, cf. [19, 7]. Instead of explicitly applying diffusion at the beginning and the end of the time intervalunder consideration, we propose to incorporate a small random perturbation continuouslyin time, i.e. instead of considering a deterministic differential equation we now use thestochastic differential equation dx = b ( t, x ) dt + εdB (5)in order to define the flow map Φ. Here, ( B t ) t ≥ is d -dimensional Brownian motion and ε >
0. Since we assume b ( t, · ) to be Lipschitz, X to be bounded and b to be periodic in x , for any initial condition ξ ∈ X , (5) has a unique continuous solution x in the sense of[24], Thm. 5.2.1.The transfer operator P ε associated to this stochastic differential equation is given bythe solution operator of the parabolic Fokker-Planck equation ∂ t u = L ε u := ε ∆ u − div( ub ) . (6)4ppropriate boundary conditions are chosen (e.g., periodic or homogeneous Neumannboundary conditions), so that for all u, w ∈ L ( X ) in the domain of L ε holds: (cid:104) w, L ε u (cid:105) = − (cid:90) X ∇ u · ∇ w + (cid:90) X u b · ∇ w. (7)More precisely, P ε u = u ( t , · ), where u is the solution to (6) with initial condition u ( t , · ) = u . Lemma 1. If (cid:107) b (cid:107) C = sup s ∈ [ t ,t ] sup x ∈ X max {| b ( s, x ) | , | ∂ x b ( s, x ) | , . . . , | ∂ x d b ( s, x ) |} < ∞ , then P ε : L ( X ) → L ( X ) is compact. A proof of this claim is given in the appendix. Note that with Schauder’s theoremalso P ε, ∗ is compact. Since we assumed b to be divergence free, we also have: Lemma 2. P ε : L → L and P ε, ∗ : L → L are stochastic.Proof. We set ε = √ L ε X = ∆1 X − div(1 X b ) = div b =0, it follows that 1 X is a steady state of (6), and consequently P ε X = 1 X .For the adjoint operator P ε, ∗ we test with u ∈ L ( X ): (cid:104)P ε, ∗ X , u (cid:105) = (cid:104) X , P ε u (cid:105) = (cid:90) X P ε u = (cid:90) X u = (cid:90) X u X = (cid:104) X , u (cid:105) , where we have used that P ε is integral conserving: ∂ t (cid:90) X u = (cid:90) X L ε u = (cid:104) X , L ε u (cid:105) = − (cid:90) X ∇ u · ∇ X + (cid:90) X u b · ∇ X = 0 , thanks to the integration-by-parts rule (7). Hence (cid:82) X u = (cid:82) X P ε u for all u ∈ L ( X ), cf.[21]. By the Riesz representation theorem, we can conclude that also P ε, ∗ X = 1 X .For proving positivity of P ε , we consider the evolution of the negative part u − ( s, x ) = − min( u ( s, x ) , ∂ t ( u − ) = − u − ∂ t u almost everywhere on X × ( t , t ), we have,using the integration-by-parts rule (7),12 ∂ t (cid:90) X u − = − (cid:90) X u − ∂ t u = − (cid:90) X u − L ε u = (cid:90) X ∇ u · ∇ u − − (cid:90) X u b · ∇ u − = − (cid:90) X |∇ u − | + 12 (cid:90) X b · ∇ ( u − )= − (cid:90) X |∇ u − | − (cid:90) X div( b ) u − = − (cid:90) X |∇ u − | ≤ . (8)Hence if u ( t , · ) is non-negative, u ( t, · ) is non-negative for all t > t , as the norm of itsnegative part does not increase. To show positivity of P ε, ∗ , let two non-negative functions u, w ∈ L ( X ) be given. Then (cid:104)P ε, ∗ w, u (cid:105) = (cid:104) w, P ε u (cid:105) ≥ P ε . For any fixed non-negative w , this relation holds for all non-negative u . This implies that P ε, ∗ w is non-negative. 5 emma 3. The leading singular value σ = 1 of P ε is isolated.Proof. With the same argument as in (8) we obtain for solutions to (6):12 ∂ t (cid:90) u ≤ − (cid:90) |∇ u | . This implies that (cid:107) u ( t, · ) (cid:107) is non-increasing in time. Moreover, if the initial datum u isnot constant on X , its gradient does not vanish almost everywhere, and hence (cid:107) u ( t, · ) (cid:107) decreases on a small time interval ( t , t + τ ). Thus (cid:107)P ε u (cid:107) < (cid:107) u (cid:107) , and u cannot bean eigenfunction of P ε corresponding to the eigenvalue 1. With the same argument theconstant function is the only eigenfunction of P ε, ∗ corresponding to the eigenvalue 1.Hence the leading singular value σ = 1, which is an eigenvalue of P ε, ∗ P ε is isolated.Hence P ε is compact and doubly stochastic with isolated, simple leading singular value σ . So we are in a similar setting as in section 3 and apply the same constructions. Remark 1.
These considerations also work for a non-conservative vector field b andcorresponding probability measure µ (cid:54) = m by suitably normalizing the operator P resp. P ε ,cf. [7]. In order to approximate the transfer operator P ε , we choose a finite dimensional approx-imation space V N ⊂ L ( X ) and use collocation. As the Fokker-Planck equation (6) isparabolic and since we assume the vector field b ( t, · ) to be smooth for all t , its solution u ( t, · ) is smooth for all t > t (see e.g. [6], Ch. 7, Thm. 7). To exploit this, we choose V N as the span of the Fourier basis ϕ k ( x ) = e i (cid:104) k ,x (cid:105) , k ∈ Z d , (cid:107) k (cid:107) ∞ ≤ ( N − / , N odd . Note that dim( V N ) = N d . Choosing a corresponding set { x , . . . , x M } ⊂ X of collocationpoints (typically on an equidistant grid), the entries of the matrix representation P ε of P ε are given by P εj k = P ε ϕ k ( x j ) , (9)where k ∈ Z d , (cid:107) k (cid:107) ∞ ≤ ( N − / j = 1 , . . . , M . For P ε , we then compute singularvalues and vectors via standard algorithms. Note that we might choose M ≥ N , i.e.more collocation points than basis functions. This turns out to be useful since we expectthe maximally coherent sets to be comparatively coarse structures which can be capturedwith a small number of basis functions. Solving the Fokker-Planck equation.
In order to compute P ε ϕ k in (9) for some basisfunction ϕ k ∈ V N we need to solve the Fokker-Planck equation (6) with initial condition u ( t , · ) = ϕ k . This can efficiently be done in Fourier space via integrating the Cauchyproblem ∂ t ˆ u = ε ˆ∆ˆ u − ˆdiv F ( F − (ˆ u ) b ) , ˆ u ( t , · ) = ˆ ϕ k , in time, where ˆ v = F ( v ) is the Fourier transform of v ∈ V N . Note that the differentialoperators in Fourier space reduce to multiplications with diagonal matrices, while F and F − can efficiently be computed by the (inverse) fast Fourier transform.6 liasing. One problem with this formulation is the possible occurence of aliasing. Asˆ u and ˆ b are trigonometric polynomials of degree N , the multiplication F − (ˆ u ) b in theadvection term leads to a polynomial F − (ˆ u ) b of degree 2 N that cannot be representedin our approximation space V N . The coefficients of degree ≥ N of this polynomial act onthe coefficients of lower degree, leading to unphysical contributions in these. This showsup in high oscillations and blow ups (see [2], Ch. 11) in the computed solution. One wayto deal with this problem is to use the advection termdiv ( bu ) = 12 div ( bu ) + 12 ( b ∇ u ) . The spectral discretization of the left-hand side is not skew symmetric, but the discretiza-tion of the right-hand side is [27]. This leads to purely imaginary eigenvalues of theresulting discretization matrix and hence to mass conservation. Consequently, for theunperturbed operator ( ε = 0), the resulting matrix has eigenvalues on the unit circle.However, this approach has to be used carefully as even though the solution does notblow up, it might still come with a large error, e.g. small scale structures may be sup-pressed. If the system produces such small scale structures the grid has to be chosen fineenough to resolve them. Time integration.
For low resolutions, the time integration of the space discretizedsystem can be performed by a standard explicit scheme. For higher resolutions, thestiffness of the system due to the Laplacian becomes problematic and a more sophisticatedmethod must be employed. Here, we use the exponential time differentiation scheme [3] forthe space discretized system. The etd-scheme separates the diffusion term L = ε D , where D is the discretized Laplacian, from the advection term N ( u, t ) = − div F ( F − (ˆ u ) b ( · , t )),where b is evaluated via spectral collocation. The system can hence be written as u t = L u + N ( u, t ) . (10)Via multiplying (10) with e − t L and integrating from t to t we obtain u ( t ) = e h L u ( t ) + e h L (cid:90) h e − τ L N ( u ( t + τ ) , t + τ ) dτ, (11)with h = t − t . A numerical scheme is derived by approximating the integral in (11), e.g.by a Runge-Kutta 4 type rule, resulting in scheme called etdrk4 . To this end note that D isa diagonal matrix. We use the version in [20], which elegantly treats a cancellation problemoccurring in a naive formulation of etdrk4 by means of a contour integral approximatedby the trapezoidal rule. Extraction of coherent sets
For the extraction of coherent sets several methods canbe used. For the decomposition into exactly two coherent sets a simple thresholding oran a posteriori line search can be used [7, 8], for the decompositition of the domain into n coherent sets the first n singular vectors should be considered (see [5, 25, 18] for theautonomous case) and post processed via a simple clustering heuristic, e.g. k-means, see[1, 14]. We here focus on the computation of singular vectors.7igure 3: Quadruple gyre vector field at t = 0 and t = 10 .
25. The four gyres are separatedby a horizontal and a vertical line, such that their intersection point moves on the diagonal.
The first numerical example is a two dimensional flow (cf. Fig. 3), an extension of thewell known double gyre flow, given by˙ x = − g ( t, x, y )˙ y = g ( t, y, x )on the 2-torus [0 , × [0 , g ( t, x, y ) = π sin( πf ( t, x )) cos( πf ( t, y )) ∂ x f ( t, y )and f ( t, x ) = δ sin( ωt ) x + (1 − δ sin( ωt )) x . We fix δ = 0 . ω = 2 π , t = 0, t = 10 . h = 0 .
205 (i.e. 50 time steps) and choose ε = 0 .
02 in such a way that the six largestsingular values of P ε roughly equal those obtained from Ulam’s method (without explicitdiffusion).We use M = 15 collocation points and N = 5 basis functions in each direction andcompute the first four singular values and right singular vectors of P ε . As shown in Fig. 4(top row), they nicely reveal the gyres in their sign structures. The computation of P ε takes less than a second, the computation of all singular values and -vectors less than 0 . . For comparison, in the bottom row of Fig. 4, we show the same singular vectorscomputed via Ulam’s method (without explicit diffusion) on a 32 ×
32 box grid using100 sample points per box. Here, the computation of the transition matrix takes around5 seconds, the computation of the six largest singular values resp. vectors less than 0.2seconds.
We now turn to a case where the vector field is only given on a discrete grid. Here, theapproach proposed in Section 4 is particularly appealing if we chose the grid points ascollocation points (resp. a subset of them): In contrast to methods which are based on anexplicit integration of individual trajectories (as Ulam’s method), no further interpolationof the vector field is necessary. Furthermore, depending on the initial point of a trajectory, Computation times are measured on an 2 . − . . σ = 0 . σ = 0 . σ = 0 . σ = 0 . − − . . σ = 0 . σ = 0 . σ = 0 . σ = 0 . X = [0 , π ] , ∂ v ∂t = ν ∆ v − ( v · ∇ ) v − ∇ p ∇ v = 0 , where v denotes the velocity field, p the pressure, and ν > w := ∇ × v , the Navier-Stokes equation in 2D can be rewritten as vorticityequation DwDt = ∂ t w + ( u ∇ ) w = ν ∆ w (12)∆ ψ = − w. where the pressure p cancels from the equation. We can extract the velocity field v fromthe streamline function ψ via v = ∂ y ψ and v = − ∂ x ψ .The equation can be integrated by standard methods, e.g. a pseudo spectral methodas proposed in [23]. For a first experiment, we choose an initial condition inducing threevortices, two with positive and one with negative spin, as initial condition: w (0 , x, y ) = e − (cid:107) ( x,y ) − ( π, π ) (cid:107) + e − (cid:107) ( x,y ) − ( π, − π ) (cid:107) − e − (cid:107) ( x,y ) − ( π , π ) (cid:107) . We solve (12) on a grid with 64 collocation points in both coordinate direction. For thecomputation of coherent sets we chose n = 16 basis functions and N = 32 collocationpoints in both directions, as well as t = 0 and t = 20. We hence use only every second9igure 5: Top row: vector field at t = 0 (left), second right singular value computed viaFokker-Planck (center) and via Ulam (right). Bottom row: vector field at t = 20 (left),second left singular vector computed via Fokker-Planck (center) and via Ulam (right).collocation point of the computed vector field. However, as the induced coherent struc-tures are way bigger than the grid size, this does not affect the result. The underresolutionof the vector field can be interpreted as additional diffusion. We use ε = 10 − which is ofthe same order as the grid resolution. In Fig. 5, we show the vector field at time t = 0(left) as well as the second right singular vector (center) in the first row as well as thevector field at t = 20 and the corresponding second left singular vector (center) in thesecond row. The computation took 35 seconds.For comparison, we show the same singular vectors computed via Ulam’s method(right) on a 32 ×
32 box grid using 100 sample points per box which were integrated byMatlab’s ode45 . Here, we need to interpolate the vector field between the grid pointsusing splines (i.e. using interp2 in Matlab). This computation also took 35 seconds.For a second experiment, we use a turbulent initial condition by choosing a real numberrandomly in [ − ,
1] from a uniform distribution at each collocation point. For the coherentset computation, we restrict the time domain to [ t , t ] = [20 ,
40] since then the initialvector field has smoothed somewhat, cf. Fig. 6. Here, we choose M = 64, i.e. morecollocation points than in the examples before as the vector field lives at smaller scales, N = 16 and ε = 10 − which is of the same order as ν . The results are non-obvious pairsof coherent sets, cf. Fig. 6. The maximally coherent set indicated by the second singularvectors describes the vortex in the upper left region. The computation took 100 seconds.For comparison, we show the same singular vectors computed via Ulam’s method ona 32 ×
32 box grid using 25 sample points per box. Here, we need to interpolate thevector field between the grid points using splines (i.e. using interp2 in Matlab). Sincethe vector field is turbulent, using Matlab’s ode45 for the vectorized system is infeasible.We therefore choose a fixed time-step of h = 0 .
01, such that the result does not seem tochange when further decreasing h . This computation also took roughly 100 seconds.10igure 6: Top row: vector field at t = 20 (left), second right singular vector computedvia Fokker-Planck (center left) and via Ulam (center right) and an incoherent set (right).Bottom row: vector field at t = 40 (left), second left singular vector computed via Fokker-Planck (center left) and via Ulam (center right) and the evolution of the incoherent set(right). Finally, based on the quadruple gyre, we construct the following flow in R with eightgyres, given by the equations ˙ x = g ( t, x, y ) − g ( t, x, z )˙ y = g ( t, y, y ) − g ( t, y, x )˙ z = g ( t, z, x ) − g ( t, z, y )on the 3-torus X = [0 , with g and the other parameters from Section 6.1. By con-struction the dynamics of this system exhibits eight gyres in each quadrant of the cube[0 , . The cross sections of this vector field at, e.g. { x = 0 . } and { x = 1 . } are againgiven by Fig. 3.In Fig. 7 we show the second to fourth left singular vectors, computed using 16 col-location points and 4 basis functions in each direction with ε = 0 .
1. Interestingly, noneof the obvious gyres is identified by the second and the third singular vectors, but twocoherent sets with ‘centers’ at [1 , ,
1] and [2 , , t = 10 .
25. The computation time is 30 seconds.
We proposed to compute transfer operators in time-variant flows fields by a direct inte-gration of the associated Fokker-Planck equation using spectral methods. In particular,this approach does not require to integrate Lagrangian trajectories, which is particularlybeneficial if the underlying flow field is only given by (a grid of) data. However, the11igure 7: Octuple gyre: Second to fourth right singular vectors at time t = 0 (redpositive, blue negative level set).Figure 8: Octuple gyre: Second to fourth left singular vectors at time t = 10 . Acknowledgement
The authors acknowledge the support by the DFG Collaborative Research Center SF-B/TRR 109 ”Discretization in Geometry and Dynamics”. A.D. acknowledges the supportby the ”International Helmholtz Graduate School for Plasma Physics (HEPP)”.
AppendixA Proof of Lemma 1
Proof.
We prove the lemma in the slightly more general case that div( b ) (cid:54) = 0. (cid:107) b (cid:107) C :=sup s ∈ [ t ,t ] , x ∈ X {| b ( t, x ) | < ∞ , i = 1 , . . . , d } . Let u = u ( t, x ) be the solution of (6) withinitial condition f = u ∈ L ( X ). Without loss of generality we set ε = √ t = 0.12e first note that (cid:107) u (cid:107) is bounded by (cid:107) u (cid:107) :12 ddt (cid:107) u (cid:107) = (cid:104) u, ∂ t u (cid:105) = (cid:104) u, ∆ u + div ( ub ) (cid:105) = (cid:104) u, ∆ u (cid:105) + (cid:104) u, div ( ub ) (cid:105) = −(cid:107)∇ u (cid:107) − (cid:104)∇ u, ub (cid:105) = −(cid:107)∇ u (cid:107) − (cid:104) ∇ ( u ) , b (cid:105) = −(cid:107)∇ u (cid:107) − (cid:104) u , div ( b ) (cid:105) ≤ −(cid:107)∇ u (cid:107) + 12 (cid:107) u (cid:107) (cid:107) b (cid:107) C . (13)Gronwall’s inequality thus implies that for all t > (cid:107) u (cid:107) ≤ e t (cid:107) b (cid:107) C (cid:107) u (cid:107) . (14)We now show that also (cid:107)∇ u (cid:107) is bounded by (cid:107) u (cid:107) . Because of (13) we have (cid:107)∇ u (cid:107) ≤ − ddt (cid:107) u (cid:107) + 12 (cid:107) u (cid:107) (cid:107) b (cid:107) C . Integrating from t = t = 0 to t = t we obtain (cid:90) t (cid:107)∇ u (cid:107) dt ≤ (cid:0) (cid:107) u (cid:107) − (cid:107) u ( t ) (cid:107) (cid:1) + 12 (cid:107) b (cid:107) C (cid:90) t (cid:107) u (cid:107) dt (14) ≤ (cid:107) u (cid:107) + 12 (cid:107) b (cid:107) C (cid:90) t e t (cid:107) b (cid:107) C (cid:107) u (cid:107) ds = 12 e t (cid:107) b (cid:107) C (cid:107) u (cid:107) . (15)Therefore there is at least one t ∗ ∈ [0 , t ], such that (cid:107)∇ u ( t ∗ , · ) (cid:107) ≤ t e t (cid:107) b (cid:107) C (cid:107) u (cid:107) , (16)and we finally get12 ddt (cid:107)∇ u (cid:107) = 12 ddt (cid:104)∇ u, ∇ u (cid:105) = 12 (cid:104) ddt ∇ u, ∇ u (cid:105) + 12 (cid:104)∇ u, ddt ∇ u (cid:105) = −(cid:104) ∆ u, ∂ t u (cid:105) = −(cid:104) ∆ u, ∆ u + ∇ ( ub ) (cid:105) = −(cid:107) ∆ u (cid:107) − (cid:104) ∆ u, ∇ ub (cid:105) − (cid:104) ∆ u, u div ( b ) (cid:105) C − S ≤ −(cid:107) ∆ u (cid:107) + (cid:107) ∆ u (cid:107) (cid:107)∇ u (cid:107) (cid:107) b (cid:107) C + (cid:107) ∆ u (cid:107) (cid:107) u (cid:107) (cid:107) b (cid:107) C ≤ (cid:107)∇ u (cid:107) (cid:107) b (cid:107) C + 12 (cid:107) u (cid:107) (cid:107) b (cid:107) C , where we used that ( a c + d e ) ≥ − b + bac + bde for a, b, c, d, e ∈ R .Integration from t = t ∗ to t = t yields (cid:107)∇ u ( t ) (cid:107) ≤ (cid:107)∇ u ( t ∗ ) (cid:107) + (cid:90) t t ∗ (cid:107)∇ u (cid:107) (cid:107) b (cid:107) C dt + (cid:90) t t ∗ (cid:107) u (cid:107) (cid:107) b (cid:107) C dt (14)(15)(16) ≤ t e t (cid:107) b (cid:107) C (cid:107) u (cid:107) + (cid:107) b (cid:107) C e t (cid:107) b (cid:107) C (cid:107) u (cid:107) + (cid:107) b (cid:107) C (cid:90) t t ∗ e t (cid:107) b (cid:107) C (cid:107) u (cid:107) dt ≤ (cid:18) t + 12 (cid:107) b (cid:107) C + (cid:107) b (cid:107) C (cid:19) (cid:107) u (cid:107) e t (cid:107) b (cid:107) C . (17)13ombining equations (14) and (17) we obtain (cid:107)P ε u (cid:107) H = (cid:107) u ( t , · ) (cid:107) H ≤ (cid:18) t + 12 (cid:107) b (cid:107) C + (cid:107) b (cid:107) C (cid:19) e t (cid:107) b (cid:107) C (cid:107) u (cid:107) . Hence P ε maps bounded sets in L ( X ) onto bounded sets in H ( X ). As the embeddingof H ( X ) onto L ( X ) is compact by Rellich’s theorem, P ε is a compact operator.14 MATLAB code for the quadruple gyre example % computation of coherent sets for a quadruple gyre system % % Andreas Denner and Oliver Junge, TUM, 2015 M = 15; x = 2/M*(0:M-1)’; [X,Y] = meshgrid(x); % collocation points D = 1i*ifftshift(-(M-1)/2:(M-1)/2)’; % derivative in frequency space Dy = D*ones(1,M); Dx = -ones(M,1)*D’; ep = 0.02; L = epˆ2/2*(Dy.ˆ2+Dx.ˆ2); % Laplace operator %% vector field dl = 0.25; om = 2*pi; f = @(t,x) dl*sin(om*t).*x.ˆ2 + (1-2*dl*sin(om*t)).*x; df = @(t,x) 2*dl*sin(om*t).*x + 1-2*dl*sin(om*t); g = @(t,x,y) sin(pi*f(t,x)).*cos(pi*f(t,y)).*df(t,y); w = @(t,v) -1/2*(Dx.*fft2(-g(t,X,Y).*ifft2(v))+Dy.*fft2(g(t,Y,X).*ifft2(v)) + fft2(-g(t,X,Y).*ifft2(Dx.*v))+fft2(g(t,Y,X).*ifft2(Dy.*v))); v = @(t,x) reshape(w(t,reshape(x,M,M)),Mˆ2,1); %% intial values N = 5; % number of basis functions U0 = zeros(M,M,M,M); for k = 1:M, for l = 1:M, U0(k,l,k,l) = 1; end, end Y0 = U0(:,:,[1:(N-1)/2+1,(M-1)-((N-1)/2-2):M], ... [1:(N-1)/2+1,(M-1)-((N-1)/2-2):M]); %% time integration / construction of transfer operator t0 = 0; t1 = 10.25; m = 50; h = (t1-t0)/m; P = zeros(Mˆ2,Nˆ2); for l = 1:N for k = 1:N y0 = reshape(Y0(:,:,k,l),Mˆ2,1); P(:,N*(l-1)+k) = etdrk4(t0,m,h,L,v,y0); end end [UU,S,VV] = svd(P); diag(S); % compute singular values and vectors %% plot singular vectors figure(1); clf; p = 128; [Xp,Yp] = meshgrid(2*(0:p-1)/p); for j = 2 v = reshape(VV(:,j),N,N)’; vp = zeros(p); vp([1:(N+1)/2,p-(N-1)/2+1:p],[1:(N+1)/2,p-(N-1)/2+1:p]) = v*pˆ2/Nˆ2; vp = ifft2(vp,’symmetric’); %vp = vp/(sum(abs(vp(:,j)))/pˆ2*(2)ˆ2); vp=vp/max(max(abs(vp))); surf(Xp,Yp,vp), shading flat, view(90,90), caxis([-1 1]), axis equal, hold on, axis off,axis tight % contour3(Xp,Yp,vp,[-0.5 -0.5],’linecolor’,[0.95 0.95 0.95],’linewidth’,4); % unix([’convert -trim tmp.png ’ filename]); end eferences [1] R. Banisch and P. Koltai. Understanding the geometry of transport: diffu-sion maps for lagrangian trajectory data unravel coherent sets. arXiv preprintarXiv:1603.04709 , 2016.[2] J. P. Boyd. Chebyshev and Fourier spectral methods . Courier Dover Publications,2001.[3] S. Cox and P. Matthews. Exponential time differencing for stiff systems.
Journal ofComputational Physics , 176(2):430–455, 2002.[4] M. Dellnitz, G. Froyland, and O. Junge. The algorithms behind
GAIO – Set orientednumerical methods for dynamical systems. In B. Fiedler, editor,
Ergodic Theory,Analysis, and Efficient Simulation of Dynamical Systems , pages 145–174. Springer,2001.[5] M. Dellnitz and O. Junge. On the approximation of complicated dynamical behavior.
SIAM Journal on Numerical Analysis , 36(2):491–515, 1999.[6] L. Evans.
Partial Differential Equations . Graduate studies in mathematics. AmericanMathematical Society, 2010.[7] G. Froyland. An analytic framework for identifying finite-time coherent sets in time-dependent dynamical systems.
Physica D: Nonlinear Phenomena , 250:1–19, 2013.[8] G. Froyland and M. Dellnitz. Detecting and locating near-optimal almost-invariantsets and cycles.
SIAM Journal on Scientific Computing , 24(6):1839–1863, 2003.[9] G. Froyland, O. Junge, and P. Koltai. Estimating long term behavior of flows with-out trajectory integration: the infinitesimal generator approach.
SIAM Journal onNumerical Analysis , 51(1):223–247, 2013.[10] G. Froyland, S. Lloyd, and N. Santitissadeekorn. Coherent sets for nonautonomousdynamical systems.
Physica D , 239:1527–1541, 2010.[11] G. Froyland and K. Padberg. Almost-invariant sets and invariant manifolds – con-necting probabilistic and geometric descriptions of coherent structures in flows.
Phys-ica D , 238:1507–1523, 2009.[12] G. Froyland and K. Padberg-Gehle. Almost-invariant and finite-time coherent sets:directionality, duration, and diffusion. In
Ergodic Theory, Open Dynamics, and Co-herent Structures , pages 171–216. Springer, 2014.[13] G. Froyland, N. Santitissadeekorn, and A. Monahan. Transport in time-dependentdynamical systems: Finite-time coherent sets.
Chaos , 20:043116, 2010.[14] A. Hadjighasem, D. Karrasch, H. Teramoto, and G. Haller. A spectral clusteringapproach to lagrangian vortex detection. arXiv preprint arXiv:1506.02258 , 2015.[15] G. Haller. Distinguished material surfaces and coherent structures in three-dimensional fluid flows.
Physica D , 149:248–277, 2001.1616] G. Haller. A variational theory of hyperbolic Lagrangian coherent structures.
PhysicaD , 240:574–598, 2011.[17] G. Haller and G. Yuan. Lagrangian coherent structures and mixing in two-dimensional turbulence.
Physica D , 147:352–370, 2000.[18] W. Huisinga and B. Schmidt. Metastability and dominant eigenvalues of transfer op-erators. In
New Algorithms for Macromolecular Simulation , pages 167–182. Springer,2006.[19] O. Junge, J. E. Marsden, and I. Mezic. Uncertainty in the dynamics of conservativemaps. In
Decision and Control, 2004. CDC. 43rd IEEE Conference on , volume 2,pages 2225–2230. IEEE, 2004.[20] A.-K. Kassam and L. N. Trefethen. Fourth-order time-stepping for stiff pdes.
SIAMJournal on Scientific Computing , 26(4):1214–1233, 2005.[21] A. Lasota and M. Mackey.
Chaos, Fractals, and Noise: Stochastic Aspects of Dy-namics . Number Bd. 97 in Applied Mathematical Sciences. Springer, 1993.[22] T. Y. Li. Finite Approximation for the Frobenius-Perron Operator. A Solution toUlam’s Conjecture.
J. Approx. Theory , 17:177–186, 1976.[23] J.-C. Nave. Computational science and engineering. 2008, (Massachusetts Instituteof Technology: MIT OpenCouseWare), http://ocw.mit.edu (Accessed July 13,2015). License: Creative Commons BY-NC-SA.[24] B. Oksendal. Stochastic differential equations: an introduction with applications.
Springer, New York , 2003.[25] C. Sch¨utte, A. Fischer, W. Huisinga, and P. Deuflhard. A direct approach to confor-mational dynamics based on hybrid monte carlo.
Journal of Computational Physics ,151(1):146–168, 1999.[26] S. M. Ulam.
Problems in modern mathematics . Courier Dover Publications, 2004.[27] T. A. Zang. On the rotation and skew-symmetric forms for incompressible flowsimulations.