side-by-side cylinder wake
by
Jimmy R¨ ojsel
June 2017 Technical Reports Royal Institute of Technology
Department of Mechanics SE-100 44 Stockholm, Sweden
Master’s Programme in Engineering Mechanics (120 ECTS)
Linn´e FLOW Centre, KTH Mechanics, Royal Institute of Technology SE-100 44 Stockholm, Sweden
Supervisor and examiner: Dr. Shervin Bagheri Jimmy R¨c ojsel, June 2017
side-by-side cylinder wake
Jimmy R¨ojsel
Linn´e FLOW Centre, KTH Mechanics, Royal Institute of Technology SE-100 44 Stockholm, Sweden
Abstract
In many situations, fluid flows can exhibit a wide range of temporal and spatial phenomena. It has become common to extract physically important features, called modes, as a first step in the analysis of flows with high complexity.
One of the most prominent modal analysis techniques in the context of fluid dynamics is Proper Orthogonal Decomposition (POD), which enables extraction of energetically coherent structures present in the flow field. This method does, however, suffer from the lack of connection with the mathematical theory of dynamical systems and its utility in the analysis of arbitrarily complex flows might therefore be limited. In the present work, we instead consider application of the Koopman Mode Decomposition (KMD), which is an approach based on spectral decomposition of the Koopman operator. This technique is employed for modal analysis of the incompressible, two-dimensional flow past two side- by-side cylinders at Re = 60 and with a non-dimensional cylinder gap spacing g∗= 1. This particular configuration yields a wake flow which exhibits in-phase vortex shedding during finite time, while later transforming into the so-called flip-flopping phenomena, which is characterised by a slow, periodic switching of the gap flow direction during O(10) vortex shedding cycles. The KMD approach yields modal structures which, in contrary to POD, are associated with specific oscillation frequencies. Specifically, these structures are here vorticity modes.
By studying these modes, we are able to extract the flow components which are responsible for the flip-flop phenomenon. In particular, it is found that the flip-flop instability is mainly driven by three different modal structures, oscillating with Strouhal frequencies St1= 0.023, St2= 0.121 and St3= 0.144, where it is noted that St3 = St1+ St2. In addition, we study the in-phase vortex shedding regime, as well as the transient regime connecting the two states of the flow. The study of the in-phase vortex shedding reveals—not surprisingly—the presence of a single fundamental frequency, while the study of the transient reveals a Koopman spectrum which might indicate the existence of a bifurcation in the phase space of the flow field; this idea has been proposed before in Carini et al. (2015b). We conclude that the KMD offers a powerful framework for analysis of this flow case, and its range of applications might soon include even more complex flows.
Keywords: Bluff body wake flow · Side-by-side cylinder wake · Flip-flop instability · Modal analysis · Koopman operator · Koopman mode decomposition
iii
I would first like to express my very profound gratitude to my supervisor, Dr.
Shervin Bagheri, for many enlightening discussions during the entire course of this project, and for steering me in the right direction whenever it was needed.
I would also like to thank my friends and family here in Stockholm as well as in Dalarna for constant support throughout the entire process. Thank you.
Stockholm, June 2017 Jimmy R¨ojsel
iv
Abstract iii
Acknowledgements iv
Chapter 1. Introduction 1
1.1. Modal description of fluid flows 1
1.2. Thesis overview 4
Chapter 2. Theoretical background 5
2.1. The side-by-side cylinder wake 5
2.2. The Koopman operator 7
Chapter 3. Computational framework 17
3.1. Flow simulation 17
3.2. Modal analysis 19
Chapter 4. Results 29
4.1. Flow simulation results 29
4.2. Modal analysis results 36
Chapter 5. Conclusions and outlook 47
Bibliography 51
v
Introduction
In the year 1687, Sir Isaac Newton introduced the idea of modeling the motion of physical systems with equations in his pioneering work, Philosophiæ Naturalis Principia Mathematica. Indeed, today we know that the time evolution of many natural systems can be described by a set of differential equations: one example of a class of such systems, which may exhibit very complex phenomena over a wide range of scales in both space and time, is the evolution of fluid flow, governed by the Navier–Stokes equations. Although this set of equations seems to describe the motion of fluids to near perfect accuracy in many cases, theoretical understanding of the mathematical properties of solutions in three spatial dimensions is still one of the greatest unsolved mysteries in physics.1
To deal with such an intriguing problem, various methods have been devel- oped for fluid flow analysis. It has become common practise to employ modal decomposition techniques in the pursuit of extracting information regarding the
dynamics of any flow configuration of interest.
1.1. Modal description of fluid flows
Traditionally, one of the most prominent methods for analysing fluid flows is Proper Orthogonal Decomposition (POD), which relies on the Karhunen–Lo`eve theorem (Lo`eve 1945). The method was introduced to the field of fluid dynamics by Lumley (1970) and was later further developed to incorporate a snapshot based framework by Sirovich (1987) for efficient analysis of large datasets. In applications of the original POD method in fluid flow analysis, one starts out with a vector field q(x, y, z, t), e.g. the velocity field, and assumes it can be decomposed as
q(x, y, z, t) =X
j
aj(t)φj(x, y, z), (1.1) where the aj(t):s are time-dependent scalar expansion coefficients and φ(x, y, z) are the POD modes, defined on a spatial domain with coordinates x, y, z.
1The Millennium Prize Problems are seven problems in the field of mathematics that were stated by the Clay Mathematics Institute in 2000—one of which is global regularity problem for the Navier–Stokes equations. If a correct solution to any of these problems is presented, the discoverer(s) will be awarded a $1 million dollar prize. Website: www.claymath.org
1
To compute the POD modes in practise, one starts out with a dataset of snap- shots {xk}Nk=1, where each xk ∈ Rn. These xk:s corresponds to measurements of the field variable q, sampled at different times, with equal timestep between each measurement, and with the time-averaged mean of the entire dataset subtracted. Typically one chooses to sample measurements of the velocity field (associated with a particular flow configuration). The measurements can be obtained from either experiments or numerical simulations. Upon constructing the so-called correlation matrix R ≡ XXT ∈ Rn×n, where the columns of the matrix X ∈ Rn×N are the measurements xk, one solves the eigenvalue problem
Rφj = λjφj, λ1≥ . . . ≥ λn ≥ 0.
The eigenvectors φjfound from the above eigenvalue problem are then identified as the POD modes, which are mutually orthonormal vectors representing the field variable data optimally in the sense thatPr
j=1kPxj− xjk22 is minimized for each r (Lo`eve 1945); here P denotes projection of the data onto the subspace spanned by the r first POD modes. In this sense, the POD modes represent energetically coherent structures in the flow field and the study of such structures may be useful to characterise the dynamics of the flow case under consideration.
For the ease of presentation, we here suppose that the measurements have been captured on a uniform grid; in any other case the expression for the correlation matrix would read R ≡ XXTW , where the entries of W represents spatial weights. Although fluid systems are typically high-dimensional, the dynamics may be restricted to a low-dimensional subspace. In general, one only retains r number of modes to represent the field variable q, where r is determined by requiring
r
X
j=1
λj/
n
X
j=1
λj≈ 1.
This is an instant of what is known as model-reduction: the high-dimensional flow field is represented by a small number of spatial modes. The discretized POD coefficients aj ∈ RN, corresponding to modes φj, are obtained by considering (1.1) in a discrete setting; one can deduce that
[a1, . . . , ar]T = [φ1, . . . , φr]TX.
In applications where the measurements are flow simulation data, the velocity field may have been discretized on a computational grid with millions of grid points. In such cases, we typically have N n. Therefore, the method of snapshots (Sirovich 1987) may be much more efficient: instead of solving the above eigenvalue problem, one considers
XTXϕj= λjϕj,
where XTX ∈ RN ×N is typically much smaller than XXT ∈ Rn×n. The nonzero eigenvalues of XXT and XTX are equivalent, and the corresponding
a) b)
c) d )
e) f )
Figure 1.1. Modal decomposition of the flow past a circular cylinder at Re = 100. Colors represent vorticity, with red as positive vorticity and blue as negative vorticity. a): Flow snapshot. b): Mean flow. c), d ): First and second POD mode.
e), f ): Real part of fundamental oscillating Koopman mode, shown at phases φ = 0 and φ = π, respectively.
POD modes may be recovered through φj= Xϕj
1 pλj
.
The snapshot POD method thus offers an intriguing framework for modal analysis and model reduction of fluid flow data; we refer to Semeraro (2013) for a further exposition on modal analysis of high-dimensional fluid simulations.
POD does, however, lack connection with the strong mathematical foundation of dynamical systems. Keeping this in mind, POD cannot be trusted to provide a reliable framework for analysis and model reduction of flow configurations exhibiting arbitrarily complex dynamics.
The Koopman Mode Decomposition (KMD) is another technique for data analysis, which has caught a great deal of attention within the fluids community during recent years. Here, the time evolution of observables, i.e., time-dependent measurements of the flow field, is connected to the evolution of trajectories in the infinite-dimensional state space of the flow, and the time evolution of each observable may be represented as a linear expansion of invariants of the so-called Koopman operator. One obtains a flow field representation with two essential components: complex-valued Koopman modes which are the global structures in the flow field that grow, decay or oscillate in time and Koopman eigenvalues
which correspond to the time scales of the growth (or decay) and oscillations of Koopman modes. This procedure thus allows for extraction of coherent structures in fluid flows which, in contrary to POD modes, are associated with specific oscillation frequencies. A variety of algorithms for computation of Koopman invariants have been proposed in the literature; most of these methods are based on a snapshot type framework, much like the snapshot POD method described above, and thus they do not require the knowledge of governing equations—only snapshot data of some field variable.
In some cases, the POD and Koopman modes may appear to be quite similar, in spite of the different underlying theoretical framework. This occurs in particular for periodic flows, such as the incompressible, two-dimensional flow past a circular cylinder at Re = 100, which exhibits the famous von K`arm`an vortex street instability; see Figure 1.1. Because there is only one frequency present in this flow (that of the vortex shedding), one should expect that there are both POD and Koopman modes associated with this frequency. Indeed, such modes exist, and the agreement between the POD and Koopman modes is striking. It is, however, important to note that if the input data is real- valued, the POD modes are real-valued. If the modes correspond to oscillating structures, a pair of POD-modes are thus required to represent these structures, while the Koopman modes are complex-valued, and may be phase-shifted to exhibit a different real part. A single Koopman mode, along with the mean flow, is thus sufficient to represent the dynamics. For thorough POD and Koopman mode based investigations of the wake flow past a single circular cylinder at low Reynolds numbers, we refer the reader to Noack et al. (2003) and Bagheri (2013), respectively.
Other than POD and KMD, there exists a range of other, less common, methods for modal analysis of fluid flows; see Taira et al. (2017) for a recent review. In the present work, we employ the Extended Dynamic Mode Decompo- sition (EDMD) algorithm (Williams et al. 2015) to obtain the KMD. Thismethod offers a powerful framework for accurately computing Koopman invariants. We consider application to the two-dimensional flow past two side-by-side circular cylinders at low Reynolds numbers, which is known to exhibit surprisingly complex dynamics in spite of the simple configuration. Specifically, this flow can exhibit the so-called flip-flopping phenomena. This flow is characterised by quasi-periodic vortex shedding, and thus contains multiple frequencies. Hence, the KMD is particularly well-suited for analysis of this flow case.
1.2. Thesis overview
The thesis is organised as follows. In Chapter 2, we review the theoretical back- ground of the present investigation, with emphasis on the Koopman operator.
Chapter 3 give details on the computational framework for the flow simulation, as well as a review on various algorithms for computation of Koopman invariants.
In Chapter 4, we present and discuss our results. Finally, in Chapter 5, we conclude our findings and discuss potential future investigations.
Theoretical background
The objective of this chapter is to introduce the physical and mathematical background used in the present investigation. In section 2.1, the flow case under consideration is briefly described along with historical notes. Section 2.2 introduces the Koopman operator along with some of its properties.
2.1. The side-by-side cylinder wake
Bluff-body wake interaction, such as the flow past two or more circular cylinders in various configurations, has enjoyed considerable attention within the fluids community in the past. These flows plays important roles in several industrial applications such as in the off-shore industry, or in the design of high-rise buildings. A simple, yet interesting model of these flow configurations is the incompressible, two-dimensional flow past two circular cylinders placed in a uniform free-stream. The flow is essentially governed by two non-dimensional parameters: the Reynolds number, Re = U∞D/ν, and the non-dimensional gap spacing between the cylinders, g∗ = d/D. Here U∞ is the free-stream velocity, D is the cylinder diameter, ν is the kinematic viscosity, and d is the (dimensional) gap spacing between the two cylinders.
As the gap spacing g∗and the Reynolds number Re are varied, the wake flow has been observed to change significantly in past investigations; Zdravkovich (1977) summarized the state of knowledge up to that point. Kang (2003) conducted a comprehensive numerical parameter study of this flow case in the laminar regime; here six different wake patterns were identified: in-phase synchronized (g∗≤ 1.5), anti-phase synchronized (g∗≥ 2), ‘single bluff body’
(g∗ ≤ 0.4), deflected oscillatory (50 ≤ Re ≤ 110 and 0.2 ≤ g∗ ≤ 1), steady symmetric (Re ≤ 40 and g∗≥ 0.5) and finally the so-called flip-flopping wake flow for 0.4 ≤ g∗≤ 1.5. Examples of these wake patterns are shown in Figure 2.1.
In the present work, we are mainly interested in the flip-flop phenomenon.
This wake flow is characterised by an asymmetric unsteady wake in which the gap flow is deflected alternatively towards one of the cylinders with a frequency about one order of magnitude lower than the vortex shedding frequency; this phenomenon mainly arises as a secondary instability associated with the in-phase synchronized regime (Kang 2003).
5
-1.5 0 1.5
a) b)
c) d )
Figure 2.1. Wake patterns. a): ‘Single bluff body’ flow;
g∗= 0.3, Re = 60. b): Anti-phase synchronized flow; g∗= 2.5, Re = 80. c): In-phase synchronized flow; g∗= 1, Re = 50. d ):
Flip-flopping flow; g∗= 1, Re = 60. Colors represent vorticity levels, cut off at -1.5 and 1.5 for the ease of presentation.
Recently, Carini et al. (2014) discovered that the onset of the flip-flop phenomon may be attributed to to a bifurcation in the state space of the flow, which occurs for certain values of g∗ and Re. The authors performed a Floquet stability analysis based on the centre-manifold reduction technique (see Carini et al.
(2015a) for details). The method was applied to flow data aquisited as the flow exhibited the in-phase wake pattern, right before the onset of the flip-flopping wake, and it was found that a pair of complex-conjugate Floquet modes become unstable for Re > 61.74 if the gap size is fixed at g∗= 0.7. Later, Carini et al.
(2015b) conducted a parameter study and provided a stability map (with g∗ and Re as coordinates) and found that in the low Reynolds number regime, the flip-flopping wake arises for two different gap size intervals; g∗ . 1.4 and g∗ & 1.6, and argued that different physical mechanisms may be the cause of this. We have reproduced the stability map from Carini et al. (2015b), see Figure 2.2. The aforementioned investigations focused on low Reynolds numbers (Re < 110). For slightly higher Reynolds numbers, the flip-flopping has been investigated experimentally: Williamson (1985) found that the flip-flop is observed when 0.6 < g∗< 0.9 for 50 < Re < 150, and Peschard & Gal (1996) concluded that the flip-flop is non-chaotic but quasiperiodic for 1.6 < g∗< 3.2 and 110 < Re < 150. Interestingly, the flip-flop is found at much higher number as well; Zdravkovich & Pridden (1977) considered the turbulent regime at Re = O(105) in an experimental investigation and observed that the gap flow switched direction at irregular time intervals. More recently, Afgan et al. (2011) performed a large-eddy simulation for Re = 3900 and observed that also in this regime the gap flow switches direction intermittently.
present study
In-phase
Anti-phase
Flip-flop Flip-flop
Figure 2.2. Neutral stability curves: the different areas rep- resent various vortex shedding regimes for the side-by-side cylinder wake. Reproduced from Carini et al. (2015b).
In the present work, we shall consider the case g∗ = 1 and Re = 60, which, according to Carini et al. (2015b), results in an in-phase vortex shedding wake followed by the flip-flopping wake flow, as shown in Figure 2.1 c) and d ), respectively. Specifically, we shall perform Koopman mode analysis based on flow simulation data to characterise the pre-flip-flop regime, the transient regime, and the post-transient flip-flopping regime, respectively.
2.2. The Koopman operator
Traditionally, analysis of dynamical systems rely on a geometrical picture which focuses on “dynamics of states”, which means that the trajectories of the system in the state space are analysed essentially by means of basic notions in geometry.
This framework was introduced by Henri Poincar´e in the 19th century as he was working on the famous three-body problem in classical mechanics, by employing a concept we today know as Poincar´e maps. This approach has been proven useful in the analysis of low-dimensional dynamical systems (Mainieri et al.
2015). This viewpoint is, however, clearly inherently ill-suited for analysis of high-dimensional systems, such as the evolution of fluid flow. To overcome this issue, data-driven, operator-based techniques have been developed. One of the protagonists in the field of fluid dynamics is the so-called Koopman operator, which governs the evolution of observables: functions defined on the state space on which the dynamical system of interest evolves. The Koopman operator was discovered in the 1930’s by Bernard Koopman, who presented this operator-theoretical approach to describe the evolution of Hamiltonian (measure-preserving) flows in the state space (Koopman 1931; Koopman & von
Neumann 1932).
The Koopman operator is linear, and in general infinite-dimensional. Because of its linearity, it admits a spectral decomposition, even for dissipative systems, which was proven only recently in Mezi´c & Banaszuk (2004) and Mezi´c (2005), who obtained the Koopman mode decomposition by applying harmonic averaging.
Analogous to e.g. Fourier analysis, it is possible to extract coherent structures in a dynamical system, oscillating with their respective frequencies, by using the Koopman operator theory. When considering linear systems, Koopman analysis essentially reduces to global stability analysis (Chen et al. 2012).
2.2.1. General setting
To properly boil down the Koopman theory to applications in fluid dynamics, we shall start out in a more general setting. To this end, we consider a possibly nonlinear dynamical system
∂
∂tx = f (x), (2.1)
where we identify the state x with a point on a state space M ⊂ Rn, and f is, formally, a smooth vector field. For any given time increment τ > 0, we may consider the flow map Fτ: M → M, which maps the state x(t) forward into the future:
Fτ(x(t)) = x(t + τ ) = x(t) + Z t+τ
t
f (x(η))dη.
This induces the discrete-time dynamical system
xk+1= Fτ(xk), (2.2)
where xk ≡ x(kτ ). Historically, discrete-time dynamical systems are more general than continuous time systems, and was considered by Koopman (1931), but we start with continuous time for illustrative purposes. We shall, in addition, also introduce the notion of real-valued observables g(x) : M → R, which are elements of the infinite-dimensional Hilbert space F . In the context of Koopman analysis, this Hilbert space is typically given by the Lebesque square-integrable functions on M; that is, F = L2(M). Other measure spaces are also valid (Budiˆsi´c et al. 2012), but we may, and shall, be far more restrictive in our analysis. For a given τ > 0, the associated Koopman operator Kτ is an infinite-dimensional linear operator that acts on observable functions g as:
Kτg = g ◦ Fτ where ◦ denotes composition. That is,
Kτg(xk) = g(Fτ(xk)) = g(xk+1).
The Koopman operator Kτthus defines an infinite-dimensional linear dynamical system that advances the observable gk= g(xk) to the next timestep:
gk+1= Kτgk. (2.3)
Note that (2.3) holds for any observable g and for any xk∈ M, why it makes sense that the Koopman operator is infinite-dimensional. Furthermore, (2.3) holds for any τ > 0, so whenever we speak of the Koopman operator, we have
implicitly assumed that there is a corresponding time increment τ . In Koopman (1931), measure-preserving fluid systems were investigated. In such cases, the Koopman operator Kτ is unitary, and forms a one-parameter family of unitary transformations in the Hilbert space F , parametrized by τ . In the context of fluid dynamics, a measure would define a notion of volume in the (infinite- dimensional) state space of all possible velocity fields, and measure-preserving dynamics would preserve this notion of volume. In general, this is the case when dynamics are given by the Euler equations for incompressible fluid flows.
Here, however, we do not restrict ourselves to measure-preserving dynamical systems, as we shall consider the (dissipative) Navier–Stokes equations.
It is possible to describe the continuous-time version of the Koopman dynamical system (2.3) with the infinitesimal generator K of the one-parameter Koopman group Kτ:
d
dtg = Kg. (2.4)
The linear dynamical systems (2.3) and (2.4) are analogous to the dynamical systems in (2.1) and (2.2), respectively. It is important to note that the above reasoning also holds for vector-valued observables g = (g1, . . . , gm), and thus even the state x itself may be the observable, and the infinite-dimensional operator Kτ will thus advance the state.
An eigenfunction of the Koopman operator is a complex-valued observable φ(x) ∈ F which satisfies
Kτφ(x) = exp(λτ )φ(x). (2.5)
In what follows, we shall refer to λ ∈ C as the continuous-time Koopman eigenvalue (corresponding to φ(x)), while µ ≡ exp(λτ ) is refered to as the corresponding discrete-time eigenvalue. The eigenfunctions may be interpreted as “coordinates” along which the evolution of the observable g is linear. Given a vector-valued observable g : M → Rm, the Koopman mode, vj, corresponding to λj is the vector of expansion coefficents required to represent g in a basis spanned by the eigenfunctions, and these coefficients may be complex numbers.
In what follows, we shall refer to the tuple (λ, φ, v) of Koopman invariants as the Koopman tuple. Following Mezi´c (2005), the Koopman mode decomposition (KMD) of any smooth and bounded1 observable may be written as
g(xk) =
∞
X
j=0
φj(x)vj+ R, (2.6)
where R represents any contribution from continuous or residual parts of the Koopman spectrum. If the dynamical system measure preserving, then the eigenfunctions form an orthogonal basis of F , which means that R vanishes; see Mezi´c (2005) for further details. In what follows, we will only be interested in the point, or singular, part of the spectrum since we will aim to compute the
1Because our ultimate goal is analysis of a fluid flow continuum, these assumptions on the nature of the observables are clearly justified. Note, however, that the Koopman theory generalizes to e.g. state spaces without structure (Budiˆsi´c et al. 2012).
x1 x2
x3
xN
xN
Fτ
M
g
y1
y2
y3
yN Rm
Kτ
Figure 2.3. The action of the Koopman operator is equivalent to evolving the nonlinear dynamics. A similar cartoon can be found in Brunton et al. (2016).
Koopman tuples by means of (inherently discrete) data-driven methods, and we thus neglect the residual term R for now. Thus, by (2.6), the interpretation of the Koopman modes vj from a physical perspective is spatial structures that oscillates with frequencies ωj = Im{λj} and grows exponentially with σj = Re{λj}. Note that j = 0 corresponds to the (ergodic) time-averaged mean of the observable. Let τ > 0 be given. Then, by combining (2.3) and (2.6), and assuming all elements of g lie in the span of the eigenfunctions, we have
g(xk+1) = Kτg(xk) =
∞
X
j=0
exp(λjτ )φj(xk)vj. (2.7) From (2.7) it is clear that we may either choose to evolve the dynamics through our nonlinear mapping or through the Koopman operator. To further elucidate this, let yk ≡ g(xk). Then, evolving the state xk through Fτ is equivalent to evolving yk through Kτ. Additionally, if the mapping g : xk→ yk is injective, then we have determined coordinates in which the evolution is linear; that is, the components of g are Koopman eigenfunctions. See Figure 2.3. The above expansion can be applied to fields of observables as well, such as a velocity field, or a vorticity field. In such case, the Koopman modes becomes fields of expansion coefficients, and one can think of such Koopman modes as a steady velocity or vorticity field.
2.2.2. Applied Koopmanism
To build an intuition for the Koopman mode decomposition (KMD), we shall consider two examples of applied Koopman analysis.
A linear system
Consider the linear system (Mezi´c 2013)
˙ x = Ax,
where A ∈ Rn×n. Let λj, vj and w∗j denote the eigenvalues, the right eigenvec- tors and the left eigenvectors of A, respectively, scaled so that hvi, wji = δij.
Then, the λj:s are Koopman eigenvalues with vj as the Koopman modes and φj(x) = hx, wji are Koopman eigenfunctions. This follows directly from
φ = h ˙˙ x, wi = hAx, wi = hx, A∗wi = λφ and thus
φ(x˙ 0) = Kτφ(x0) = exp(λτ )φ(x0).
Assuming A has a full set of eigenvectors with distinct eigenvalues, we have x =
n
X
j=1
hx, wjivj =
n
X
j=1
φj(x)vj
and consequently, by (2.7), we have Kτx0=
n
X
j=1
hvj, wjivj =
n
X
j=1
exp(λjτ )φj(x0)vj.
The φj:s are, however, not the only Koopman eigenfunctions, and the λj:s are not the only eigenfunctions: if φ1 and φ2 are two Koopman eigenfunctions with associated eigenvalues λ1and λ2, then φ1· φ2is an eigenfunction with eigenvalue λ1+ λ2. This holds in general (Budiˆsi´c et al. 2012) and not only for the above linear system, and follows directly from (2.5):
Kτ(φ1(x)φ2(x)) = φ1(Fτ(x)) · φ2(Fτ(x))
= exp(λ1τ )φ1(x) exp(λ2τ )φ2(x)
= exp((λ1+ λ2)τ )φ1(x)φ2(x).
Thus, by induction, one can easily show that if φ is a Koopman eigenfunction with eigenvalue λ, then φn is a Koopman eigenfunction with eigenvalue nλ if n ∈ Z+; the eigenfunctions of the Koopman operator form an Abelian group.
Burgers’ equation
An interesting example of applied Koopman analysis in the context of fluid dynamics is that of the viscous Burgers’ equation (Kutz et al. 2016b):
∂u
∂t + u∂u
∂x = ε∂2u
∂x2, ε > 0, x ∈ (−∞, ∞). (2.8) This equation models the evolution of e.g. shock waves; in such applications the solution variable u(x, t) is the fluid flow velocity at the coordinate x and time t. The evolution of the system, as illustrated in Figure 2.4, indeed gives rise to a shock wave formation, which is regularized by the viscous term on the right-hand side of Eq. (2.8). For ε > 0, the solution is continuous for all t > 0, given that the initial condition is continuous. Burgers’ equation is one of few nonlinear partial differential equations whose analytic solution can be derived, here by means of the Cole–Hopf transformation (Hopf 1950, Cole 1951). Let
u = −2ε1 v
∂v
∂x,
10 8 6
t
0 4
-10 0.5
u(x;t)
-5
x
2 1
0
5 0
10
Figure 2.4. Solution of Burgers’ equation (2.8) with initial condition u(x, 0) = exp(−x2/2). Here ε = 0.05.
then, the transformation to the new variable v(x, t) results in the linear heat equation
∂v
∂t = ε∂2v
∂x2. (2.9)
We have here assumed that ε > 0 to ensure well-posedness of (2.9). We now assume that v(x, t) can be expanded in a Fourier basis:
v(x, t) =X
k∈Z
ˆ
vk(t) exp(ikx). (2.10)
Combining (2.9) and (2.10) yields
∂ ˆvk
∂t = −εk2vˆk, k ∈ Z, with the solution
ˆ
vk= ˆv0exp(−εk2t), (2.11) where ˆv0is the Fourier transform of the inital condition. In the spatial domain, the solution is found to be
v(x, t) = exp − Rx
−∞u(ξ, t)dξ 2ε
! .
For each k, we may consider the observable g = ˆvk. By (2.11), we identify the Koopman operator as
Kt= exp(−εk2t).
For most PDE:s, it is however very difficult, if not impossible, to construct the Koopman operator analytically, why one is typically forced to resort to data-driven methods for computing Koopman invariants.
2.2.3. Koopman analysis of fluid flows
When analysing fluid flows by means of data-driven Koopman algorithms, one must bear in mind that that each experiment or numerical simulation will only provide us with a single trajectory in the state space. That is, if u0is the initial velocity state, (at time t = 0), and g(u0) is an observable of this velocity field, then the KMD in (2.7) may be expressed as
g(u(τ )) = Kτg(u0) =
∞
X
j=0
exp(λjτ )φj(u0)vj, (2.12) assuming the observable lie in the span of the eigenfunctions. For limit-cycling fluid systems, the growth rate of each mode is zero; that is Re{λ} = 0. Moreover, in the case of quasi-periodic limit cycling flow, the dynamics exhibits multi- ple frequencies. Let Ω = (ω1, . . . , ωp)T denote a vector of incommensurable frequencies. Then, (2.12) can be expressed as (Arbabi & Mezi´c 2017)
g(u(τ )) = Ktg(u0) = X
j∈Zp
exp(ij · Ωτ )φj(u0)vj, (2.13) again assuming g lie in the span of the eigenfunctions. Note that this is essentially nothing but a generalized Fourier-series expansion of the observable g. This means that the time-series of a quasi-periodic velocity field u can be expressed as a sum of the periodic signals with frequencies that are integral combinations of the frequencies in Ω. Above, we have assumed that the observables are in the span of the eigenfunctions. In particular, this is true for the observable g(u) = u under the assumption that the underlying dynamics is non-chaotic and ergodic (Arbabi & Mezi´c 2017), which is assumed to be the case for the flip-flopping cylinder wake.
In what follows, we shall employ data-driven, snapshot type algorithms for computation of Koopman tuples. To this end, we must consider the choice of input data. Since the state of the fluid flow is fully defined through the velocity field u, the natural choice would be to collect a series of snapshots of the velocity field obtained from a numerical simulation of the flow. However, we may cut the computational cost in half by performing the modal analysis on snapshots of the z-component of the associated vorticity field
ω = ∇ × u
instead, as opposed to use data corresponding to both the x and y-components of the full velocity field. The reason for this is that we would like to exhibit the vorticity of the modal structures since this is the most common approach in the literature; instead of computing the Koopman tuples using the full velocity field and then compute the vorticity of these modes, we use vorticity snapshots as input data. To use vorticity snapshots in lieu of velocity snapshots is quite
uncommon in the literature, but the following recent result from Arbabi &
Mezi´c (2017) asserts that this approach is valid.
Proposition 1 (Arbabi & Mezi´c 2017). Let g, h : Ω × M → R be two fields of observables defined on the the flow domain Ω and the state space of the flow, denoted by M. Assume the observables are bounded over their domain and related to each other through
DΩg = h,
where DΩdenotes the partial derivative with respect to a spatial coordinate in the flow domain Ω. Let vgj and vjh denote the Koopman modes of observables associated with the Koopman frequencies ωj, j = 1, 2, . . .. Then,
DΩvjg= vhj for every x, y, z ∈ Ω.
Proof. Let us define the finite-time harmonic average of g at the frequency ωj by
gjT = 1 T
Z T 0
g(x, y, z, t) exp(−iωkt)dt. (2.14) The Koopman mode of g associated with ωk is then given by (Mezi´c 2005)
vkg= lim
T →∞gkT.
The Leibniz rule implies that the spatial derivative commutes with the finite-time averaging in (2.14), i.e.,
hTk = DΩgTk
for any finite T . Let us now fix the spatial coordinates x, y, z. If both g(x, y, z, t) and h(x, y, z, t) are integrable functions, the existence of vgk(x, y, z) implies exis- tence of vkh(x, y, z), and since the scalar functions hTk(x, y, z) and DΩgTk(x, y, z) are equal up to any finite T , their limits as T → ∞ must be equal, i.e.,
vhk(x, y, z) = Dvkg(x, y, z).
The above proposition thus implies that we may use vorticity snapshots instead of velocity snapshots when computing Koopman tuples by means of data-driven methods given that the vorticity field is bounded—which is the case for incompressible, viscous two-dimensional flow with no-slip conditions (Leray 1934)—and given that the system evolves on a periodic attractor. In particular, applying KMD to vorticity snapshots yields, in principle, the same eigenvalues as KMD applied to velocity snapshots.
If a general dynamical system is evolving on an intertial manifold, such as the transient regime between the periodic in-phase vortex shedding regime and the flip-flopping limit cycle, Koopman eigenfunctions might not even exist; and the KMD expansion (2.6) reduces to the residual term R. However, it has been shown (Gaspard & Tasaki 2001) that for systems undergoing a Hopf bifurcation (Hopf 1942), in which the state of a system changes from an unstable fixed point to a periodic attractor, the Koopman operator has both a continuous spectrum as well as a point spectrum, where the (continous-time) eigenvalues form a
pyramid-shaped lattice in the stable half-plane. A typical example of a Hopf bifurcation in the context of fluid mechanics is the von K`arm`an instability in the two-dimensional flow past a circular cylinder, which arises as the Re & 47;
see Sreenivasan et al. (1987). Data-driven Koopman analysis of transient states may, however, be extremely hard; if possible at all. Assuming, however, that a point spectrum is present, the question of what observables to use arises;
data-driven Koopman analysis will only reveal the part of the spectrum which is contained within the span of the observables. We shall discuss this issue further in Section 3.2.
As a final remark in this section, we may briefly consider the KMD of an observable of a chaotic, ergodic system, such as the turbulent high Reynolds number flip-flopping cylinder wake, which is given by
g(u(τ )) = Kτg(u0) =
∞
X
j=0
exp(iωjτ )φj(u0)vj+ Z ∞
−∞
exp(iατ )˜v(α)ρ(α)dα.
The second term on the right-hand side represents the contribution of the continuous spectrum, where iα, α ∈ R denotes a continuum of eigenvalues, and the ρ(α) represents the spectral density which indiciates how strong the contribution of α is in the evolution of the observable. The above expression for the KMD is clearly suitable for turbulent flows with mixed spectra, but we shall not consider turbulence here; instead we refer to Arbabi & Mezi´c (2017) for a thorough exposition on Koopman analysis of such flows. For a review on applications of the Koopman theory in fluid flow analysis in general, we refer to Mezi´c (2013) and references therein.
Computational framework
In the following sections, we describe the computational framework. Specifically, in section 3.1, details on the direct numerical simulation (DNS) of the side- by-side cylinder wake are given. In section 3.2, we describe our approach to compute Koopman tuples by means of collecting data snapshots from the DNS.
3.1. Flow simulation
The fluid flow under consideration is governed by the two-dimensional, incom- pressible Navier–Stokes equations in non-dimensional form, given by
˙
u = −(u · ∇)u + 1
Re∇2u − ∇p + λ 0 = ∇ · u,
where u(x, y, t) is the fluid velocity, p(x, y, t) is the fluid pressure, and λ may be interpreted as a Lagrange multiplier to enforce incompressibility on solid boundaries (the cylinder walls). Moreover, Re = U∞D/ν is the Reynolds number, where U∞ is the freestream velocity, D is the cylinder diameter and ν is kinematic viscosity. Following Carini et al. (2015b), we choose a rectangular computational domain of size 125D × 100D. Specifically, both cylinders are positioned on the y-axis symmetrically about the x-axis, and spaced a distance d apart. The inflow boundary Γin is positioned a distance Lin= 50D upstream of the cylinder centers while the outflow boundary Γout is at Lout = 75D downstream of the cylinders, and the half-width of the domain is Ls= 50D.
The grid resolution is Nx= 430D times Ny = 450D, with progressive refinement towards the region [−1D, 1D] × [−2.5D, 2.5D], enclosing both cylinders, where the resolution is D/50. We impose convective boundary conditions at the outflow boundary Γout, symmetric velocity boundary conditions at Γtop and Γbottom, and uniform inflow condition (U∞, 0) at the inflow boundary Γin. See Figure 3.1 for details. The governing equations (in primitive variabels) are discretized on a staggered grid using the finite volume, second order projection-correction method; see Perot (1993) for details. We employ second-order semi-implicit time integration using a time step ∆tDNS = 1/150. Choosing U∞ = D = 1, the CFL number based on the refined part of the grid thus amounts to 1/3.
Moreover, a probe is placed at (0.5D, 0) for measuring the vertical velocity.
Since the cylinder boundaries does not conform with the grid, we employ the
17
d x y
Lin
Lout
probe
D
Ls 2
1
Γout Γin
Γbottom
Γtop
U∞
Figure 3.1. Computational domain, purposely not to scale.
Immersed Boundary Method (Taira & Colonius 2007); here the fluid pressure is projected onto Lagrangian points placed on the cylinder surfaces to enforce incompressibility and no-slip conditions at these boundaries; 96 Lagrangian points are positioned on each cylinder boundary, uniformly spaced. A natural by-product of this apporach is that we may measure the lift (CL) and drag (CD) coefficients, respectively, which provides a means of validating the grid by benchmarking the lift- and drag coefficents in the case of a single cylinder against values found in the literature. Although the present DNS code has been used in previous works (Bagheri et al. 2012; Bagheri 2013; L¯acis et al. 2016), we perform a simple comparison study for validation. To this end, we only use the top cylinder and consider the time-averaged drag coefficent ¯CL, as well as CL0; the maximum amplitude of the lift coefficient. Additionally, we consider the Strouhal number St = ω0D/2πU∞, where ω0is the vortex shedding frequency.
As observed in Table 3.1, our grid yields acceptable results.
Table 3.1. Validation of numerical method: comparison study for flow over a single cylinder (the top cylinder in Figure 3.1).
Re C¯D CL0 St
Present 40 1.52 – –
100 1.34 0.33 0.164
Kang (2003) 40 1.51 – –
100 1.33 0.32 0.165
Park et al. (1998) 40 1.51 – –
100 1.33 0.33 0.165
Williamson (1989) 100 – – 0.164
3.2. Modal analysis
As stressed in the introduction, we have chosen to work within the Koopman framework for our modal analysis of the side-by-side cylinder flow because of the extensive mathematical theory that comes with it. Computation of Koopman invariants is, however, a challenging problem and is an active area of research.
A variety of algorithms for applied Koopman analysis have been proposed in the past. Examples include harmonic averaging (Mezi´c & Banaszuk 2004;
Mezi´c 2005), generalized Laplacian analysis (Budiˆsi´c et al. 2012), diffusion maps (Giannakis 2016), Dynamic Mode Decomposition (Schmid & Sesterhenn (2008);
Schmid (2010); Tu et al. (2014)), and Extended DMD (Williams et al. 2015) and its kernel variant (Williams et al. 2016).
Because of the high-dimensional nature of fluid physics, one is typically limited to the use of data-driven methods for applied Koopman analysis. Much like the snapshot POD method as described in the introduction, DMD, Ex- tended DMD (EDMD) and Kernel DMD (KDMD) are methods to produce approximations of Koopman tuples solely by providing a collection of data snapshots, i.e., no governing equations are required. These snapshots are, again, typically measurements of the flow field obtained from either experiments or numerical simulations. DMD was introduced by Schmid & Sesterhenn (2008) and Schmid (2010), while the the actual connection to the Koopman operator was established in Rowley et al. (2009). A variety of improvements to the original DMD algorithm have been proposed during recent years. Noticeable examples are the Sparsity Promoting DMD (Jovanovi´c et al. 2014) for selection of a small subset of modes best describing the dynamics, or Total DMD by Hemati et al. (2015) for mitigating the effect of measurements noise, which may contaminate the snapshot data.
The EDMD and KDMD methods are developed with the aim of more accurately determining Koopman tuples from nonlinear data, and we shall further review these methods in what follows. At this point, it may be appro- priate to stress that, in general, one seeks to compute Koopman eigenfunctions as accurately as possible, since the Koopman modes—the oscillating spatial structures that will give us insight in the flow physics—are merely projections of observables, such as a velocity field, onto these eigenfunctions.
3.2.1. Dynamic Mode Decomposition
For consistency, we start out by reviewing the DMD algorithm, as given by Tu et al. (2014); dubbed Exact DMD. As shown in that work, this formulation of the algorithm is, under certain circumstances, equivalent with that of Schmid (2010); which in turn is arguably the most widely used Koopman algorithm in the context of fluid flow analysis up to this date. The reason why we do not consider the original algorithm is that the Exact DMD has a more obvious connection with EDMD and KDMD: while Schmid & Sesterhenn (2008) only consider a single collection {xk}Nk=1of snapshot data sampled from a dynamical system, where the xk∈ Rn, as a basis for computing Koopman parameters, Tu
Algorithm 1 Exact DMD (Tu et al. 2014) Input: Snapshot pairs {(xk, x0k)}Nk=1
Output: Modes vj and corresponding discrete-time eigenvalues µj
1: Compute the reduced SVD of X: X = U ΣVT.
2: If desired: truncate the SVD by only considering the first r columns of U and V , and the first r rows and columns of Σ, to obtain Ur, Σr and Vr.
3: Let ˜A = UrTAUr= UrTX0VrΣ−1r (an r × r matrix).
4: Find the eigenvalues µj and eigenvectors ˜vj of ˜A by solving ˜A˜vj= µjv˜j; every nonzero µj is a DMD eigenvalue.
5: The DMD mode corresponding to µj is given by vj = µ−1j X0VrΣ−1r v˜j.
et al. (2014) generalized the idea to incorporate a snapshot pair based framework.
This is an idea also used in the EDMD and KDMD methods. Suppose we have pairs of snapshots {(xk, x0k)}Nk=1 where x0k = Fτ(xk) (following our previous notation) and that there is a linear relationship between these snapshot pairs, i.e.
x0k+1= Axk.
The Schmid DMD is a special case with xk+1 = x0k, but this generalization enables the use of multiple time series. We shall, however, not pursue this possibility in our analysis. As with POD, these measurements are most often an entire velocity field—or, as in our case, an entire vorticity field, sampled from DNS data. Upon forming the matrices X and X0, the columns of which are the snapshot collections {xk}Nk=1 and {x0k}Nk=1, one forms the matrix
A = X0X†, (3.1)
where†denotes the Moore-Penrose pseudoinverse: if X = U ΣVT is the reduced singular value decomposition (SVD) of X, then the pseudoinverse is given by X† = V Σ−1UT. We refer to Trefethen & Bau (1997) for an exposition on computational matrix algebra.
The DMD modes and eigenvalues are defined as the eigenvectors and eigenvalues of A. If the snapshots {xk}Nk=1 are linearly independent, the linear relationship is satisfied exactly. Otherwise, A is the least-square solution to (3.1). Note that the if the snapshots consist of DNS data, typically sampled from millions of grid points, i.e. n N , then we have an overconstrained problem and, in addition, the matrix A could be extremely large. However, one does not need to form A explicitly. Instead, a low-rank approximation of it can be formed; see Algorithm 1.
In order to explain the connections between the DMD algorithm(s) and the Koopman operator, we first recall the definition of an observable from Chapter 2:
an observable g : M → R is a scalar-valued function of the state x, which may be a nonlinear mapping. As before, we consider continuously differentiable, real valued observables, even though generalizations are permissible in some cases.
Again, we let g = (g1, . . . , gm) denote a vector of such observables, and we let
yk ≡ g(xk). We assume a time-step ∆t > 0 is given, and that yk and yk+1 correspond to data captured time ∆t apart for all k. Now, we consider applying the DMD algorithm to the yk:s. Suppose the governing dynamics are nonlinear.
The DMD procedure is now performing a linear fit to the data, which arose from a nonlinear system. Under certain conditions however, the eigenvalues obtained from this calculation are true Koopman eigenvalues. Furthermore, the eigenvectors may be used to find the corresponding Koopman eigenfunctions, and consequently, the computed modes will correspond to Koopman modes.
The following theorem from Tu et al. (2014), as stated by Rowley & Dawson (2017), gives conditions under which this relationship holds.
Theorem 1 (Tu et al. 2014). Suppose φ is eigenfunction of the Koopman operator K∆t with discrete-time eigenvalue µ = exp(λ∆t), and suppose φ ∈
span{g1, . . . , gm}, so that
φ(x) = w∗g(x)
for some constants w = (w1, . . . , wm) ∈ Cm. Suppose further that w lies in the range of the data matrix Y , whose columns are yk, k = 1, . . . , N . Then µ is an eigenvalue of A = Y0Y†, and w is a left eigenvector of A: that is, w∗A = µw∗.
Proof. Because K∆tφ = µφ, we have w∗g(F (x)) = µw∗g(x) for all x ∈ M.
In particular, this is true for x = xk, so because yk = g(xk), this becomes w∗yk0 = µw∗yk. Arranging these scalar equations into a row vector given w∗Y0 = µw∗Y , and thus w∗A = µw∗Y Y† = µw∗, where the last equality
holds because w is in the range of Y .
Rowley & Dawson (2017) concludes the implications of the above theorem as follows: to compute Koopman eigenfunctions accurately, we need a sufficiently rich set of observables so that the Koopman eigenfunctions lies in the span of these observables; moreover, we need a sufficiently rich set of data, so that w is in the range of Y .
3.2.2. Extended Dynamic Mode Decomposition
Since the eigenfunctions of the Koopman operator are universal, i.e., they are not related to any observable, there may be betters ways of computing them (via DMD) than using measurements of the state itself as input data.
This is, however, the most commonly encountered approach. In the analysis of fluid flows, one typically uses measurements of the velocity field as input data, i.e. the state itself. The corresponding modes will then essentially correspond to fluctuating parts of the velocity field. The EDMD method, proposed by Williams et al. (2015), provides a means of using DMD with additional observables rather than the input data only. As before, we consider a dataset of snapshot pairs {(xk, x0k)}Nk=1, x0k= Fτ(xk), xk, x0k ∈ Rn. EDMD is a Galerkin weighted residual approach which uses a so-called dictionary of observables to compute Koopman eigenfunctions and eigenvalues. The Koopman modes are then, as usual, obtained by projecting the input data onto these
(hopefully more accurately computed) eigenfunctions. Recall that we have assumed that all observables belong to the vector space F . Let FD ⊂ F be a subspace spanned by a dictionary D ≡ {ψ1, . . . , ψD}, where ψj : M → R.
Then any θ ∈ FD can be expressed as θ(x) = ψ(x)Ta for some a ∈ RD, where ψ(x) = (ψ1(x), . . . , ψD(x))T. Under the action of the Koopman operator, we have
Kτθ(x) = (ψ ◦ F (x))Ta = ψ(x)TKa + r(x),
where K is a finite-dimensional approximation of Kτ, and r(x) is the residual as FD may not be invariant under action of Kτ. One can then formulate a least-squares problem of minimizing the square of the residual according to
N
X
k=1
|r(xk)|2=
N
X
k=1
|((ψ(x0k)T − ψ(xk)TK)a|
to obtain
K = Ψ†xΨx0,
where, again,† denotes the Moore-Penrose pseudoinverse and
Ψx=
ψ(x1)T
... ψ(xN)T
, Ψx0 =
ψ(x01)T
... ψ(x0N)T
.
Here Ψx, Ψx0 ∈ RN ×D, and thus K ∈ RD×D. Typically, we would want to use the state itself as input data when analysing fluid systems, e.g. a discretized velocity field obtained via a DNS. Additionally, we would want more basis functions than the dimension of the state space itself to actually gain from this method as compared to the usual DMD. Moreover, in this case, we typically have D N , so the storage required for the Koopman operator approximation K grows quadratically with the number of basis functions; in the field of machine learning, this is known as the curse of dimensionality.
To circumvent this issue, Williams et al. (2016) proposes the following approach, known as Kernel DMD (KDMD): since the range of K is contained within that of Ψx, any eigenvector of K with a nonzero corresponding eigenvalue µ may be written as v = Z ˆv for some ˆv ∈ CN, where Z is obtained from the reduced SVD of Ψx:
Ψx= QΣZT. The eigenproblem for K can thus be expressed as
µv = Kv
⇔ µZ ˆv = ZΣ†QTΨx0Z ˆv
= Z(Σ†QT)(Ψx0Z)ˆv
= Z(Σ†QT)(Ψx0ΨTx)(QΣ†) ˆv.
An eigenvector of K can thus be computed by first forming K = (Σˆ †QT) ˆA(QΣ†),
where ˆA ≡ Ψx0ΨTx, compute an eigenvector ˆv of ˆK, and set v = Z ˆv. Note that ˆK ∈ RN ×N, so the computational cost is determined by the number of snapshots N rather than the number of observables D in our dictionary D. In addition to ˆA, we define ˆG = ΨxΨTx, and we have
Gˆij = ψ(xj)Tψ(xi); Gˆij= ψ(xj)Tψ(x0i).
Note that ˆG = QΣ2QT, so we can obtain Q and Σ through eigen-decomposition of ˆG. Thus, we can form ˆK in O(N2D) time, which is a large improvement over ‘standard EDMD’, but still unpractical as D can be extremely large. The cure for this is called the kernel trick : the entries of ˆG and ˆA can be computed directly without forming ψ explicitly. Let k(z, w) : M × M → R be a kernel function, then
Gˆij = k(xi, xj), Aˆij= k(x0i, xj). (3.2) The kernel function k implicitly defines FD, the subspace spanned by the elements of ψ(x), and evaluates the inner products in FD in O(N ), rather than O(D), time. An illuminating example is the polynominal kernel, k(z, w) = (1 + zTw)d, where d ∈ N. In the case d = 2, and with z, w ∈ R2, we have
k(z, w) = (1 + z1w1+ z2w2)2
= 1 + 2z1w1+ 2z2w2+ 2z1z2w1w2+ z12w12+ z22w22
= ψ(z)Tψ(w), if ψ(z) =1,√
2z1,√ 2z2,√
2z1z2, z12, z22T
. In general, a polynomial kernel of the form
k(z, w) = 1 + zTwd
corresponds to a dictionary that can represent all polynomials up to and including terms of degree d, and takes only O(n) time (the dimension of the input data) to evaluate. The polynomial kernel was used in Williams et al.
(2016) to illustrate the utility of the kernel approach. However, there are many other choices of kernel functions to consider; we refer the reader to Bishop (2006) (chapter 6) for an in-depth review.
The kernel trick has its roots in machine learning; kernel principal component analysis (PCA) to be precise. The most commonly encountered kernel function in such applications is the Gaussian kernel:
k(z, w) = ψ(z)Tψ(w) = exp
−kz − wk2 2σ02
, (3.3)
where σ0 > 0 is refered to as the bandwidth of the kernel. This is a typical instance of what is known as a radial basis function (RBF), since the value of k only depends on the (squared) distance between x and z. Furthermore, the Gaussian kernel has the appealing property that the corresponding dictionary of observables, or feature map, ψ, is infinite-dimensional. For example, consider
Algorithm 2 Kernel DMD (Williams et al. 2016)
Input: Dataset with snapshot pairs {(xk, x0k)}Nk=1, kernel function k(z, w) Output: Koopman tuple {(µj, φj, vj)}
1: Compute ˆG and ˆA via (3.2), using the kernel k.
2: Let ˆG = QΣ2QT be the eigendecomposition of ˆG, where Σ2 is a diagonal matrix with eigenvalues of ˆG, with diagonal entries ordered according to descending modulus.
3: Let r be the rank of ˆG. Let Qr be the matrix containing the r first columns of Q, and let Σr be the square matrix containing the r first diagonal elements of Σ =p(Σ2).
4: Let ˆK = (Σ−1r QTr) ˆA(QrΣ−1r ) be the Koopman operator approximation.
5: Let ˆK = ˆV ˆD ˆVT be the eigendecomposition of ˆK, where ˆD is a diagonal matrix with diagonal entries µj ordered according to descending modulus;
these are our Koopman eigenvalues.
6: The matrix Φ with Koopman eigenfunctions as columns is now computed through Φ = QrΣrV .ˆ
7: Finally, the Koopman modes are computed by projecting the input data onto the eigenfunctions. That is, [v1, . . . , vN] = [x1, . . . , xN]Φ†.
the case x ∈ R1. Then, by Taylor expansion, it can be shown that
ψ(x) = exp
−y2 2σ0
"
1, s
1 1!σ20x,
s 1 2!σ40x2,
s 1
3!σ60x3, . . .
#
. (3.4)
Other commonly used RBF:s with this property are available, such as multi- quadric kernels: f (r) =p1 + (εr)2 or inverse multiquadric kernels: f (r) =
1 + (εr)2−1
, where r ≡ kz − wk and ε is a positive scaling parameter. The implicit choice of dictionary, through the choice of kernel function, is obviously extremely important. Notwithstanding, we choose to consider with the Gaussian RBF kernel in what follows, for simplicity.
The steps for the Koopman tuple computation based on KDMD are sum- marized in Algorithm 2; for further details on the derivation of the algorithm, we refer to Williams et al. (2016). Being weighted residual Galerkin approaches, EDMD and KDMD converges to the Koopman operator as N → ∞, as shown by Korda & Mezi´c (2017). The total computational cost for the KDMD approach is O(N2max(n, N )), which is the same as for the Exact DMD method described in Algorithm 1.
In what follows, we shall only consider the case xk+1= x0k, and we assume that our snapshots are columns of the matrix X ∈ Rn×N +1. Then, when using any RBF as kernel function, matrices ˆA and ˆG can be assembled efficiently (as opposed to using a for-loop): first, we construct (in partial Matlab notation) Gˆfull≡ [sum(X2, 1)]T1T + 1T[sum(X2, 1)] − 2XTX, (3.5)
