Estimating long term behavior of flows without trajectory integration: the infinitesimal generator approach
EEstimating long term behavior of flows without trajectoryintegration: the infinitesimal generator approach
Gary Froyland ∗ Oliver Junge † P´eter Koltai ‡ November 10, 2018
Abstract
The long-term distributions of trajectories of a flow are described by invariant densities,i.e. fixed points of an associated transfer operator. In addition, global slowly mixing struc-tures, such as almost-invariant sets, which partition phase space into regions that are almostdynamically disconnected, can also be identified by certain eigenfunctions of this operator.Indeed, these structures are often hard to obtain by brute-force trajectory-based analyses. Ina wide variety of applications, transfer operators have proven to be very efficient tools for ananalysis of the global behavior of a dynamical system.The computationally most expensive step in the construction of an approximate transferoperator is the numerical integration of many short term trajectories. In this paper, wepropose to directly work with the infinitesimal generator instead of the operator, completelyavoiding trajectory integration. We propose two different discretization schemes; a cell baseddiscretization and a spectral collocation approach. Convergence can be shown in certaincircumstances. We demonstrate numerically that our approach is much more efficient thanthe operator approach, sometimes by several orders of magnitude.
Analysis of the long-term behavior of flows can be broadly classified into geometric methodsand statistical methods. Geometrical methods include the determination of fixed points, periodicorbits, and invariant manifolds. Invariant manifolds of fixed points or periodic orbits act as barriersto transport as trajectories may not cross the manifolds transversally. Statistical methods includedetermining the distribution of points in very long trajectories of a very large set of initial points(i.e. a physical invariant measure [29], often possessing an invariant density) and the identificationof meta-stable or almost-invariant sets [8, 15, 14]. Almost-invariant sets partition the phase spaceinto almost dynamically disconnected regions and are important for revealing global dynamicalstructures that are often invisible to an analysis of trajectories. These metastable dynamics alsogo under the names of persistent patterns, or strange eigenmodes, both of which are realisable aseigenfunctions of a transfer operator. Frequently, the boundaries of maximally almost-invariantsets (those sets which are locally closest to invariant sets) coincide with certain invariant manifolds[16].Our focus in the present paper is on statistical methods, although we demonstrate via casestudies the relationships to geometrical methods. The most commonly used tool for statisticalmethods is the transfer operator (or Perron-Frobenius operator). Fixed points of the transferoperator correspond to invariant densities, while eigenfunctions corresponding to real positiveeigenvalues strictly less than one provide information on almost-invariant sets. In practice, onetypically constructs a finite-rank numerical approximation of a transfer operator and computes ∗ School of Mathematics and Statistics, University of New South Wales, Sydney, NSW 2052, Australia. † Faculty of Mathematics, Technische Universit¨at M¨unchen, 85748 Garching, Germany. ‡ Faculty of Mathematics, Technische Universit¨at M¨unchen, 85748 Garching, Germany. P´eter Koltai was partlysupported by the TopMath PhD program within the Elite Network of Bavaria and the TUM Graduate School. a r X i v : . [ m a t h . NA ] J a n arge spectral values and eigenfunctions for this finite-rank operator. The construction of thefinite-rank approximation requires the integration of many relatively short trajectories with initialpoints sampled over the domain of the flow. It is this use of short trajectories that gives thetransfer operator approach additional stability and accuracy when compared with computationsbased upon very long trajectories. Long trajectories continually accumulate small errors fromimperfect numerical integration and finite computer representation of numbers; these small errorsquickly grow in chaotic flows. While the transfer operator approach is very stable, it still requiresthe computation of many small trajectories which can be very time consuming in some systems.The approach we describe in the present work obviates the need for any trajectory integration atall and works directly with the vector field.Our approach exploits the fact that the evolution of probability densities u = u ( t ) can be de-scribed by generalized solutions of the abstract Cauchy problem ∂ t u = A u . The Perron–Frobeniusoperator is the evolution operator of this equation, and has the same eigenfunctions as the opera-tor A . The operator A is an unbounded hyperbolic (if the underlying dynamics is deterministic)or elliptic (if the deterministic dynamics is perturbed by white noise) partial differential operator.Standard techniques allow us to approximate the eigenmodes we are interested in: finite difference,finite volume, finite element and spectral methods yield such discretizations; see [25, 19, 21, 4, 5]and the references therein.An outline of the paper is as follows. In Section 2 we provide background on the infinitesi-mal operator arising from smooth vector fields and describe conditions under which the operatorgenerates a semigroup of transfer operators. In order to obtain formal results, we will requirethe addition of a small amount of diffusion to the deterministic flow. In the small diffusion set-ting we discuss existence of invariant densities and spectral results for the associated infinitesimaloperator. Section 3 describes a spectral Galerkin method for the approximation of the infinitesi-mal operator. We apply results from the numerical analysis of advection-diffusion PDEs to showthat our Galerkin method approximates the true eigenfunctions of the infinitesimal operator asour Galerkin basis becomes increasingly refined. We can also show that the convergence rate isspectral; that is, faster than any polynomial. Section 4 describes an Ulam-based Galerkin methodfor approximating the infinitesimal operator. This Ulam-based approach is new and shares somesimilarities with finite-difference schemes. While we cannot show convergence of this approxima-tion scheme, the numerical results obtained are extremely fast and accurate. Sections 5–8 detailthe practical application of our two approximation methods to flows in one-, two-, and three-dimensional domains. We demonstrate the spectral accuracy of our spectral Galerkin method andcompare with the accuracy of the Ulam-based Galerkin method. We also compare the accuracy vs.computational effort of our two new approaches with standard transfer operator approaches. Thenatural relationships between the outputs of our statistical methods and geometric objects suchas the vector field and invariant manifolds are also elucidated in each case study. We conclude inSection 9. Let the domain M ⊂ R d of our flow be a smooth compact manifold and m the (normalized)Lebesgue measure on M . Denote by F : M → R d the vector field generating the flow and byΦ t : M → M , t ∈ R , the flow, i.e. Φ t ( x ) represents the location of a trajectory beginning at x ∈ M after flowing for t ∈ R time units. One has that d Φ t ( x ) /dt = F ( x ). Note that – providedthat the components F i , i = 1 , . . . , d have continuous derivatives – the function Φ t : M → M is adiffeomorphism for every t ∈ R .Invariant sets are structures of dynamical interest. A set A ⊂ M is called invariant , iffΦ − t A = A for all t . Also, one asks how the flow changes probability measures. Sample x accordingto a probability measure µ ; the distribution of Φ t x is then given by µ ◦ Φ − t . Special attentionis to be drawn to invariant measures , which do not change under the dynamics ( µ = µ ◦ Φ − t ).Invariant measures µ are called ergodic if invariant sets have either zero or full measure, i.e. if A ⊂ M satisfies Φ − t A = A then µ ( A ) ∈ { , } . An even more restricted class of ergodic measures2re the physically relevant (or natural ) ones, satisfying (cid:90) ψ d µ = lim T →∞ T (cid:90) T ψ (cid:0) Φ t x (cid:1) d t for all ψ : M → R continuous and x ∈ U ⊂ M with m ( U ) >
0. One has that absolutelycontinuous ergodic measures are natural invariant measures. The density of an invariant measureis called the invariant density . When looking for invariant densities, one can rephrase the action of the flow on measures as anaction on densities. If f denotes the density of µ , then f is evolved by the flow as P t f ( x ) = f (cid:0) Φ − t x (cid:1) | D Φ − t x | , where | B | denotes | det B | for a matrix B . The linear operator P t : L ( M ) (cid:9) is known as a transferoperator or the Perron-Frobenius operator associated with the flow Φ. Note that invariant densitiesare fixed points of P t . For each t ∈ R we have(i) P t is linear,(ii) P t f ≥ f ≥ (cid:107)P t f (cid:107) = (cid:107) f (cid:107) for all f ≥
0, where (cid:107) · (cid:107) is the L norm.i.e. P t is a Markov operator . Moreover, these properties also imply (cid:107)P t (cid:107) ≤
1, hence P t is acontraction and its eigenvalues lie inside the unit disc. In fact, since the flow Φ t is a bijection, wehave (cid:107)P t f (cid:107) = (cid:107) f (cid:107) for all f ∈ L . To see this, note that (cid:90) A P t f dm = (cid:90) Φ − t A f dm (1)for t > A ⊂ M measurable. Now write f = f + − f − , with f ± ( x ) = max {± f ( x ) , } , then (cid:107) f (cid:107) = (cid:82) ( f + − f − ) = (cid:82) S + f + + (cid:82) S − f − , where S ± is the support of f ± . Now use properties (ii),(iii) and (1) with A = S ± . A consequence of this is that the eigenvalues of P t lie on the unitcircle. The transfer operator also inherits some (semi)group properties of the flow Φ t . Definition 2.1.
Let ( X, (cid:107) · (cid:107) ) be a Banach space. A one parameter family {T t } t ≥ of boundedlinear operators T t : X → X is called a semigroup on X , if(a) T = I ( I denoting the identity on X ),(b) T t + s = T t T s for all t, s ≥ (cid:107)T t (cid:107) ≤
1, the family is called a semigroup of contractions.If lim t → (cid:107)T t f − f (cid:107) = 0 for every f ∈ X , then T t is a continuous semigroup ( C semigroup ).The transfer operator P t is a C semigroup of contractions on L , see [20] for a proof (inparticular Remark 7.6.2 for the continuity). µ is absolutely continuous (with respect to m , if not specified further), if there is an 0 ≤ f ∈ L ( m ) such thatfor all m -measurable B ⊂ M we have µ ( B ) = (cid:82) B f d m . The function f is called the density of µ . efinition 2.2. For a semigroup T t we define the operator A : D ( A ) → X by A f = lim t → T t f − ft , f ∈ D ( A ) , with D ( A ) ⊂ X being the linear subspace of X where the above limit exists. The operator A iscalled the infinitesimal generator of the semigroup.For P t , the infinitesimal generator turns out to be (provided the F i are continuously differen-tiable) A P F f = − div( f F ) , (2)see [20]. The following result (see eg. Theorem 2.2.4 [23]) shows the connection between theeigenvalues of the semigroup operators and their infinitesimal generator: Theorem 2.3 (Spectral mapping theorem) . Let T t be a C semigroup and let A be its infinitesimalgenerator. Then e tσ ( A ) ⊂ σ (cid:0) T t (cid:1) ⊂ e tσ ( A ) ∪ { } , where σ ( · ) denotes the point spectrum of the operator. The corresponding eigenvectors are identical. This has important consequences for invariant densities:
Corollary 2.4.
The function f is a invariant density of P t for all t ≥ if and only if A P F f = 0 . According to the discussion at the end of Section 2.1, we have
Corollary 2.5.
The eigenvalues of A P F lie on the imaginary axis.
In this section we discuss almost-invariant sets and escape rates via a spectral analysis of theinfinitesimal generator. As remarked near the end of Section 2.1, the L spectrum of P t lies onthe unit circle, and lacks a spectral gap. Applying the spectral mapping theorem, we see that thespectrum of A must be pure imaginary. In the following discussion, to prove formal results, wewill add a small random perturbation to P t . Later, in the numerics section, we will see that ournumerical methods introduce a numerical diffusion that plays the role of creating a spectral gap. In many real world situations, a deterministic model of some physical system is not appropriate.Rather, one should account for the fact that certain external perturbations are present which mightbe unknown or for which a detailed model would be overly complicated. Often, it is appropriateto account for these influences by incorporating a small random perturbation into the description.From a theoretical point of view, this even facilitates the analysis of the system: Under certainassumptions, the transfer operator becomes compact.We therefore now leave the deterministic setting and assume that our dynamical system de-scribed by the ordinary differential equation dx/dt = F ( x ) is slightly stochastically perturbed by white noise . An exact mathematical treatment of this topic would require tools which are beyondthe scope of this work (for an exact derivation see [20], chapter 11). We are primarily interested inthe time evolution of probability densities. The following material highlights the relevant formalstatements and attempts to point out the intuition behind them. Instead of an ordinary differentialequation, we now deal with a stochastic differential equation (SDE) dXdt = F ( X ) + ε dWdt , (3)4ith ε > W being a d -dimensional Wiener process. The solutions are time-dependentrandom variables X ( t ) with values in R d . Just as in the deterministic case, we look at theevolution of density functions; now the density functions represent the distribution of randomoutcomes of the variables X ( t ). The density functions f satisfyProb( X ( t ) ∈ A ) = (cid:90) A f ( x, t ) d x. Assuming the existence of such a density f ( x, t ), where f ( · ,
0) = f is given, we may again definethe transfer operator as P tε f := f ( t, · ). If the vector field F is smooth enough we have followingcharacterization of the density function, the Fokker-Planck or Kolmogorov forward equation , cf.[20], 11.6: ∂f∂t = ε f − div( f F ) =: A ε f. (4) Proposition 3.1 ([1, 22]) . The operator A ε (with Neumann boundary conditions ) is the in-finitesimal generator of C semigroups P tε, on L , P tε, on L , and P tε, on C . We note that the operator P tε, is identical to the transfer operator P tε . Invariant densities
Again, we are particularly interested in distributions (densities) which donot change under the evolution. Those are again given by fixed points of P tε, or alternatively, byfunctions in the null space of A ε . One can show [30] that P tε, is compact and thus the null spaceof A ε is finite dimensional. Furthermore, due to the white noise, the support of the stochastictransition function associated to P tε, is unbounded, and the null space of A ε is one-dimensional,i.e. there is a unique invariant density (see again [30], Theorem 1), characterizing the long termdynamical behavior of (3). We call a set A ⊂ M almost-invariant with respect to a (not necessarily invariant) probabilitymeasure ν (cf. [15, 14]), if ρ tν ( A ) := ν (Φ − t ( A ) ∩ A ) ν ( A ) ≈ t . The analogous expression for Φ perturbed by a small random perturbation is ρ tν,ε ( A ) := Prob ν ( X (0) ∈ A, X ( t ) ∈ A )Prob ν ( X (0) ∈ A ) ≈ t .We can alternatively characterize this property of a set using an infinitesimal representationof P tε . To this end let f ∈ L and consider a measurable set A ⊂ M . We define the functional A ε,A : L ⊃ D A ( A ε ) → R by A ε,A f := lim t → (cid:90) A P tε f − ft dm, f ∈ D A ( A ε ) , (7)where D A ( A ε ) is the linear subspace of L where the above limit exists. Let f A := f χ A / ( (cid:82) f χ A dm ). Proposition 3.2.
Let ≤ f ∈ D A ( A ε ) be the density of the probability measure ν . Then ρ tν,ε ( A ) = 1 + ( A ε,A f A ) t + o ( t ) (8) for measurable A as t → . The dynamical reason why Neumann boundary conditions are chosen here is discussed later. For two functions f, g : R → R we say “ f ( t ) = o ( g ( t )) as t → t → f ( t ) g ( t ) = 0. roof. We have ρ tν,ε ( A ) − t = Prob ν ( X (0) ∈ A, X ( t ) ∈ A ) − ν ( A ) t Prob ν ( X (0) ∈ A )= (cid:82) A P tε ( f χ A ) dm − (cid:82) A f dmt ν ( A ) = (cid:90) A P tε f A − f A t dm, where the second equation follows from the fact that X (0) is distributed according to ν (withdensity f ), and X (0) ∈ A ; which gives the density f χ A for X (0), see Section 3.1. Thus for t → A will be almost invariant with respect to the measure with density f ∈ L , if A ε,A f A ≈ Theorem 3.3.
Suppose that the generator A ε possesses a real eigenvalue λ < with correspondingeigenfunction f . Then (cid:82) f dm = 0 . Define A + = { f ≥ } , A − = { f < } . Let f A + , f A − ∈ D A ( A ) .Then A ε,A + | f A + | + A ε,A − | f A − | = λ. Proof.
By Proposition 5.7 in [8] and the spectral mapping theorem we have that ρ tν,ε ( A + ) + ρ tν,ε ( A − ) = exp( tλ ) + 1 , where ν is the probability measure with density | f | . Using Proposition 3.2 we obtain 1 + A ε,A + | f A + | t + A ε,A − | f A − | t = exp( tλ ) + o ( t ), i.e. A ε,A + | f A + | + A ε,A − | f A − | = exp( tλ ) − t + o (1)and for t → λ ≈ A ε,A ± | f A ± | ≈
0, which means that the sets A + and A − will bealmost-invariant with respect to the probability measure with density | f | . Other techniques forextracting sets A + , A − from the eigenfunction f may be found in [15, 14]. The papers [15, 14]also discuss almost-invariance with respect to physical invariant probability measures, a propertythat is particularly meaningful when studying typical dynamical behavior. In all cases, the basisfor these methods are eigenfunctions of A ε corresponding to (real) eigenvalues close to 0.If A ε has a complex eigenvalue with real part close to zero, then the corresponding complexeigenfunction may also be used to construct almost-invariant sets. Lemma 3.4.
Let A ε f = λ A f with λ A ∈ C and let f re and f im denote the real and imaginarypart of f , respectively. Let t > be such that e tλ A = λ P F ∈ R . Then P tε f re = λ P F f re and P tε f im = λ P F f im .Proof. From the proof of Theorem 2.2.4 [23] we have P tε f = λ P F f . Note, that P tε : L ( M, R ) (cid:9) .By linearity we have P tε f = P tε f re + i P tε f im . Thus, λ P F f re (cid:124) (cid:123)(cid:122) (cid:125) ∈ R + i λ P F f im (cid:124) (cid:123)(cid:122) (cid:125) ∈ R = P tε f re (cid:124) (cid:123)(cid:122) (cid:125) ∈ R + i P tε f im (cid:124) (cid:123)(cid:122) (cid:125) ∈ R . The claim follows immediately.Hence, if for a t > ≈ e tλ A ∈ R , then the real (and imaginary) part of thecorresponding eigenfunction yields a decomposition of the phase space into almost invariant setsin the sense of Theorem 3.3. 6 .3 Escape rates We have seen in the previous section that there is a connection between eigenvalues of A ε close tozero and almost-invariant sets. There is also a strong connection between the escape rate of thesets A + , A − and the corresponding eigenvalue of A ε .In the deterministic setting ( ε = 0) the upper (resp. lower) escape rate of a set A under thetime- t map of the flow Φ t is defined as:¯ E ( A ) = − lim inf k →∞ k log m ( A ∩ Φ − t A ∩ · · · ∩ Φ − kt A ) (9)E( A ) = − lim sup k →∞ k log m ( A ∩ Φ − t A ∩ · · · ∩ Φ − kt A ) (10)If both limits are equal, we call E ( A ) = lim k →∞ k log m ( A ∩ Φ − t A ∩ · · · ∩ Φ − kt A ) the escaperate from A . The escape rate is the asymptotic exponential rate of loss of m -mass from the set A . One might expect that almost-invariant sets have low escape rate and vice-versa, but simplecounterexamples in [17] show that this is not the case. This is because the notion of almost-invariance as defined above is a finite-time property, while escape rate is an asymptotic quantity.The ε > E ε ( A ) := − lim inf k →∞ k log Prob m ( X (0) ∈ A, X ( t ) ∈ A, . . . , X ( kt ) ∈ A ) (11)= − lim inf k →∞ k log (cid:90) A ( P tε,A ) k (1) dm. (12)where P tε,A ( f ) := P tε ( f χ A ) is the standard restriction of the operator P tε to A [24]. The followingtheorem provides a recipe for constructing sets with low escape rate from eigenfunctions of A ε with real eigenvalues. Theorem 3.5.
Suppose that A ε f = λf for some λ < and f ∈ L ∞ ( m ) . Then ¯ E ε ( A + ) ≤ − tλ and ¯ E ε ( A − ) ≤ − tλ, where A + = { f ≥ } and A − = { f < } .Proof. Let f be scaled such that (cid:82) | f | dm = 2, and define the (signed) measure ν as ν ( A ) = (cid:82) A f dm . Then ν + = ν (cid:12)(cid:12) A + and ν − = − ν (cid:12)(cid:12) A − are both probability measures. For an event E defineProb ν ( E ) := Prob ν + ( E ) − Prob ν − ( E ). Since f is an eigenfunction of P tε with eigenvalue e λt , wehave e λkt = e λkt (cid:90) A + f dm = (cid:90) A + P ktε f dm = Prob ν (cid:0) X ( kt ) ∈ A + (cid:1) = Prob ν (cid:0) X ( (cid:96)t ) ∈ A + , (cid:96) = 0 , . . . , k (cid:1) + k − (cid:88) n =0 Prob ν (cid:0) X ( nt ) ∈ A − , X ( (cid:96)t ) ∈ A + , (cid:96) = n + 1 , . . . , k (cid:1)(cid:124) (cid:123)(cid:122) (cid:125) =: p n . It follows for the summands p n in the above sum: p n = Prob P ntε ν (cid:0) X (0) ∈ A − , X ( (cid:96)t ) ∈ A + , (cid:96) = 1 , . . . , k − n (cid:1) = e λnt Prob ν (cid:0) X (0) ∈ A − , X ( (cid:96)t ) ∈ A + , (cid:96) = 1 , . . . , k − n (cid:1) = − e λnt Prob ν − (cid:0) X ( (cid:96)t ) ∈ A + , (cid:96) = 1 , . . . , k − n (cid:1) , where in the first equation the action of the transfer operator P tε on the measure ν is definedthrough its action on the density f of ν . The second equation follows from ν being an eigenmeasure.Clearly, p n is non-positive, and thusProb ν (cid:0) X ( (cid:96)t ) ∈ A + , (cid:96) = 0 , . . . , k (cid:1) ≥ e λkt holds. The rest of the proof follows the lines of the one of Theorem 2.4 in [17].7he recipe of Theorem 3.5 is in fact the same as in Theorem 3.3. The difference is the measureused: in Theorem 3.3, almost-invariance is computed with respect to the particular measure ν withdensity | f | where f is the eigenfunction in question, whereas in Theorem 3.5 escape is computedwith respect to Lebesgue measure. Remark 3.6.
A more natural notion of escape rates in the time-continuous case would be thefollowing: ¯ E ε ( A ) := − lim inf t →∞ t Prob m (cid:0) X ( s ) ∈ A, s ∈ [0 , t ) (cid:1) . However, for this definition it is more complicated to obtain similar results to the one in Theo-rem 3.5, and a fuller discussion will appear elsewhere. One can view (11) with t = 1 as the “timesampled version” of this continuous-time definition. Having seen that certain eigenpairs of the transfer operator (resp. infinitesimal generator) carrythe information we are seeking, we describe in the sequel our proposed approximation of theseeigenpairs. To this end, we define finite dimensional approximation spaces and consider the eigen-value problem projected onto these spaces. Throughout this section we assume that the underlyingdeterministic vector field F is smooth. We describe here the “standard” Ulam approach; see the surveys [13, 6] for more details. Wepartition M into d -dimensional connected, positive volume subsets { B , . . . , B n } . Typically, each B i will be a hyperrectangle or simplex to simplify computations. As an approximation space weconsider the space ∆ n = sp { χ B , . . . , χ B n } of functions which are piecewise constant on the cellsof the partition. Let π n : L → ∆ n , π n f = (cid:80) ni =1 1 m ( B i ) (cid:82) B i f dm χ B i , be the L -orthogonalprojection onto ∆ n . We let P tn : ∆ n → ∆ n , P tn := π n P t , be the approximate Frobenius-Perronoperator. Note that P tn χ B i = π n P t χ B i = (cid:80) nj =1 1 m ( B j ) (cid:82) B j P t χ B i dm χ B j , i.e. the matrix repre-sentation P tn ∈ R n × n of P tn with respect to the basis χ B , . . . , χ B n and multiplication on the leftis ( P tn ) ij = 1 m ( B j ) (cid:90) B j P t χ B i dm = m ( B i ∩ Φ − t B j ) m ( B j ) . This matrix is easily constructed numerically using eg. GAIO [6].In the stochastic setting, letting P tε,n : ∆ n → ∆ n , P tε,n := π n P tε , one obtains( P tε,n ) ij = 1 m ( B j ) (cid:90) B j P tε χ B i dm as above. In principle one can use Monte-Carlo integration to compute these entries, however,this is not particularly efficient. We partition M as in the standard Ulam’s method. We will see that the numerical scheme itselfintroduces some diffusion and we therefore consider the deterministic generator A (i.e. ε = 0).We wish to construct an operator A n : ∆ n → ∆ n that is close in some sense to the operator A . Motivated by Ulam’s method, one would like to form π n A , which unfortunately does not exist,because ∆ n (cid:42) D ( A ), cf. (2). Instead of differentiating w.r.t. time and then doing the projection,we swap the order of these operations. Let us build the Ulam approximation P tn first, which will8 ot be a semigroup any more, but for fixed t it approximates P t . Taking the time derivative, ourcandidate approximate operator is A n f := lim t → (cid:18) π n P t π n f − π n ft (cid:19) . (13)The following lemma emphasizes the intuition behind this definition: if P tn is a sufficiently goodapproximation (for small t ) of the Markov jump process generated by the dynamics on the sets B i , then A n will be the generator of this process. To be exact, the following is the case: Lemma 4.1.
The matrix representation of A n : ∆ n (cid:9) with respect to the basis χ , . . . , χ n undermultiplication on the left is ( A n ) ij = lim t → m ( B i ∩ Φ − t B j ) t · m ( B j ) , i (cid:54) = j ; lim t → m ( B i ∩ Φ − t B i ) − m ( B i ) t · m ( B i ) , otherwise. (14) Proof.
We consider the action of P t on χ B i .lim t → π n P t χ B i − χ B i t = lim t → n (cid:88) j =1 m ( B j ) (cid:32)(cid:90) B j P t χ B i − χ B i t dm (cid:33) χ B j = lim t → (cid:88) j (cid:54) = i m ( B j ) (cid:32)(cid:90) B j P t χ B i t dm (cid:33) χ B j + lim t → m ( B i ) (cid:18)(cid:90) B i P t χ B i − χ B i t dm (cid:19) χ B i = lim t → (cid:88) j (cid:54) = i m ( B j ) (cid:32)(cid:90) Φ − t B j χ B i t dm (cid:33) χ B j + lim t → m ( B i ) (cid:18)(cid:90) Φ − t B i χ B i t dm − (cid:90) B i χ B i t dm (cid:19) χ B i = (cid:88) j (cid:54) = i lim t → m ( B i ∩ Φ − t ( B j )) t · m ( B j ) χ B j + lim t → m ( B i ∩ Φ − t B i ) − m ( B i ) t · m ( B i ) χ B i Thus under left multiplication we obtain (14).
Remark 4.2.
Lemma 4.1 states that ( A n ) ij is the outflow rate of mass from B i into B j .Let R tn := exp ( t A n ) = I − π n + exp( t A n | V n ) π n denote the semigroup generated by A n . Then R tn and P t are near to each other in the following sense, cf. [18]: Proposition 4.3. As t → R tn f − π n P t f = O ( t ) (15) for all f ∈ ∆ n . The following lemma allows us to construct A n without the computation of the flow Φ t . Lemma 4.4.
For i (cid:54) = j , define n ij to be the the unit normal vector pointing out of B i into B j if B i ∩ B j is a d − -dimensional face, and the zero vector otherwise. The matrix representation of A n : ∆ n (cid:9) with respect to the basis χ , . . . , χ n under multiplication on the left is ( A n ) ij = m ( B j ) (cid:90) B i ∩ B j max { F ( x ) · n ij , } dm d − ( x ) , i (cid:54) = j ; − (cid:88) j (cid:54) = i m ( B j ) m ( B i ) ( A n ) ij , otherwise. (16) For two functions f, g : R → R we say “ f ( t ) = O ( g ( t )) as t → t → | f ( x ) || g ( x ) | < ∞ . roof. From (14) we have for i (cid:54) = j that A n,ij = lim t → m ( B i ∩ Φ − t B j ) t · m ( B j ) . Denoting M ij ( t ) = m ( B i ∩ Φ − t B j ) we have that A n,ij = M (cid:48) ij (0) /m ( B j ) where the prime denotes differentiation with respectto t . The quantity M (cid:48) ij (0) is simply the rate of flux out of B i through the face B i ∩ B j into B j and so M (cid:48) ij (0) = (cid:82) B i ∩ B j max { F ( x ) · n ij , } dm d − ( x ).For the diagonal elements A n we have A n,ii = lim t → m ( B i ∩ Φ − t B i ) − m ( B i ) t · m ( B i ) . Note that m ( B i ) − m ( B i ∩ Φ − t B i ) = m ( B i \ Φ − t B i ). Clearly B i \ Φ − t B i = B i ∩ (cid:83) j (cid:54) = i Φ − t B j = (cid:83) j (cid:54) = i B i ∩ Φ − t B j modulo sets of Lebesgue measure zero. Thus, m ( B i ) − m ( B i ∩ Φ − t B i ) = (cid:80) j (cid:54) = i m ( B i ∩ Φ − t B j ).Now, by (14), A n,ii = − lim t → (cid:80) j (cid:54) = i m ( B i ∩ Φ − t B j ) m ( B i ) = − (cid:80) j (cid:54) = i m ( B j ) m ( B i ) A n,ij . In one dimension, (16) has a particularly simple form.
Corollary 4.5.
Let M = T . Assume (without loss ) that F ( x ) ≥ . Let x < x < · · ·
0. The discretized FPO P tn can be related to a non-deterministic dynamical system, which, after mapping the initial point, adds some uncertainty toproduce a uniform distribution of the image point in the box where it landed; see [12]. Thus, thisuncertainty resulting from the numerical discretization, equivalent to the numerical diffusion inthe upwind scheme, can be viewed as the reason for robust behavior — stability.Finally, we show that our constructions (16) and (17) always provide a solution to the system A n f = 0 for some f ∈ ∆ n . If F (cid:3) F (cid:2)
0, we have one or more stable fixed points, and every trajectory converges to one of them.Hence, there is no interesting statistical behavior to analyze. emma 4.7. There exists a nonnegative, nonzero f ∈ ∆ n so that A n f = 0 .Proof. Let D n,ij = m ( B i ) δ ij and note that Q n := D − n A n D n satisfies( Q n ) ij = (cid:40) m ( B j ) m ( B i ) A ij , i (cid:54) = j ; − (cid:80) j (cid:54) = i m ( B j ) m ( B i ) ( A n ) ij , otherwise. (18)Note that all row sums of Q n equal zero. Let c = (cid:80) j (cid:54) = i Q n,ij . The matrix ˆ Q n := Q n + cI is nonnegative with all row sums equal to c . By the Perron-Frobenius theorem [3], the largesteigenvalue of ˆ Q n is c (of multiplicity possibly greater than 1) and one of the corresponding lefteigenvectors pair u n is nonnegative. Clearly u n is a left eigenvector of Q n corresponding to theeigenvalue 0 and thus D n u n is a nonnegative left eigenvector corresponding to 0 for A n . The eigenfunctions of A ε are smooth. This motivates the use of smooth approximation functions,e.g. polynomials. We here outline the general principles of spectral collocation methods; for amore thorough presentation we refer to [5], [4] or [27].Choose a family of approximation spaces { V n } n ∈ N , such that V n ⊂ C ∞ ( M ) for all n . Dependingon the type of the phase space, we use two different approximation spaces. We introduce them inone dimension; the multidimensional ones can then be constructed by tensor products. In bothcases, the approximation space comes with an associated set of collocation nodes: • Periodic domain/uniform grid.
We have M = T and restrict ourselves to odd values of n .Then the basis we choose for V n is (cid:8) e ikx (cid:9) − n/ − ≤ k ≤ n/ . The associated collocation nodes are the uniform grid { , /n, . . . , ( n − /n } . • Standard domain/Chebyshev grid.
Here, M = [ − , V n is spanned by themonomials of order 0 to n . We use Chebyshev polynomials as basis functions: { cos ( k arccos ( x )) } ≤ k ≤ n , together with the Chebyshev grid {− cos(2 πj/n ) , j = 0 , . . . , n } , as collocation nodes.Let f ∈ V n and I n : C ∞ → V n be the interpolation operator for the given collocation nodes. Wedefine the approximate generator by A ε,n f := I n A ε f. For both cases we have following:
Theorem 4.8 (Spectral accuracy, [5]) . For f ∈ C ∞ ( M ) let f n be the best approximation of f in V n w.r.t. the supremum norm (cid:107) · (cid:107) ∞ . Then for each k ∈ N there is a c k > such that (cid:107) f − f n (cid:107) ∞ ≤ c k n − k for all n ∈ N . (19)Convergence then follows from standard results on the analysis of Galerkin methods for ellipticdifferential operators and spectral approximation (cf. also [18]): Theorem 4.9.
Let M be a compact tensor product domain with smooth boundary. Let F be asmooth vector field on M and A ε , A ε,n defined as above, with Neumann (resp. periodic) boundaryconditions. Then the eigenvalues and eigenfunctions of A ε are approximated by the correspondingeigenvalues and eigenfunctions of A ε,n with spectral accuracy. If ˆ Q n is aperiodic (there exists k such that ˆ Q kn >
0) then the eigenvalue c is simple and the correspondingeigenvector is positive (see [3]). hy Neumann boundary conditions? The objects that we are approximating are densities,so there should be no “loss of mass” as time goes on. In a closed physical system, the flow normalto the boundary is zero. Hence, there is no loss of mass through advection. The physical meaningof diffusion is a flow of mass in the direction opposite to the gradient and proportional to itsmagnitude. Thus, no loss of mass via diffusion translates into Neumann boundary conditions forthe densities: their normal gradients have to be zero at the boundary. Equivalently, one may write(4) as a continuity equation, ∂ t f = div( I ) , with I := ε ∇ f − f F being the probability flow or probability current. The condition I | ∂M = 0leads precisely to Neumann boundary conditions. Algorithm 1 (Ulam’s method) .
1. Partition M into connected sets { B , . . . , B n } of positive volume. Typically, each B i will bea hyperrectangle or simplex.2. Choose t > P tn or P tε,n as described in Section 4.1.3. Estimates of invariant densities for Φ t are given by fixed points of P tn (or P tε,n ): Let vP tn = v .Then f := (cid:80) ni =1 v i χ B i satisfies P tn f = f .4. Similarly, eigenvectors of P tn corresponding to real eigenvalues λ ≈ P tn typically is a sparse matrix and the number of nonzero entries per column willbe determined by Lipschitz constants of Φ t . Algorithm 2 (Ulam’s method for the generator) .
1. Partition M into connected sets { B , . . . , B n } of positive volume. Typically each B i will bea hyperrectangle or simplex.2. Compute ( A n ) ij = m ( B j ) (cid:90) B i ∩ B j max { F ( x ) · n ij , } dm d − ( x ) , i (cid:54) = j , − (cid:88) j (cid:54) = i m ( B j ) m ( B i ) ( A n ) ij , otherwise,where one uses standard quadrature rules to estimate the integral.3. Estimates of invariant densities lie in the null space of A n . Let vA n = 0 (the existence ofsuch a v is guaranteed by Lemma 4.7), then f := (cid:80) ni =1 v i χ B i satisfies A n f = 0.4. Similarly, eigenvectors of A n corresponding to large real eigenvalues λ < A n is a sparse matrix since A n,ij = 0 if B i and B j do notshare a boundary. Algorithm 3 (Spectral collocation for the generator) .
12. Set up a grid x , . . . , x n , which is the Chebyshev grid x j = − cos (2 πj/n ), j = 0 , . . . , n , if M = [ − , x j = j/n , if M = T . (If the state space M is an affinetransformation of the above ones, the grid is transformed analogously.) For multidimensionaltensor product spaces the grid is the tensor product of the one dimensional grids.2. Let { (cid:96) , . . . , (cid:96) n } denote the Lagrange basis on the grid, i.e. all (cid:96) i are polynomials of maximaldegree n with (cid:96) i ( x j ) = δ ij . Denote the interpolation on the grid by I n , and define thediscretization matrix obtained by collocation as( A ε,n ) ij = ( I n A ε (cid:96) j )( x i ) . Note, that the computation of these matrix entries is very simple. If M = T , we may switchbetween the evaluation space and the frequency space simply by the fast Fourier transfor-mation. Multiplication by the vector field F is pointwise multiplication in the evaluationspace, computing the derivative is a diagonal scaling in the frequency space. This can besimply extended to M = [0 ,
1] (cf. [27] Chapter 8), as well as for a multidimensional phasespace.3. Compute the left eigenvector v of A ε,n at the eigenvalue of smallest magnitude. Then, (cid:80) ni =0 v i (cid:96) i approximates the invariant density.4. Similarly, eigenvectors of A ε,n corresponding to large real eigenvalues λ < Once the matrix approximation of the operator or the generator has been computed, we then haveto solve the corresponding eigenvalue problem. In Ulam’s method one seeks dominant eigenvalues,and standard Arnoldi type methods easily provide the interesting part of the spectrum (we usethe eigs function in Matlab).For the approximate generator, we look for eigenvalues with small modulus. Here, the Arnoldiiteration in eigs uses inverse iteration, i.e. repeatedly solves linear systems of the type A n u = b .Since the matrix A n is singular, this is not a well posed problem and we might expect some diffi-culties with these iterative methods. A possibility for overcoming them is to solve the eigenvalueproblem for ¯ A n := A n − σI n × n , with a suitable shift σ . This merely introduces a shift of thespectrum and the eigenfunctions stay the same. For n big enough , the spectrum of ¯ A n is goingto be well separated from zero. The size n of the linear system may be large, so a direct solutionvia LU factorization may not be feasible. However, in Ulam’s method for the generator, A n is asparse matrix, hence iterative methods (e.g. GMRES) can be used.For spectral collocation, in general, the matrix A n will not be sparse. Nevertheless, in thecase where the domain is periodic, the matrix-vector product A n x can be computed fast using theFFT. For d dimensional systems in general, (16) shows us that setting up the approximate generatorrequires the computation of ( d − d -dimensional integrals. Of course, the gain is biggest in one dimension (zero dimensionalintegrals are one-point-evaluations). Equally importantly, (16) does not require integration of thevector field.For Ulam’s method and Ulam’s method for the generator, the matrix entries should be com-puted to an accuracy O ( n − /d ); otherwise we cannot expect the approximate eigenfunctions to If A n is obtained from the Ulam-based approach, the spectrum of A n , which is the generator of a Markov jumpprocess, is in the left complex half plane. This does not hold in general for the spectral collocation approach, herewe rely on the fast convergence of eigenvalues. O ( n − /d ) for n partitionelements in a d dimensional space). Assuming the quadrature rule used to compute the entries toalso suffer from the curse of dimension, we expect to need O ( n ) and O ( n ( d − /d ) function evalu-ations for each matrix entry in the Ulam-matrix, and the matrix arising from Ulam’s method forthe generator, respectively. This explains the third column of Table 1 below. A collection of examples is presented in the following. The examples range over three phasespace dimensions, different boundary types, and chaotic and regular dynamics. We find thatthe generator approaches are typically superior to the standard Ulam approach, producing moreaccurate results with less computation time. If the vector field is infinitely differentiable, thespectral collocation approach can be extremely accurate. In our final example, we find that forflows on complicated attractors, the standard Ulam method can take advantage of the attractingdynamics to produce better estimates of invariant densities than the generator approaches. In allcases, the operator based methods are far superior to direct trajectory integration for estimatingthe invariant density. Other structures such as almost-invariant sets can only be identified usingthe operator approaches.
We start with a one dimensional example, a flow on the unit circle. The vector field is given by F ( x ) = sin(4 πx ) + 1 . ,x ∈ T = [0 ,
1] with periodic boundary conditions, and we wish to compute the invariant densityof the system. Recall that an invariant density f ∈ L ( T ) needs to satisfy A f = 0, where A f ( x ) = − d( f F )( x ) / d x . The unique solution to this equation is f ∗ ( x ) = C/F ( x ), C being anormalizing constant (i.e. such that (cid:107) f (cid:107) = 1). The honours thesis [26] numerically investigatedthe estimation of the invariant density for this flow, and other two-dimensional flows, by findingapproximate eigensolutions of the infinitesimal generator using finite difference and finite elementmethods. Here we use the three methods discussed in the previous section to approximate f ∗ andcompare their efficiency relative to one another and a histogram of a long simulation:1. the classical method of Ulam for the Frobenius-Perron operator,2. Ulam’s method for the generator and3. spectral collocation for the generator.As we have a periodic domain and the vector field is infinitely smooth, we expect spectral colloca-tion to perform very well. Figure 1 shows the true invariant density (dashed line), together withits approximations by the four methods. In this 1D case, the three methods can each be realized in a few lines of Matlab. For illustrationpurposes, we include the code here.
Ulam’s method
F = @(x) sin(4*pi*x)+1.1; % definition of the vector fieldn = 32; x = 1/(2*n^2):1/(n^2):1-1/(2*n^2); % nodes for spatial integration[t,y] = ode23(@(t,x) F(x),[0 2/n],x’); % time integration of the nodesI = ceil(max(n*mod(y(end,:),1),1)); % cell indices of image pointsJ = reshape(ones(n,1)*(1:n),1,n*n); % construction of transition matrixP = sparse(I,J,1/n,n,n); v,d] = eig(full(P)); % spectrum of transition matrixf_n = @(x) abs(v(ceil(max(x*n,1)),1)); % plot approximate inv. densityL1 = quadl(f_n,0,1); fplot(@(x)f_n(x)/L1,[0,1]) The error in Ulam’s method decreases like O ( n − ) for smooth invariant densities [10]. Thus,we need to compute the transition rates between the intervals to an accuracy of O ( n − ) (sinceotherwise we cannot expect the approximate density to have a smaller error). To this end, we usea uniform grid of n sample points in each interval. This leads to O ( n ) evaluations of the vectorfield. For the numbers in Figure 1 we only counted each point once, i.e. we neglected the fact thatfor the time integration we have to perform several time steps per point. den s i t y den s i t y den s i t y den s i t y Figure 1: True invariant density (dashed line), approximation by histogramming a long simulation(top left, 1024 iterates), Ulam’s method (top right, 32 sample points per partition element, 32partition elements), Ulam’s method for the generator (bottom left, 32 partition elements) andspectral collocation (bottom right, 32 collocation points).
Ulam’s method for the generator
F = @(x) sin(4*pi*x)+1.1; % definition of the vector fieldn = 32; Fx = F(0:1/n:1-1/n); % evaluation on the boundary nodesA = n*sparse(1:n,1:n,-Fx) + sparse(1:n,[2:n 1],Fx); % assembling the discrete generator[v,d] = eig(full(A’)); % spectrum[d,I] = sort(diag(d));f_n = @(x) abs(v(ceil(max(x*n,1)),I(1))); % plot approximate inv. densityL1 = quadl(f_n,0,1); fplot(@(x)f_n(x)/L1,[0,1])
Here, only one evaluation of the vector field per interval is needed as in one dimension theintegration on the boundaries reduces to the evaluation of a single boundary point. On a partitionwith n intervals, this method then yields an accuracy of O ( n − ). Note that from Corollary 4.5 itfollows that the vector v with v i = 1 /F ( x i ) is a left eigenvector of the transition matrix (17) forthe generator at the eigenvalue 0. This fact proves pointwise convergence of the invariant densityof the discretization towards the real one. Thus Ulam’s method for the generator is very accurateand stable; the L error is O (1 /n ) as this is the rate at which approximants created from a basisof characteristic functions converge in L to regular functions.15 pectral collocation F = @(x) sin(4*pi*x)+1.1; % definition of the vector fieldn = 32; Fx = F(0:1/n:1-1/n); % evaluation in collocation nodesE = ifft(eye(n)); % basisFE = fft(diag(Fx)*E); % multiplication by vector fieldI = [0:n/2-1 -n/2:-1]; % frequenciesD = (2i*pi)*I’; % differentiation matrixA = -D*ones(1,n).*FE; % discrete generator in frequency space[v,lambda] = eig(A); % spectrumf_n = @(x) real(exp(x’*I*2i*pi)*v(:,end)); % approximate densityL1 = quadl(f_n,0,1); fplot(@(x)f_n(x)/L1,[0,1]) % plot
Here, the vector field is evaluated once per grid point (cf. the second line of the code). Aspredicted by Theorem 4.8, the accuracy increases exponentially with n (cf. Figure 2). In Figure 2 (left) we compare the efficiency of the four methods in terms of how the L -error ofthe computed invariant density depends on the number of evaluations of the vector field. The L -error was computed using an adaptive Lobatto quadrature as implemented in, e.g., Matlab’s quadl command. These errors are due to the approximation space, the approximate numericalcomputation of matrix entries, and due to solving the discretized eigenvalue problem. To illus-trate how these effects accumulate, in Figure 2 (right) we compare the errors of three differentapproximations, for different (uniform) partitions: the approximate invariant densities obtainedfrom the two Ulam type methods, and the projection of the true invariant density. We concludethat the error due to the approximation space dominates the total error. −15 −10 −5 L − e rr o r histogramUlam’s methodUlam for the generatorspectral collocation −3 −2 −1 L − e rr o r projectionUlam’s methodUlam for the generator Figure 2: Left: L -error of the approximate invariant density as a function of the number ofevaluations of the vector field (2 - 2 partition elements, resp. collocation points, 2 - 2 numberof iterations). Right: L -error of the approximate invariant densities obtained by projection ofthe true invariant density on the underlying partition, by Ulam’s method, and by Ulam’s methodfor the generator.The number of grid sets (32) chosen in Figure 1 is tiny, and chosen merely for illustrationpurposes. Likewise, the maximum number of evaluations used in both generator schemes of 512is also obviously tiny, and in practice one could very cheaply increase the number of grid sets andconcomitantly the number of evaluations. The main message from this example is that in the rightsetting: low dimension, periodic domain, and infinitely smooth vector field, spectral collocationcan significantly outperform standard Ulam and generator Ulam. Ulam’s method for the generatoris clearly outperforming standard Ulam in this example, and we will see that it continues to do soacross a variety of dimensions, domains, and vector fields.16able 1: Sources of computational cost. The dimension of the approximation space is denotedby n , d denotes the dimension of state space. In contrast to the other two methods, spectralcollocation produces a full matrix. One can reduce the cost for solving the eigenproblem in thiscase to O ( n log( n )) · approximationerror flops to set up matrix /time integration dimensionof spatialintegrals flops to solve EVP /vector iteration type Ulam’s method O ( n − /d ) O ( n ) / yes d O ( n ) / forw. Ulam’s methodfor the generator O ( n − /d ) O (cid:0) n (2 d − /d (cid:1) / no d − O ( n ) / backw. spectral collocationfor the generator O ( n − k/d ) , k ∈ N O ( n ) / no 0 O ( n ) / backw. We consider an area-preserving flow on the cylinder, defined by interpolating a numerically givenvector field as shown in Figure 3, which is a snapshot from a quasi-geostrophic flow, cf. [28]. Thedomain is periodic with respect to the x coordinate and the field is zero at the boundaries y = 0and y = 8 · . Again, we apply the three methods discussed in Section 4 in order to computeapproximate eigenfunctions of the transfer operator and the generator. x y Figure 3: Vector field of the area-preserving cylinder flow. Shown are also two saddle equilibria(black dots) together with their stable and unstable manifolds.
Since the flow is area-preserving, we have a continuum of closed orbits, and each of them is aninvariant set. On each, an invariant measure can be supported, i.e. the eigenspace at 1 (resp. 0)is infinite-dimensional. Both Ulam type methods can be interpreted as a small random perturba-tion of the original system (albeit different ones), yielding a unique stationary eigenvector of theresulting matrix. However, simultaneously with decreasing the cell size, the size of this “built-in”perturbation shrinks and it is therefore not immediate to which invariant measure the discretiza-tion schemes converge to. In order to make the problem well posed and the results comparableamong the three methods we therefore need to artificially add a random perturbation (in form ofa diffusion term) to the model.We therefore choose the noise level ε such that the resulting diffusion coefficient ε / ≈
120 (see also [18]), we choose ε = √ ·
500 here.In
Ulam’s method , for the simulation of the SDE (3) a fourth order Runge-Kutta methodis used, where in every time step a properly scaled (by a factor √ h · ε , where h is the time17tep) normally distributed random number is added. We use 1000 sample points per box and theintegration time T = 5 · , which is realized by 20 steps of the Runge–Kutta method. Note thatthe integrator does not know that the flow lines should not cross the lower and upper boundariesof the state space. Points that leave phase space are projected back along the y axis into thenext boundary box. An adaptive step–size control could resolve this problem, however at thecost of even more right hand side evaluations. The domain is partitioned into 128 ×
128 boxes.For
Ulam’s method for the generator , again, we employ a partition of 128 ×
128 boxes andapproximate the edge integrals by the trapezoidal rule using three nodes (which turns out to besufficiently accurate in this case). For spectral collocation , we employ 51 Fourier modes inthe x coordinate (periodic boundary conditions) and the first 51 Chebyshev polynomials in the y coordinate, together with Neumann boundary conditions. Due to spectral convergence, thissmaller number of modes already yields an accuracy in the invariant density comparable to theone of the other two approaches. In Table 2 we compare approximations for the three leading eigenvalues from our three schemes. Ifwe assume that spectral collocation gives a very accurate result, Ulam’s method should be O ( n − )from it. With n = 128 this matches well with the numbers in the table. Ulam’s method for thegenerator generates additional numerical diffusion; this scales the corresponding eigenvalues. InTable 2: Approximate eigenvalues.method λ λ λ Ulam’s method (log( λ i ) /T ) − . · − − . · − − . · − Ulam’s method for the generator − . · − − . · − − . · − spectral collocation for the generator − . · − − . · − − . · − Figure 4 we compare the approximate eigenvectors of the second, third and fourth relevant eigen-values of the transfer operator (resp. generator). Clearly, they all give the same qualitative picture.Yet, the number of evaluations of the vector field (and thus the associated cpu times) differ sig-nificantly, as shown in Table 3. The regions enclosed by the stable/unstable manifolds of the twoTable 3: Number of evaluations of the vector field and cpu times (on a 2.5 GHz Core Duo) forcomputing the approximate operator resp. generator and cpu times for computing the leading 4eigenvalues/-vectors.method eigs
Ulam’s method ≈ · ≈ ·
143 sec. 0.8 sec.spectral collocation for the generator ≈ · a) Ulam’s method(b) Ulam’s method for the generator(c) Spectral collocation for the generator Figure 4: Quasi-geostrophic flow: Eigenvectors at the second, third and fourth eigenvalue (fromleft to right). 19igure 5: Comparison of the approximate eigenvector at λ by Ulam’s method (left) with the onecomputed by Ulam’s method for the generator (right, together with the two saddle equilibria andtheir invariant manifolds).Ulam method we simulate a finite number of random trajectories, which incompletely sample thetrue noise distribution. In this section we demonstrate the efficacy of our generator approach for a three-dimensional flow.We consider the so-called ABC-flow [2], given by˙ x = a sin(2 πz ) + c cos(2 πy )˙ y = b sin(2 πx ) + a cos(2 πz )˙ z = c sin(2 πy ) + b cos(2 πx ) , on the 3 dimensional torus. The flow is volume-preserving, so in order to make the computationof the spectrum well posed we need to add some artifical diffusion again. The vector field is verysmooth, and in numerical experiments a diffusion coefficient ε = 0 .
04 turned out to work well.For a = √ b = √ c = 1, there is numerical evidence that the flow exhibits complicateddynamics and invariant sets of complicated geometry [11, 16].We approximate the leading few eigenfunctions of A ε with Ulam’s method for the generatoron a 64 × ×
64 grid. We used Gauss quadrature with 16 nodes on each face of a box in orderto approximate the surface integrals, requiring ∼ · ≈ · evaluations of the vector field.The L -error of the approximate invariant density is ≈ − . The second and fifth eigenfunctionsare shown in Figure 6. The regions highlighted in red and blue are invariant cylinders of thedeterministic flow; see [16] for more details on how to extract invariant and almost-invariant setsfrom the eigenfunctions.We also apply spectral collocation on a 11 × ×
11 grid, requiring 11 = 1331 evaluationsof the vector field. The L -error of the approximate invariant density is ≈ − . We observefast convergence (after comparing with computations on coarser resolutions) of not just the firstnontrivial eigenfunctions, but also of higher order ones.20igure 6: Second and fifth eigenvectors of the infinitesimal generator – Ulam type discretization.Left: slices, right: regions, where the absolute value of the left eigenfunctions exceed a giventhreshold – the “core” of the almost invariant sets, representing invariant cylinders in this example. As a last example we consider a system where the effective dynamics is supported on a set ofcomplicated geometry and fractional dimension, the well-known Lorenz system˙ x = σ ( y − x ) , ˙ y = x ( (cid:37) − z ) − y, ˙ z = xy − βz, where we make the standard parameter choices σ = 10, (cid:37) = 28 and β = 8 /
3. Since the effective(complicated) dynamics happens on a lower dimensional attractor we do not expect spectralcollocation to work well and we therefore apply Ulam’s method for the generator here.A decade ago, numerical techniques were introduced for computing box coverings for attractorsof complicated structure, cf. [7, 6]. These techniques make use of the fact that the set M to becomputed is an attractor, hence each trajectory starting in its vicinity will be pulled to M in afairly short time. In our approach time is not considered, we only use the vector field. Sincethe boundary of the box covering does not have to coincide with the boundary of the attractor,a tight box covering might not show the desired results, because of relatively big outflow ratesin boundary boxes. The simplest idea is to use a big enough rectangle – in our case [ − , × [ − , × [ − , × ×
128 boxes and used Gauss quadraturein 5 × u : We cut off at a threshold value 5 · − suchthat 96% of the invariant measure is supported on { u > c } . Figure 8 shows the correspondingcovering as well as the sign structure of the second and third eigenvector.Figure 8: Ulam’s method for the generator in an example with a low-dimensional attractor: TheLorenz attractor (left) and almost invariant sets: sign structure of the second eigenvector (middle)and third eigenvector (right). We introduced the infinitesimal generator as a tool for studying long term behavior of flows withoutrequiring any trajectory integration. The (typically unit dimensional) null space of the generatoris spanned by the invariant density of the flow; the invariant density describes the asymptoticdistribution of trajectories. We showed that eigenfunctions of the generator corresponding to22igenvalues close to 0 had strong connections with almost-invariant or metastable behavior, andcould be used to identify regions in the flow from which the escape rate is very low.We proposed two new numerical approaches for the estimation of the generator; one based onUlam’s method and the other based on spectral collocation. We tested our numerical approacheson “full phase space” flows in one, two, and three dimensions and found that both approachessignificantly outperformed the standard Ulam approach based on the Perron-Frobenius operator.Moreover, all operator/generator based methods were superior to standard trajectory integra-tion for estimating the invariant density; the identification of almost-invariant sets is essentiallyimpossible without an operator/generator approach. The spectral collocation approach workedparticularly well in periodic domains when the vector field is infinitely smooth.
Acknowledgements
We thank Yuri Latushkin for pointing us to useful books, and Nick Trefethen for careful proof-reading and helpful questions and suggestions.
References [1]
H. Amann , Dual semigroups and second order linear elliptic boundary value problems , IsraelJournal of Mathematics, 45. (1983), pp. 225–254.[2]
V. Arnol’d , Sur la topologie des ´ecoulements stationnaires des fluides parfaits , C. R. Acad.Sci. Paris, 261 (1965), pp. 17–20.[3]
A. Berman and R.J. Plemmons , Nonnegative matrices in the mathematical sciences ,SIAM, 1994.[4]
J.P. Boyd , Chebyshev and Fourier Spectral Methods , Dover Publications, 2. ed., 2001.[5]
C. Canuto, M.V. Hussaini, A. Quarteroni, and T.A. Zang , Spectral methods in fluiddynamics , Springer Series in Computational Physics, Springer-Verlag, 3. ed., 1991.[6]
M. Dellnitz, G. Froyland, and O. Junge , The algorithms behind GAIO–set orientednumerical methods for dynamical systems , in Ergodic theory, analysis, and efficient simulationof dynamical systems, Bernold Fiedler, ed., Springer, Berlin, 2001, pp. 145–174, 805–807.[7]
M. Dellnitz and A,?Hohmann , A subdivision algorithm for the computation of unstablemanifolds and global attractors , Numer. Math., 75 (1997), pp. 293–317.[8]
M. Dellnitz and O. Junge , On the approximation of complicated dynamical behavior ,SIAM J. Numer. Anal., 36 (1999), pp. 491–515.[9]
M. Dellnitz, O. Junge, W.S. Koon, F. Lekien, M.W. Lo, J.E. Marsden, K. Pad-berg, R. Preis, S.D. Ross, and B. Thiere , Transport in dynamical astronomy and multi-body problems , Int. J. Bifurcation Chaos Appl. Sci. Eng., 15 (2005), pp. 699–727.[10]
J. Ding, Q. Du, and T. Y. Li , High order approximation of the Frobenius-Perron operator ,Appl. Math. Comp., 53 (1993), pp. 151–171.[11]
T. Dombre, U. Frisch, M. Henon, J. M. Greene, and A. M. Soward , Chaotic stream-lines in the ABC flows , Journal of Fluid Mechanics, 167 (1986), pp. 353–391.[12]
G. Froyland , Finite approximation of Sinai-Bowen-Ruelle measures for Anosov systems intwo dimensions , Random and Computational Dynamics, 3 (1995), pp. 251–264.[13]
G. Froyland , Extracting dynamical behaviour via Markov models , in Nonlinear dynamicsand statistics, Alistair I. Mees, ed., Birkhauser, 2001, pp. 281–321.2314]
G. Froyland , Statistically optimal almost-invariant sets , Physica D, 200 (2005), pp. 205–219.[15]
G. Froyland and M. Dellnitz , Detecting and locating near-optimal almost-invariant setsand cycles , SIAM J. Sc. Comp., 24 (2003), pp. 1839–1863.[16]
G. Froyland and K. Padberg , Almost-invariant sets and invariant manifolds – connectingprobabilistic and geometric descriptions of coherent structures in flows. , Physica D, 238 (2009),pp. 1507–1523.[17]
G. Froyland and O. Stancevic , Escape rates and Perron-Frobenius operators: Open andclosed dynamical systems , Disc. Cont. Dyn. Sys. –Series B, 14 (2010), pp. 457–472.[18]
P. Koltai , Efficient approximation methods for the global long-term behavior of dynamicalsystems – Theory, algorithms and examples , PhD thesis, Technische Universit¨at M¨unchen,2010.[19]
D. Kr¨oner , Numerical Schemes for Conservation Laws , Wiley & Teubner, 1997.[20]
A. Lasota and M.C. Mackey , Chaos, fractals, and noise , vol. 97 of Applied MathematicalSciences, Springer-Verlag, New York, second ed., 1994.[21]
R.J. LeVeque , Finite Volume Methods for Hyperbolic Problems , Cambridge UniversityPress, 2002.[22]
A. Lunardi , Analytic Semigroups and Optimal Regularity in Parabolic Problems , Birkh¨auser,1995.[23]
A. Pazy , Semigroups of linear operators and applications to partial differential equations ,vol. 44 of Applied Mathematical Sciences, Springer-Verlag, New York, 1983.[24]
G. Pianigiani and J.A. Yorke , Expanding maps on sets which are almost invariant: decayand chaos , Transactions of the American Mathematical Society, (1979), pp. 351–366.[25]
A. Quarteroni and A. Valli , Numerical Approximation of Partial Differential Equations ,no. 23 in Springer Series in Computational Mathematics, Springer, 1994.[26]
O. Stancevic , Transfer operator methods in continuous time dynamical systems . Honoursthesis, University of New South Wales, 2007.[27]
L.N. Trefethen , Spectral methods in MATLAB , vol. 10 of Software, Environments, andTools, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2000.[28]
A.M. Treguier and R.L. Panetta , Multiple zonal jets in a quasigeostrophic model of theAntarctic Circumpolar Current , J. Physical Oceanography, 24 (1994), pp. 2263–2277.[29]
L.-S. Young , What are SRB measures, and which dynamical systems have them? , J. Statist.Phys., 108 (2002), pp. 733–754.[30]