Reduced-order modelling based on Koopman operator theory
Diana A. Bistrian Gabriel Dimitriu Ionel M. Navon


Reduced-order modelling based on Koopman operator theory]Reduced-order modelling based on Koopman operator theory


 (D.A. Bistrian) University Politehnica of TimisoaraDepartment of Electrical Engineering and Industrial InformaticsRevolutiei Nr.5, 331128 Hunedoara, Romania
Corresponding author: diana.bistrian@upt.ro


 (G. Dimitriu) University of Medicine and Pharmacy “Grigore T. Popa”Department of Mathematics and InformaticsUniversitatii 16, 700115 Iasi, Romania
gabriel.dimitriu@umfiasi.ro


 (I.M. Navon) Florida State UniversityDepartment of Scientific ComputingDirac Science Library Building, Tallahassee, FL 32306-4120, USA
Inavon@fsu.edu


1520231-217-27

[2010] 93A30, 70K75,  65C20.


§ ABSTRACT 

The present study focuses on a subject of significant interest in fluid dynamics: the identification of a model with decreased computational complexity from numerical code output using Koopman operator theory. A reduced-order modelling method that incorporates a novel strategy for identifying the most impactful Koopman modes was used to numerically approximate the Koopman composition operator. 




[
    Diana A. Bistrian, Gabriel Dimitriu, Ionel M. Navon
    5 September 2024
=======================================================






§ INTRODUCTION


Despite the fact that complex nonlinear dynamical systems can appear challenging to understand, the existence of similar flow characteristics indicates that a number of different dynamic phenomena are likely governed by the same fundamental processes. Using the modal decomposition method <cit.> is an effective strategy to find a helpful low-dimensional reference frame for capturing prominent dynamical processes. Model order reduction approaches based on modal decomposition have made significant improvements in the recent decade <cit.>.

The choice of an appropriate reduced order basis to describe the system dynamics in relation to the system order reduction is the main topic of this study. 
The Koopman operator theory <cit.> provides the mathematical foundation for determining the reduced-order model of a complex nonlinear dynamical system. 

There are various advantages to transposing the nonlinear dynamics into a reduced-order model, which is a linear model by construction. These include the ability to identify the dominant frequencies, the production of a mathematical reduced-order model with higher fidelity, and a notable increase in computing speed.

The authors have made a substantial contribution to the advancement of modal decomposition techniques <cit.> and the introduction of new numerical algorithms for the modeling of nonlinear dynamical systems with lower computational complexity <cit.>  in the past years.

The present work contains a presentation of the mathematical aspects of modal decomposition technique based on Koopman operator theory <cit.>, with application to the Saint-Venant nonlinear dynamical system model.

The remainder of the article is organized as follows.
Section 2 discusses the mathematical considerations on the Koopman operator theory.
Section 3 presents the numerical method developed for reduced-order modelling.
Section 4 presents the test problem, consisting in nonlinear Saint-Venant equations dynamical model.  A qualitative study of the reduced-order model is conducted in the case of two experiments.
A summary and conclusions are given in Section 5.



§ MATHEMATICAL CONSIDERATIONS ON THE KOOPMAN OPERATOR

 
Let Ω⊂ℝ^n be a compact and non-empty space. 
Let


    L^2( Ω) = {ψ :Ω→ℝ|  ∫_Ω| ψ|.^2dΩ  < ∞}


be a Hilbert space of square integrable functions on  Ω, endowed with the inner product ⟨ψ _i,ψ _j⟩  = ∫_Ωψ _iψ _jdΩ and the norm ψ = √(⟨ψ ,ψ⟩) for ψ∈L^2( Ω).

Let us consider a nonautonomous continuous-time dynamical system on  domain Ω⊂ℝ^n governed by the nonlinear ordinary differential equation


    {[ dy/dt( x,t) = f( y,u,t),   t ∈ℝ_≥ 0;                 y( x,t_0) = y_0( x) ].


where the map f is locally Lipschitz continuous, x∈ℝ^n is the Cartesian coordinate vector, u ∈ℝ^m is the input vector, with n ≫ m. Forward invariance of the set Ω⊂ℝ^n w.r.t. system dynamics (<ref>) is assumed, i.e. any solution y( x,t) ∈Ω , t ≥ 0 holds for all y_0.


Dimensionality reduction in reduced-order modelling.
The principle of modal reduction aims to finding an approximation solution of the form


    {[ y( x,t) ≈∑_j = 1^p a_j( t )ψ _j( x) ,   t ∈ℝ_≥ 0;         da_j( t )/dt = g( t ), a_j( t_0) = a_j^0 ].


expecting that this approximation becomes exact as p →∞, assuring preservation of dynamic stability, computational stability, and a small global approximation error compared to the true solution of (<ref>).



Let us consider a scalar observable function φ :Ω→ℂ, u = φ( y ), y ∈Ω, t ∈ℝ_≥ 0  with a smooth and Lipschitz continuous flow F^t:Ω→Ω:


    F^t( y_0) = y_0 + ∫_t_0^t_0 + tf( y( τ))dτ,


which is forward-complete, i.e. the flow F^t( y ) has a unique solution on ℝ_≥ 0 from any initial condition y_0.

The system class that fits the aforementioned assumption is quite vast and encompasses a wide range of physical systems, including whirling flows, shallow water flows, convection-diffusion processes, and so on.

The Koopman operator describes the propagation of state space observables over time. An observable might be any sort of system measurement or the dynamical reaction of the system. The recurrence of a fixed time-t flow map, i.e. sequential compositions of the map with itself, is assumed to describe the dynamical development. 


Koopman operator.
For dynamical systems of type (<ref>), the semigroup of Koopman operators {𝒦^t}_t ∈ℝ_≥ 0:Ω→Ω acts on scalar observable functions φ :Ω→ℂ by composition with the flow semigroup {F^t}_t ∈ℝ_≥ 0 of the vector field f:


    𝒦^tφ  = φ( F^t).

The Koopman operator is also known as the composition operator.




Linearity of the Koopman operator.
Consider the Koopman operator 𝒦^t and two observables φ _1,φ _2∈Ω and the scalar α∈ℝ. Using (<ref>) it follows that:

    𝒦^t( αφ _1 + βφ _2) = ( αφ _1 + βφ _2)( F^t) = αφ _1( F^t) + βφ _2( F^t) = α𝒦^tφ _1 + β𝒦^tφ _2.




Infinitesimal generator.
Let us assume that there is a generator 𝒢_𝒦:ℱ→Ω, ℱ being the domain of the generator and Ω the Banach space of observables. The operator 𝒢_𝒦 stands as the infinitesimal generator of the time-t indexed semigroup of Koopman operators {𝒦^t}_t ∈ℝ_≥ 0, i.e.

    𝒢_𝒦φ  = lim_t ↘ 0𝒦^tφ  - φ/t = dφ/dt.




Koopman eigenfunction.
An observable ϕ∈Ω is called a Koopman eigenfunction if it satisfies the relation:

    𝒢_𝒦ϕ( y ) = dϕ/dt( y ) = sϕ( y ),


associated with the complex eigenvalue s ∈ℂ.



Koopman mode.
Let ϕ _i∈Ω be an eigenfunction for the Koopman operator, corresponding to eigenvalue λ _i. For an observable φ :Ω→ℂ, the Koopman mode corresponding to ϕ _i is the projection of φ onto span{ϕ _i}. 



Koopman Spectral Decomposition.
Any observable φ :Ω→ℂ  admits a Koopman spectral decomposition of the following form:

    φ( y ) = ∑_j = 1^∞a_j( φ)λ _j^tϕ _j,


where λ _j^t = e^s_jt w.r.t. s_j = σ _j + iω _j with eigen-decay/growth σ _j and eigenfrequencies ω _j.



The Koopman mode decomposition of form (<ref>), first provided in <cit.>, was later coupled with numerical approaches for modal decomposition, such as dynamic mode decomposition <cit.>.

Since for any semigroup of Koopman operators {𝒦^t}_t ∈ℝ_≥ 0, exists an infinitesimal generator 𝒢_𝒦, the following relation is satisfied for any λ ^t = e^st:


    𝒦^tϕ( y ) = ϕ( F^t( y )) = λ ^tϕ( y ).



Let us consider that the space Ω is chosen to be a Banach algebra, i.e. the set of eigenfunctions forms an Abelian semigroup under product of functions. If ϕ _1,ϕ _2∈Ω are two eigenfunctions of the composition operator 𝒦^t with eigenvalues λ _1,λ _2, then the function product ϕ _1ϕ _2 is also an eigenfunction of 𝒦^t with the eigenvalue λ _1λ _2. Thus, products of eigenfunctions are, again, eigenfunctions. It follows that, for any  observable function written in the following form:

    φ( y ) = ∑_j = 0^∞a_j( φ)ϕ _j( y ) ,


the Koopman operator acts as follows:

    𝒦^tφ  = ∑_j = 1^∞a_j( φ)( 𝒦^tϕ _j)  = ∑_j = 1^∞a_j( φ)λ _j^tϕ _j.






§ REDUCED-ORDER MODELLING BASED ON KOOPMAN OPERATOR


Dynamic Mode Decomposition (DMD) <cit.> , is a data-driven approach for estimating the modes and eigenvalues of the Koopman operator without numerically executing a Laplace transform. DMD has emerged as a popular approach for finding spatial-temporal coherent patterns in high-dimensional data, with a strong connection to nonlinear dynamical systems via the Koopman mode theory <cit.>. We present in the following an improved numerical algorithm based on dynamic mode decomposition.

Let us consider a set of observables in the following form:

    u_i( x,t) = u( x,t_i),  t_i = iΔ t,   i = 0,...,N_t


at a constant sampling time Δ t, x representing the spatial coordinates, whether Cartesian or Cylindrical. 

A data matrix whose columns represent the individual data samples, called the snapshot matrix, is constructed in the following manner:

    V = [ [   u_0   u_1   ... u_N_t ]] ∈ℝ^N_x× (N_t + 1)


Each column u_i is a vector with N_x components, representing the numerical measurements.

The Koopman decomposition theory assumes that an infinitesimal operator 𝒦^t exists that maps every vector column onto the next one:

    {u_0, u_1 = 𝒦^tu_0, u_2 = 𝒦^tu_1 = ( 𝒦^t)^2u_0,. .., u_N_t = 𝒦^tu_N_t - 1 = ( 𝒦^t)^N_tu_0}.



Our aim is to build the best numerical approximation of the Koopman operator using the DMD technique. The next step consists in forming two data matrices from the observables sequence, in the form:

    V_0 = [ [       u_0       u_1       ... u_N_t - 1 ]] ∈ℝ^N_x×N_t, V_1 = [ [   u_1   u_2   ... u_N_t ]] ∈ℝ^N_x×N_t.



Assume that over a sufficiently long sequence of snapshots, the latest snapshot may be expressed as a linear combination of preceding vectors, so that:

    u_N_t = c_0u_0 + c_1u_1 + ... + c_N_t - 1u_N_t - 1 + ℛ,


where c_i∈ℝ,i = 0,...,N - 1 and ℛ is the residual vector. The following relations are true:

    {u_1,u_2,...u_N_t} = 𝒦^t{u_0,u_1,...u_N_t - 1} = {u_1,u_2,...,V_0c} + ℛ,


where c = ( [     c_0     c_1     ... c_N - 1 ])^T is the unknown column vector. Eq.(<ref>) is equivalent to the following relation:

    𝒦^tV_0 = V_0𝒮 + ℛ,  𝒮 = ( [         0       ...         0       c_0;         1                   0       c_1;         ⋮         ⋮         ⋮         ⋮;         0         …         1 c_N_t - 1 ]),


where 𝒮 is the companion matrix.

The relationship (<ref>) is true when the residual is minimized. It follows that the vector c must be chosen such that ℛ is orthogonal to span{u_0,...,u_N_t - 1}.
The goal of dynamic mode decomposition is to solve the eigenvalue problem of the companion matrix:

    V_1 = 𝒦^tV_0 = V_0𝒮 + ℛ,


where 𝒮 approximates the eigenvalues of the Koopman operator 𝒦^t when ℛ_2→ 0.

As a direct result of resolving the minimization problem (<ref>), minimizing the residual enhances overall convergence, and so the eigenvalues and eigenvectors of 𝒮 will converge toward the eigenvalues and eigenvectors of the Koopman operator, respectively. 

The advantage of this method is that the Koopman operator has an infinite number of eigenvalues, whereas its DMD approximation is linear and has a finite number of terms. 

Model reduction is highly dependent on the selection of dynamic modes. The superposition of all Koopman modes, weighted by their amplitudes and complex frequencies, approximates the whole data sequence, but some modes contribute insignificantly. In this research, we create a reduced-order model of the data that only includes the most important modes that make a substantial contribution to the representation of the data, which we refer to as the leading modes.

The data snapshots at every time step will be represented as a Koopman spectral decomposition of the form:

    u_DMD( x,t_i) = ∑_j = 1^N_DMDa_j( t_i)λ _j^i - 1ϕ _j( x) ,   1mu i ∈{1,...,N_t},   1mut_i∈{t_1,...,t_N_t},


where N_DMD≪N_t represents the number of Koopman leading modes ϕ( x) involved in the spectral decomposition of data snapshots, λ _j are the Koopman eigenvalues, and a_j∈ℂ are the modal amplitudes of the Koopman modes, respectively.


 The leading modes indicate a subset of Koopman modes that will be chosen from all computed DMD modes using an original criterion, discussed in the following.

We define the weight of each Koopman mode as follows:

    w𝒦_j = ∫_Δ t^t_N_t∑_i = 1^N_ta_j( t )λ _j^i - 1 dt,


where λ _j are the Koopman eigenvalues, and a_j∈ℂ are the modal amplitudes of the Koopman modes, respectively.


Let

    Er_DMD = u( x) - u_DMD( x)_2/u( x)_2,


be the relative error of the difference between the variables of the full model and approximate DMD solutions over the exact one, where u( x) represents the full solution of the model and u_DMD( x) represents the reduced order solution. 

The leading dynamic modes and their related frequencies are chosen in descending order of the modal entropy, until a minimal relative error of the reduced-order model is obtained. To produce the reduced-order model amounts to finding the solution to the following minimisation problem:


    {[ Find N_DMD∈ N, w.r.t.  u_DMD( x,t_i) = ∑_j = 1^N_DMDa_jϕ _j( x)λ _j^i - 1 ,;                                       i ∈{1,...,N_t},  t_i∈{t_1,...,t_N_t},;              Subject to  min_N_DMD{w𝒦_1 > w𝒦_2 > ... > w𝒦_N_DMD, Er_DMD≤ε}. ].



As a consequence, the modes and frequencies with the highest effect on approximation accuracy are selected to be included in the model with a reduced computational complexity.



§ REDUCED-ORDER MODELLING OF SAINT-VENANT EQUATIONS MODEL


The test problem used in this paper consists of the nonlinear Saint-Venant equations (also called the shallow water equations <cit.>) in a channel on the rotating earth:



    ∂( ũh)/∂ t + ∂( ũ^2h + gh^2/2)/∂ x + ∂( ũṽh)/∂ y = h( fṽ - g∂ H/∂ x),




    ∂( ṽh)/∂ t + ∂( ũṽh)/∂ x + ∂( ṽ^2h + gh^2/2)/∂ y = h(  - fũ - g∂ H/∂ y),




    ∂h/∂ t + ∂( ũh)/∂ x + ∂( ṽh)/∂ y = 0,


where ũ and ṽ are the velocity components in the x̃ and ỹ axis directions respectively, h̃ represents the
depth of the fluid, H( x,y) is the the orography field, f̃ is the Coriolis factor and g is the acceleration of gravity.

The reference computational configuration is the rectangular 2D domain Ω = [ 0,L_max] ×[0,D_max]. Subscripts represent the derivatives with respect to time and the streamwise and spanwise coordinates.

The Coriolis parameter is modelled as varying linearly in the spanwise direction, such that


    f = f_0 + β (y - D_max),


where f_0,β are constants, L_max,D_max are the dimensions of the rectangular domain of integration.

The height of the orography is given by the fixed two-dimensional field


    H( x,y) = αe^y^2 - x^2.



The model (<ref>)-(<ref>) is associated with periodic boundary conditions in the  x̃-direction and solid wall boundary condition in
the ỹ-direction:


    ũ( 0,ỹ,t̃) = ũ( L_max,ỹ,t̃), ṽ( x̃,0,t̃) = ṽ( x̃,D_max,t̃) = 0,


and also with the initial Grammeltvedt type condition <cit.> as the initial height field, which propagates the energy in wave number one, in the
streamwise direction:


    h_0( x̃,ỹ) = H_0 + H_1tanh( 10(D_max/2 - ỹ)/D_max) + H_2sin( 2πx̃/L_max)cosh ^ - 2( 20(D_max/2 - ỹ)/D_max).



Using the geostrophic relationship ũ =  - h̃_ỹ( g/f̃), ṽ = h̃_x̃(
g/f̃), the initial velocity fields are derived as:


    u_0( x̃,ỹ) =  - g/f̃10H_1/D_max( tanh^2( 5D_max - 10ỹ/D_max) - 1) -



    18g/f̃H_2sinh( 10D_max - 20ỹ/D_max)sin( 2πx̃/L_max)/D_maxcosh^3( 10D_max - 20ỹ/D_max),



    v_0( x̃,ỹ) = 2πH_2g/f̃L_maxcos 1mu( 2πx̃/L_max)cosh ^ - 2( 20(D_max/2 - ỹ)/D_max).



The constants used for the test problem are


    f_0 = 10^ - 4s^ - 1,   α  = 4000,   β  = 1.5 ×10^ - 11s^ - 1m^ - 1,   g = 9.81ms^ - 1,



    D_max = 60×10^3m,  L_max = 265×10^3m,



    H_0 = 10 ×10^3m,  H_1 = -700m,  H_2 = -400m.




The error of the numerical algorithm is set to be less than ε=10^-7. A non-dimensional analysis was performed to assess the performances of the reduced-order shallow water model. 
Reference quantities of the dependent and independent variables in the shallow water model are considered, i.e. the length scale
L_ref = L_max and the reference units for the height and velocities, respectively, are given by the initial conditions h_ref =
h_0, u_ref = u_0. A typical time scale is also considered, assuming the form t_ref = L_ref/u_ref. 

In order to make the
system of equations (<ref>)-(<ref>) non-dimensional, the non-dimensional variables

    ( t,x,y) = ( t̃/t_ref,x̃/L_ref,ỹ/L_ref),  ( h,u,v) = ( h̃/h_ref,ũ/u_ref,ṽ/u_ref)

are introduced. 

The numerical results are obtained employing a Lax-Wendroff finite difference discretization scheme <cit.> and used in further numerical experiments in dimensionless form. The training data
comprises a number of 289 unsteady solutions of the two-dimensional shallow water equations model (<ref>)-(<ref>), at regularly spaced time
intervals of Δ t = 1800s for each solution variable.

The numerical results of two tests illustrating the computing performance of the approach are presented below. In the first experiment, the threshold is set to be ε  = 10^ - 3 for solving the optimization problem (<ref>). In the second experiment, the threshold is set at ε  = 10^ - 4 for solving the optimization problem (<ref>).

Figures <ref>–<ref> present the spectrum of Koopman
decomposition eigenvalues, of geopotential height field h, streamwise field u and spanwise field v, respectively, in the case of two experiments, and the leading Koopman modes selected by resolving the optimization problem (<ref>). In the second experiment, extra modes are selected (darker colored dots) to improve the reduced-order model precision.











The representation of the height field compared to its reduced-order model is displayed in Figures <ref>–<ref>, in the case of both experiments.







The vorticity field compared to its reduced-order model is illustrated in Figures <ref>–<ref>, in the case of both experiments, at different time instances.






Table <ref> presents the percentage reduction of the computational complexity of the reduced-order model, in the two experiments performed.












§ CONCLUSIONS


The current study concentrated on a topic of significant interest in fluid dynamics: the identification of a model of reduced computational complexity from numerical code output, based on Koopman operator Theory.
The full model consisted in the Saint-Venent equations model, that have been computed using a Lax-Wendroff finite difference discretization scheme. The Koopman composition operator have been numerically approximated with the algorithm of reduced-order modelling, endowed with a novel criterion of selection of the most influential Koopman modes, based on the modes weights. It automatically selects the most representative Koopman modes, even if they exhibit rapid development with lower amplitudes or are composed of high amplitude fast damped modes.

Two tests were carried out in order to evaluate the algorithm's computing efficiency in order to enhance the reduced-order model precision. It was demonstrated that the model rank may be decreased by up to 92% without compromising model accuracy.

This approach is a useful tool for creating reduced-order models of complex flow fields characterized by non-linear models.



99



Holmes1996 Holmes, P., Lumley, J., Berkooz, G., Turbulence, coherent structures, dynamical systems and symmetry, Cambridge University Press, 1996.

Kaiser2018 Kaiser, E., Kutz, J., Brunton, S., Sparse identification of nonlinear dynamics for model predictive control in the low-data limit, Proceedings Of The Royal Society A, 474, 2018.
    
ChKutz2019Champion, K., Brunton, S., Kutz, J., Discovery of nonlinear multiscale systems: Sampling strategies and embeddings, SIAM Journal On Applied Dynamical Systems, 18, 312-333, 2019.

Arcucci2023 Arcucci, R., Xiao, D., Fang, F., Navon, I.M., Wu, P., Pain, C.C., Guo, Y.K., A reduced order with data assimilation model: Theory and practice, Computers & Fluids 257, 105862, 2023.

Mou2023 Mou, C., Merzari, E., San, O., Iliescu, T., An energy-based lengthscale for reduced order models of turbulent flows, Nuclear Engineering and Design 412, 112454, 2023.
    
Iliescu2022Iliescu, T., ROM Closures and Stabilizations for Under-Resolved Turbulent Flows, 2022 Spring Central Sectional Meeting, 2022.
    
Daba2023 Dabaghian, P.H., Ahmed, S.E., San, O.,  Nonintrusive Reduced Order Modelling of Convective Boussinesq Flows, International Journal of Computational Fluid Dynamics 36 (7), 578-598, 2022.
     
Sanfilippo2023 Sanfilippo, A., Moore, I., Ballarin, F., Iliescu, T., Approximate deconvolution Leray reduced order model for convection-dominated flows, Finite Elements in Analysis and Design 226, 104021, 2023.

Koopm1931 Koopman, B., Hamiltonian systems and transformations in Hilbert space, Proc. Nat. Acad. Sci., 17, 315-318, 1931.

BistrianDimitriuNavon2019 Bistrian, D.A., Dimitriu, G., Navon, I.M., Processing epidemiological data using Dynamic Mode Decomposition method, AIP Conference Proceedings 2164, 080002, 2019.

BistrianDimitriuNavon2020a Bistrian, D.A., Dimitriu, G., Navon, I.M. Modeling dynamic patterns from COVID-19 data Using Randomized Dynamic Mode Decomposition in predictive Mode and ARIMA, AIP Conference Proceedings 2302, 080002, 2020.

BistrianDimitriuNavon2020b Bistrian, D.A., Dimitriu, G., Navon, I.M., Application of deterministic and randomized dynamic mode decomposition in epidemiology and fluid dynamics, An. Ştiinţ. Univ. Al. I. Cuza Iaşi. Mat. (N.S.), Tomul LXVI, F. 2, 2020.

BistrianNavon2017b Bistrian, D.A., Navon, I.M., Randomized Dynamic Mode Decomposition for non-intrusive reduced order modelling, International Journal for Numerical Methods in Engineering, ISSN 0029-5981, Volume: 112, Issue: 1, Page:3-25, 2017.
    
Bistrian2022a Bistrian, D.A., High-Fidelity Digital Twin Data Models by Randomized Dynamic Mode Decomposition and Deep Learning with Applications in Fluid Dynamics, Modelling, Volume 3, Issue 3, pp. 314-332, 2022.

Bistrian2022b Bistrian, D.A., Mathematical considerations on Randomized Orthogonal Decomposition method for developing twin data models, Transylvanian Journal of Mathematics and Mechanics, Volume 14, Number 2, pp. 105-115, 2022.

Mezic2005 Mezic, I., Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics, Volume 41, pp. 309-325, 2005.

SchmidSesterhen2008 Schmid, P.J., Sesterhenn, J., Dynamic Mode Decomposition of Numerical and Experimental Data, 61st Annual Meeting of the APS Division of Fluid Dynamics, San Antonio, Texas, Vol. 53(15), American Physical Society, 2008.

Schmid2010 Schmid, P.J., Dynamic Mode Decomposition of Numerical and Experimental Data, Journal of Fluid Mechanics 656: 5–28, 2010.

Mezic2013 Mezic, I., Analysis of fluid flows via spectral properties of the Koopman operator, Annual Review of Fluid Mechanics, Volume 45, pp. 357-378, 2013.

Tu2014 Tu, J., Rowley, C., Luchtenburg, D., Brunton, S., Kutz, J., On dynamic mode decomposition: Theory and applications, Journal Of Computational Dynamics, 1, 391-421, 2014.
    
Bistrian2015Bistrian, D., Navon, I., An improved algorithm for the shallow water equations model reduction: Dynamic Mode Decomposition vs POD, International Journal For Numerical Methods In Fluids, 78, 552-580, 2015.

chenoptimal2012Chen, K., Tu, J., Rowley, C., Variants of dynamic mode decomposition: boundary condition, Koopman and Fourier analyses, Nonlinear Science, 22, 887-915, 2012.

SaintVenant Saint-Venant, A.J.C., Barré de, J.C.,  Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et a l'introduction de marées dans leurs lits, Comptes rendus de l'Académie des Sciences 73, 237-240, 1871.

Grammeltvedt1969 Grammeltvedt, A., A survey of finite-difference schemes
for the primitive equations for a barotropic Fluid, Monthly Weather Review 97 (5): 384–404, 1969.

Brass2011Brass, H., Petras, K., Quadrature Theory: The Theory of Numerical Integration on a Compact Interval, American Mathematical Soc., 2011.


