The eigenspectrum of the Koopman operator enables the decomposition of nonlinear dynamics into a sum of nonlinear functions of the state space with purely exponential and sinusoidal time dependence. For a limited number of dynamical systems, it is possible to find these Koopman eigenfunctions exactly and analytically. Here, this is done for the Korteweg–de Vries equation on a periodic interval using the periodic inverse scattering transform and some concepts of algebraic geometry. To the authors’ knowledge, this is the first complete Koopman analysis of a partial differential equation, which does not have a trivial global attractor. The results are shown to match the frequencies computed by the data-driven method of dynamic mode decomposition (DMD). We demonstrate that in general, DMD gives a large number of eigenvalues near the imaginary axis and show how these should be interpreted in this setting.

Topics
Non linear dynamics, Algorithms and data structure, Computational methods, Operator theory, Partial differential equations, Riemann surfaces, Algebraic geometry, Inverse scattering, Korteweg-de Vries equation

The Koopman operator was introduced by Koopman1 to describe the nonlinear behavior of a dynamical system as the linear evolution of nonlinear observables of that system. It is well known that a finite dimensional nonlinear system can be converted into an infinite dimensional linear system; Koopman analysis extends this to infinite dimensional nonlinear systems. Recent interest in the Koopman operator was initiated by Mezić,2,3 and it has become popular for its close connection to the computational method of dynamic mode decomposition.

DMD was originally invented by Schmid4 as a method for simultaneously extracting the important spatial and temporal features of a timeseries. It has been successfully applied to a wide array of numerical, observational, and experimental data, most notably in fluid dynamics.5–7 Under certain conditions,8 the results of DMD can be seen as a numerical approximation to the Koopman modes and eigenvalues of the underlying dynamical system, and therefore, understanding the Koopman operator aids interpretation of the results of DMD. In particular, to understand the spatial patterns called DMD modes, it would be helpful to have analytic results on nonlinear partial differential equations (PDEs).

Several authors9–11 have successfully performed Koopman analysis on the Burgers equation, a dissipative nonlinear PDE, which can be transformed into a linear heat equation with a suitable change of variables. Nakao and Mezić12 considered the Burgers equation as well as a non-trivial transformation of it to the phase-diffusion equation. Though insightful, these PDEs admit only a single steady global attractor and are strongly dissipative, which precludes much interesting nonlinear behavior relevant to applying DMD to situations involving sustained waves, a common use-case.

Parker and Page13 considered the Korteweg–de Vries (KdV) equation, an integrable partial differential equation of one variable. This Hamiltonian dynamical system behaves very differently from the dissipative systems mentioned above. In that work, some Koopman eigenfunctions were found for soliton solutions of the KdV equation on the real line, which excluded a large class of solutions; in particular, this excluded the spatially periodic solutions, which are the natural choice for computer simulations of solitons. It was argued that the simplest periodic solutions, cnoidal waves, give purely imaginary Koopman eigenvalues, in sharp contrast to the isolated solitons, which have purely real Koopman eigenvalues, despite being a naturally limiting case of the former.

A periodic domain is the natural setting for numerical solutions of one-dimensional PDEs, such as the KdV equation. Indeed, early pioneering numerical simulations of the KdV equation on a periodic domain14 were used to shed light on the famous Fermi–Pasta–Ulam–Tsingou problem,15 in which the initial state of a system recurs to arbitrary precision after complex nonlinear dynamics. It was later proven16 that this is because the solutions lie on quasiperiodic invariant tori. We will show that it is possible to use DMD to determine the underlying frequencies. In fact, when the dynamics are confined to an invariant torus, DMD is equivalent to a Fourier transform in time, though it generalizes when other modes grow or decay.

In the present paper, we study a particular but very general class of solutions of the KdV equation on a periodic interval, for which we are able to define Koopman eigenfunctions. These eigenfunctions require the evaluation of contour integrals on Riemann surfaces, which must be performed numerically in all but the simplest of cases. This method is not, therefore, recommended as a general approach for nonlinear PDEs, where DMD could easily be applied. However, the semi-formal mathematical treatment presented here gives instructive results: we will see that in this case, Koopman eigenvalues necessary for the decomposition of the state of the system densely fill the imaginary axis, and therefore, the results of DMD are subtle to interpret. This paper proceeds as follows: in Sec. II, the KdV equation is introduced, and its history and significance are briefly described; in Sec. III, the Koopman operator and its spectrum are defined; in Sec. IV, we then derive Koopman eigenfunctions for the KdV equation; Sec. V presents the results of applying these to an example solution, which is compared to the numerical results of DMD in Sec. VI. Concluding remarks are given in Sec. VII.

The KdV equation is given by
4 ∂ t u ( x , t ) = 6 u ( x , t ) ∂ x u ( x , t ) + ∂ x 3 u ( x , t ) , u ( x , 0 ) = u 0 ( x ) .
(1)
Note that many different conventions are employed in the literature, and notably, solitons can be either positive or negative depending on this choice. Here, we follow Belokolos et al.17 The equation was derived by Korteweg and de Vries18 to describe the weakly nonlinear evolution of shallow water waves propagating in one direction. The pioneering computational results of Zabusky and Kruskal14 demonstrated the existence of soliton solutions when this equation is solved on a finite periodic domain. For an infinite domain, the celebrated inverse scattering method of Gardner et al.19 gives a straightforward procedure to solve the equation and provides intuitive interpretations for the existence of solitons as conserved quantities. Analytical results on a periodic interval have proven much more complicated, despite the early computation successes and the well known cnoidal wave solution.

In addition to shallow water waves, the KdV equation naturally arises from weakly nonlinear theory in many physically relevant flows.20–22 Our interest, however, derives from the fact that the equation admits complex, nonlinear but non-chaotic solutions amenable to analytic treatment. Unlike other PDEs for which Koopman spectra have been derived, it is a Hamiltonian system with an infinite number of conserved quantities, rather than having one unique attractor.

The solution to (1) is well-posed on T = R / 2 π Z⁠, the periodic domain of length 2 π⁠, for initial conditions in Sobolev spaces H s ( T , R ) with s ≥ − 1⁠,23 with a well-defined evolution operator S t : H s ( T , R ) → H s ( T , R ) for each time t⁠. A single periodic wave of sufficient amplitude breaks down into a spatially periodic solution, with quasiperiodic behavior in time. That is to say, solutions lie on an invariant torus, as is usual for integrable Hamiltonian systems, and these invariant tori foliate phase space. In fact, any solution of (1) can be approximated to arbitrary precision as an invariant M-torus,16, M ∈ N⁠. These are the so-called finite gap solutions, which will be the focus of the present study. Practically, the evolution of an arbitrary initial condition can be approximated by truncating the scattering data at a judiciously chosen M⁠.24 

Let us, therefore, define our solution space Ω M to be the subset of L 2 ( T , R ) for which there are M-gap solutions (this terminology should become clearer in Sec. IV). This is a well-posed invariant subspace. For convergent Koopman decompositions, it will be necessary to further restrict this space, in a manner analogous to that of Balabane et al.,11 in Sec. V.

It will be useful later to make the so-called Hirota transformation u ↦ ϑ defined by25,
u ( x , t ) = 2 ∂ x 2 log ⁡ ϑ ( x , t )
(2)
so that (1) becomes
4 ϑ ∂ x ∂ t ϑ − 4 ∂ x ϑ ∂ t ϑ − 3 ( ∂ x 2 ϑ ) 2 + 4 ∂ x ϑ ∂ x 3 ϑ − ϑ ∂ x 4 ϑ = 0.
(3)
This is analogous to the Cole–Hopf transformation exploited by previous authors for the Burgers equation,9–11 though in this case, the transformation seems at first glance to have made the equation more complicated. The utility comes from the fact that ϑ can be expressed as a Riemann theta function, as explained in Sec. IV.
Let O M be the space of continuous maps Ω M → C⁠, which are called observables of the system. The Koopman operator, a composition operator for dynamical systems, is defined for each t ≥ 0 by
K M t : O M → O M , [ K M t ϕ ] ( u 0 ) = ϕ ( u ( ⋅ , t ) )
(4)
for any observable ϕ ∈ O M⁠, where the evolution of u ( x , t ) is governed by (1). The Koopman operator is a linear operator, amenable to spectral theory. The Koopman eigenvalues λ ∈ C and Koopman eigenfunctions φ ∈ O M satisfy
[ K M t φ ] ( u 0 ) = e λ t φ ( u 0 )
(5)
or, equivalently,
φ ( u ( ⋅ , t ) ) = e λ t φ ( u 0 ) .
(6)
A simple example of a Koopman eigenfunction would be any conserved quantity of the dynamics, with eigenvalue λ = 0⁠. More generally, they are any observable for which the temporal behavior is purely (complex) exponential as the state evolves. Since the system we study is Hamiltonian, we expect only purely imaginary eigenvalues λ = i ω⁠, giving purely oscillatory behavior. In certain circumstances, it may be possible that the Koopman eigenfunctions form a basis for O M⁠, in which case we can decompose all other observables as a sum over Koopman eigenfunctions. In particular, we are interested in whether it is possible to write the state of the system u as a convergent sum,
u ( x , t ) = ∑ ν c ν ( x ) φ ν ( u 0 ) e i ω ν t .
(7)
Here, the c ν : T → C are called Koopman modes, first defined by Mezić,2 which encode spatial information for each eigenvalue and are independent of the particular choice of initial condition u 0⁠, whose contribution is included in the value of φ ν ( u 0 )⁠. If this is possible, it means that the dynamics of (1) can be decomposed as a sum over nonlinear functions whose temporal behavior is purely oscillatory in time.

It would take a whole textbook to fully explain the periodic inverse scattering transform. We refer readers to Novikov et al.,26 Belokolos et al.,17 and Osborne27 for accessible introductions, including the necessary background of Riemann surfaces and theta functions, though note the differing notations and conventions between these (we follow the notation of Belokolos et al.17). Here, we give a summary of the relevant results for the KdV equation, which are implemented in the Mathematica notebook given in the supplementary material.

One of the key results in the solution of the KdV equation was the discovery of a Lax pair:28 a pair of linear operators L ( u ) , A ( u ) : L 2 ( R ) → L 2 ( R ) such that
d L d t = L ( t ) ° A ( t ) − A ( t ) ° L ( t ) ,
(8)
where L ( t ) := L ( u ( ⋅ , t ) ) etc. In the case of the KdV equation (1),
L ( u ) = − d 2 d x 2 − u ,
(9)
A ( u ) = d 3 d x 3 + 3 4 ( u d d x + d d x u ) .
(10)
The operator (10) is skew-adjoint. The operator (9) is the well known self-adjoint Schrödinger operator, with potential u⁠. From (8), it can be shown that the spectrum of L is independent of t when u ( x , t ) satisfies (1). Finding eigenvalues λ ∈ R and eigenfunctions ψ ∈ C ∞ reduces to the Sturm–Liouville problem
ψ ″ + u ψ + λ ψ = 0.
(11)
In the case of periodic potential u ( x + 2 π ) = u ( x )⁠, Eq. (9) is known as Hill’s operator and has been widely studied.29 The admissible eigenvalues λ for a bounded eigenfunction ψ reside in intervals [ E 1 , E 2 ]⁠, [ E 3 , E 4 ]⁠, …, where − ∞ < E 1 < E 2 ≤ E 3 < E 4 ≤ E 5 < ⋯⁠. Outside these regions, only unbounded solutions are possible, and these are termed forbidden gaps. The E k are the values of λ for which F ( λ )⁠, defined as half the trace of the monodromy matrix of (11), is ± 129 (see Fig. 1). Though the monodromy matrix is not invariant under the dynamics (1), its trace, and, therefore, also the E k are invariant. We explicitly consider only the case when there is a finite number g of non-degenerate forbidden gaps ( E 2 , E 3 )⁠, ( E 4 , E 5 )⁠, …, so that E 2 k = E 2 k + 1 for all k > g⁠.
FIG. 1.
F ( λ ), half the trace of the monodromy matrix, against real λ, for u 0 ( x ) = sin ⁡ x. Allowable values of λ for a bounded eigenfunction are when | F ( λ ) | ≤ 1. The forbidden regions ( − ∞ , E 1 ), ( E 2 , E 3 ), …, are called gaps.
View largeDownload slide

F ( λ )⁠, half the trace of the monodromy matrix, against real λ⁠, for u 0 ( x ) = sin ⁡ x⁠. Allowable values of λ for a bounded eigenfunction are when | F ( λ ) | ≤ 1⁠. The forbidden regions ( − ∞ , E 1 )⁠, ( E 2 , E 3 )⁠, …, are called gaps.

FIG. 1.
F ( λ ), half the trace of the monodromy matrix, against real λ, for u 0 ( x ) = sin ⁡ x. Allowable values of λ for a bounded eigenfunction are when | F ( λ ) | ≤ 1. The forbidden regions ( − ∞ , E 1 ), ( E 2 , E 3 ), …, are called gaps.
View largeDownload slide

F ( λ )⁠, half the trace of the monodromy matrix, against real λ⁠, for u 0 ( x ) = sin ⁡ x⁠. Allowable values of λ for a bounded eigenfunction are when | F ( λ ) | ≤ 1⁠. The forbidden regions ( − ∞ , E 1 )⁠, ( E 2 , E 3 )⁠, …, are called gaps.

Close modal
The hyperelliptic curve
μ 2 = ∏ j = 1 2 g + 1 ( λ − E j )
(12)
defines a Riemann surface of two sheets, with a branch point at each E j (see Fig. 2). The genus of this surface is simply the number of gaps g⁠. It is then possible to define a basis of contours a j and b j (⁠ 1 ≤ j ≤ g⁠) for the Riemann surface such that any contour can be expressed, up to continuous deformations, as a sum of a j and b j⁠. Such a choice of basis is not unique and will have implications for the final results—see discussions of the wave basis and the soliton basis in Osborne.27 Our convention is shown in Fig. 3. We also define a basis of holomorphic differentials on this surface,
ω ˚ k = λ g − k d λ μ ,
(13)
and then make a linear transformation to the canonical basis ω k = ∑ j c k j ω ˚ j such that
∫ a j ω k = 2 π i δ j k , j , k = 1 , … , g .
(14)
In this new basis, we define the period matrix of the Riemann surface
B j k = ∫ b j ω k .
(15)
It can be shown that this matrix is symmetric B j k = B k j with all entries having a strictly negative real part. It is then the case that the variable ϑ for the transformed KdV equation can be written as
ϑ ( x , t ) = θ ( U x + W t + D , B ) ,
(16)
where we define the Riemann theta function
θ ( z , B ) = ∑ m ∈ Z g exp ⁡ ( 1 2 m T B m + z T m ) .
(17)
For B real and z imaginary, this gives real values by symmetry. The wavenumber vector U and the frequency vector W are calculated as
U j = 2 i c j 1 ,
(18)
W j = − 2 i ( c j 2 + 1 2 c j 1 ∑ k = 1 2 g + 1 E k ) .
(19)
By construction, U j must be integers, but the frequencies W j will be incommensurate, in general. Both U and W are purely imaginary. Since they are calculated only from the E k⁠, all of B⁠, U⁠, and W are constant as the system evolves. The value of the vector of phases D⁠, conversely, depends on the particular state at time t = 0 (and its evolution is absorbed into W⁠).
FIG. 2.
Riemann surface of genus 2 for the initial condition u 0 ( x ) = sin ⁡ x, truncated to a two-gap solution. The horizontal axes show the real and imaginary parts of λ, and the vertical axis shows the real part of μ, as defined by (12).
View largeDownload slide

Riemann surface of genus 2 for the initial condition u 0 ( x ) = sin ⁡ x⁠, truncated to a two-gap solution. The horizontal axes show the real and imaginary parts of λ⁠, and the vertical axis shows the real part of μ⁠, as defined by (12).

FIG. 2.
Riemann surface of genus 2 for the initial condition u 0 ( x ) = sin ⁡ x, truncated to a two-gap solution. The horizontal axes show the real and imaginary parts of λ, and the vertical axis shows the real part of μ, as defined by (12).
View largeDownload slide

Riemann surface of genus 2 for the initial condition u 0 ( x ) = sin ⁡ x⁠, truncated to a two-gap solution. The horizontal axes show the real and imaginary parts of λ⁠, and the vertical axis shows the real part of μ⁠, as defined by (12).

Close modal
FIG. 3.
Schematic of the basis of contours for a g = 2 Riemann surface, where dashed lines show the parts of the contour taken on the second sheet. The colors are as for Fig. 2.
View largeDownload slide

Schematic of the basis of contours for a g = 2 Riemann surface, where dashed lines show the parts of the contour taken on the second sheet. The colors are as for Fig. 2.

FIG. 3.
Schematic of the basis of contours for a g = 2 Riemann surface, where dashed lines show the parts of the contour taken on the second sheet. The colors are as for Fig. 2.
View largeDownload slide

Schematic of the basis of contours for a g = 2 Riemann surface, where dashed lines show the parts of the contour taken on the second sheet. The colors are as for Fig. 2.

Close modal
To find the phases D⁠, it is necessary to define a second set of eigenvalues for the Sturm–Liouville problem (11), now with the boundary conditions ψ ( 0 ) = ψ ( 2 π ) = 0⁠. This discrete set of eigenvalues λ 1 , λ 2 , … lies in the gaps so that E 2 < λ 1 < E 3⁠, E 4 < λ 2 < E 5⁠, etc. These eigenvalues are not constant as the state evolves and depend on the time of measurement. The formula for D is then given by17 
D j = − ∑ k = 1 g ∫ ∞ λ k ω j − ∑ k = 1 g B j k + i π j ,
(20)
where here, λ k represents a point on the surface with λ = λ k⁠, with care be taken to evaluate the integral on the correct sheet. D is also purely imaginary.
From (16), we can view D as evolving in time as D + W t⁠, and so with the goal of writing a Koopman decomposition for ϑ via (17), we define a Koopman eigenfunction for (1) for each m ∈ Z g⁠,
φ m ( u 0 ) := exp ⁡ ( D T m ) ,
(21)
which then evolves as
φ m ( u ( ⋅ , t ) ) = exp ⁡ ( W T m t ) exp ⁡ ( D T m )
and, thus, has Koopman eigenvalue W T m⁠, an integer linear combination of the fundamental frequencies. The Koopman eigenfunction is written as a function of the initial condition u 0⁠; in this case, the dependence arises through D from using u 0 as the potential in the Sturm–Liouville problem. The eigenvalues, by contrast, depend only on W and, therefore, only on the invariant Riemann surface, rather than the precise choice of an initial condition.

As the frequencies W j are incommensurate, for genus g = 2 and greater, these Koopman eigenvalues densely fill the imaginary axis. This is a significant complication over previously studied PDEs. Note that φ m 1 ( u ) φ m 2 ( u ) = φ m 1 + m 2 ( u )⁠, and φ 0 ( u ) = 1⁠. Indeed, rather generally, integer linear combinations of Koopman eigenvalues are eigenvalues and products of Koopman eigenfunctions are Koopman eigenfunctions,30 though these eigenfunctions are not necessarily useful for decompositions in all cases.

Clearly, the expression (16) is directly a convergent Koopman decomposition for ϑ ( x , t )⁠, as it is a sum of terms whose time dependence is purely exponential. We can write it as
ϑ ( x , 0 ) = ∑ m ∈ Z g e U T m x exp ⁡ ( 1 2 m T B m ) φ m ( u 0 ) .
(22)
It is somewhat more involved to obtain a decomposition of u⁠, but this is still possible so long as ϑ is sufficiently small, via (2),
u 0 ( x ) = 2 ∂ x 2 log ⁡ [ 1 + ∑ m ∈ Z g ∖ { 0 } exp ⁡ ( 1 2 m T B m + U T m x ) φ m ( u 0 ) ] = − 2 ∂ x 2 [ ∑ n = 1 ∞ ( − 1 ) n n ∑ m 1 , … , m n ∈ Z g ∖ { 0 } exp ⁡ { ∑ q = 1 n ( 1 2 m q T B m q + U T m q x ) } ∏ q = 1 n φ m q ( u 0 ) ] = ∑ m ∈ Z g ( − 2 ( U T m ) 2 e U T m x ∑ n = 1 ∞ ( − 1 ) n n ∑ m 1 , … , m n ∈ Z g ∖ { 0 } ∑ q m q = m exp ⁡ { ∑ q = 1 n 1 2 m q T B m q } ) φ m ( u 0 ) .
(23)
This complicated series is absolutely convergent when 0 < ϑ ( x , t ) < 2⁠. For larger ϑ⁠, other expansions could be found, using only the Koopman eigenfunctions given in Sec. IV. To summarize this expression, we have found a Koopman decomposition for u using Koopman eigenfunctions φ m ( u 0 ) with Koopman eigenvalues W T m⁠. The corresponding Koopman modes are
− 2 ( U T m ) 2 e U T m x ∑ n = 1 ∞ ( − 1 ) n n ∑ m 1 , … , m n ∈ Z g ∖ { 0 } ∑ q m q = m exp ⁡ { ∑ q = 1 n 1 2 m q T B m q } .
(24)
Notice that these are purely sinusoidal in x⁠. The Koopman modes depend on U and B⁠, which are functions of the Riemann surface and, therefore, of the g-torus to which the dynamics are constrained in phase space, but the Koopman modes do not depend on the choice of initial conditions beyond this.
As a concrete example, we consider the initial condition u 0 ( x ) = sin ⁡ x⁠. Despite the simplicity of this choice, it gives an apparently infinite number of non-degenerate gaps (see Fig. 1), but using only g = 2 or g = 3 results in good agreement. With g = 4⁠, not shown here, the reconstructed solution is virtually indistinguishable from the initial condition. Only two solitons are visible per spatial period in a numerical solution; the genus of the Riemann surface is not the number of solitons. For g = 2⁠, we find numerically that
U = ( − 1 i , − 2 i ) , W ≈ ( − 0.036 , 1.930 i ) ,
(25)
B ≈ ( − 3.171 − 1.930 − 1.930 − 7.247 ) ,
(26)
and for g = 3⁠,
U = ( − 1 i , − 2 i , − 3 i ) , W ≈ ( − 0.036 i , 1.931 i , 6.741 i ) ,
(27)
B ≈ ( − 3.171 − 1.929 − 1.321 − 1.929 − 7.247 − 3.190 − 1.321 − 3.190 − 12.810 ) .
(28)

The reconstructed u ( x , t ) given by a finite truncation of the Koopman decomposition (23) is shown in Figs. 4 and 5. Including terms higher than n = 4 in the series may increase the accuracy of these, but the Koopman modes become prohibitively expensive to evaluate numerically.

FIG. 4.
Numerical solution u (solid) and Koopman reconstruction (dashed) at t = 0 (black), t = 0.5 (blue), and t = 1 (red) using g = 2, i.e., assuming the dynamics are constrained to a two-dimensional quasiperiodic invariant torus. The Koopman decomposition (23) is truncated after n = 4 with m , m q ∈ { − 3 , … , 3 } 2 ∖ 0.
View largeDownload slide

Numerical solution u (solid) and Koopman reconstruction (dashed) at t = 0 (black), t = 0.5 (blue), and t = 1 (red) using g = 2⁠, i.e., assuming the dynamics are constrained to a two-dimensional quasiperiodic invariant torus. The Koopman decomposition (23) is truncated after n = 4 with m , m q ∈ { − 3 , … , 3 } 2 ∖ 0⁠.

FIG. 4.
Numerical solution u (solid) and Koopman reconstruction (dashed) at t = 0 (black), t = 0.5 (blue), and t = 1 (red) using g = 2, i.e., assuming the dynamics are constrained to a two-dimensional quasiperiodic invariant torus. The Koopman decomposition (23) is truncated after n = 4 with m , m q ∈ { − 3 , … , 3 } 2 ∖ 0.
View largeDownload slide

Numerical solution u (solid) and Koopman reconstruction (dashed) at t = 0 (black), t = 0.5 (blue), and t = 1 (red) using g = 2⁠, i.e., assuming the dynamics are constrained to a two-dimensional quasiperiodic invariant torus. The Koopman decomposition (23) is truncated after n = 4 with m , m q ∈ { − 3 , … , 3 } 2 ∖ 0⁠.

Close modal
FIG. 5.
As for Fig. 4 but with g = 3, and m , m q ∈ { − 3 , … , 3 } 3 ∖ 0.
View largeDownload slide

As for Fig. 4 but with g = 3⁠, and m , m q ∈ { − 3 , … , 3 } 3 ∖ 0⁠.

FIG. 5.
As for Fig. 4 but with g = 3, and m , m q ∈ { − 3 , … , 3 } 3 ∖ 0.
View largeDownload slide

As for Fig. 4 but with g = 3⁠, and m , m q ∈ { − 3 , … , 3 } 3 ∖ 0⁠.

Close modal
Given a discrete time series of snapshots { v j } j from some dynamical system, DMD seeks to find a linear map A such that v j + 1 = A v j⁠. In practice, DMD finds an eigendecomposition of A in which each mode of the decomposition has an associated amplitude a k⁠, spatial pattern u k⁠, and growth rate λ k⁠. As in the case of the Koopman eigenvalues, we expect λ k = i ω k to be purely imaginary because the system is Hamiltonian. We can use our expectation of purely imaginary eigenvalues as a heuristic for a well-resolved mode: if | R e ( λ k ) | ≫ 0⁠, we infer that ω k is an inaccurate guess. Then, if the time between snapshots is τ⁠, we can reconstruct the evolution of the system as
v j = ∑ k a k e i ω k j τ u k .

In the original and most basic form of the algorithm, the number of DMD modes that come from the spectral decomposition will be equal to the spatial dimension. We can increase both the robustness and number of discovered of DMD modes found with delay embedding, so-called Hankel-DMD,31 a higher-order extension in which temporal resolution is substituted for spatial resolution.7 We found much better results when employing this method; we used 20 delays.

As DMD is designed to detect the important temporal frequencies of the dynamics, it should be possible to use it to reconstruct an approximation W ~ for W from a time series. However, as discussed in Sec. IV, the Koopman eigenvalues densely fill the imaginary axis, and therefore, the results of DMD are obscured. For example, if we believe that the solution will be well-represented by an invariant 2-torus—and the DMD eigenfrequencies ω k are sufficiently well-resolved—we expect to see frequencies i ω = m 1 W 1 + m 2 W 2 for many different m 1 , m 2 ∈ Z⁠.

However, due to the exponential term in Eq. (24), we expect the amplitudes to be very low for large values of m k (in absolute terms), and thus, we concentrate on the largest amplitude DMD modes, associating these with small values of m 1 and m 2⁠. We can guess the two fundamental frequencies W ~ 1 and W ~ 2 to be high amplitude modes, and furthermore, the smaller of these should be equal to the gap between the modes. See Fig. 6 for a visual representation of this. We can extend this argument in an obvious way for g > 2⁠.

FIG. 6.
DMD eigenvalues λ from the analysis of a timeseries of u, where points are colored to represent the amplitude | a | of each DMD mode (the colorbar saturates at 10). Note the very small horizontal axis: we plot only eigenvalues with the real part less than 0.1 in magnitude, which excludes a number of spurious modes. High amplitudes are associated with simple integer combinations of the fundamental frequencies, and in particular, we take W ~ 1 = 0.03558 i, W ~ 2 = 1.9312 i, and W ~ 3 = 6.7394 i as the DMD approximations of W—marked with horizontal lines. Prominent eigenvalues are also evident near 2 W 1, 2 W 2, W 1 + W 2, W 3 − W 2, etc.
View largeDownload slide

DMD eigenvalues λ from the analysis of a timeseries of u⁠, where points are colored to represent the amplitude | a | of each DMD mode (the colorbar saturates at 10). Note the very small horizontal axis: we plot only eigenvalues with the real part less than 0.1 in magnitude, which excludes a number of spurious modes. High amplitudes are associated with simple integer combinations of the fundamental frequencies, and in particular, we take W ~ 1 = 0.03558 i⁠, W ~ 2 = 1.9312 i⁠, and W ~ 3 = 6.7394 i as the DMD approximations of W—marked with horizontal lines. Prominent eigenvalues are also evident near 2 W 1⁠, 2 W 2⁠, W 1 + W 2⁠, W 3 − W 2⁠, etc.

FIG. 6.
DMD eigenvalues λ from the analysis of a timeseries of u, where points are colored to represent the amplitude | a | of each DMD mode (the colorbar saturates at 10). Note the very small horizontal axis: we plot only eigenvalues with the real part less than 0.1 in magnitude, which excludes a number of spurious modes. High amplitudes are associated with simple integer combinations of the fundamental frequencies, and in particular, we take W ~ 1 = 0.03558 i, W ~ 2 = 1.9312 i, and W ~ 3 = 6.7394 i as the DMD approximations of W—marked with horizontal lines. Prominent eigenvalues are also evident near 2 W 1, 2 W 2, W 1 + W 2, W 3 − W 2, etc.
View largeDownload slide

DMD eigenvalues λ from the analysis of a timeseries of u⁠, where points are colored to represent the amplitude | a | of each DMD mode (the colorbar saturates at 10). Note the very small horizontal axis: we plot only eigenvalues with the real part less than 0.1 in magnitude, which excludes a number of spurious modes. High amplitudes are associated with simple integer combinations of the fundamental frequencies, and in particular, we take W ~ 1 = 0.03558 i⁠, W ~ 2 = 1.9312 i⁠, and W ~ 3 = 6.7394 i as the DMD approximations of W—marked with horizontal lines. Prominent eigenvalues are also evident near 2 W 1⁠, 2 W 2⁠, W 1 + W 2⁠, W 3 − W 2⁠, etc.

Close modal
We apply DMD to a numerical solution of (1) with initial condition u 0 ( x ) = sin ⁡ x as in Sec. V. Our numerical simulation has a length of 450 time units, with a time resolution τ = 0.1⁠. We find very good agreement in the identification of W as computed analytically to those found in DMD,
W ~ ≈ ( 0.036 i , 1.931 i , 6.739 i ) .
(29)
Note that the differing signs represent degeneracy of the formulation, and these could be recovered in the analytic method by using a different basis of integration contours. Figure 6 shows the DMD eigenvalues as well as their relative amplitudes.
Additionally, given that we know the Riemann theta function form (17), we can exploit the Hirota transform (2). By applying DMD to a time series of ϑ rather than u⁠, we can recover all parameters for the Riemann theta function, i.e., U⁠, W⁠, and B⁠. The details of this procedure are given in the  Appendix. We find
B ~ ≈ ( − 3.077 − 1.907 − 1.907 − 7.150 ) ,
(30)
which is approximately consistent with the values computed analytically, given in (28).

We have performed a Koopman decomposition of the periodic KdV equation. This is almost immediate once the convoluted but well-defined process of periodic inverse scattering is performed. Additionally, we have shown how this result relates to the output of DMD for such a system. DMD gives a very large number of near-imaginary eigenvalues associated with the different harmonics of the nonlinearly interacting waves, which correspond to the Koopman eigenvalues found analytically, which densely fill the imaginary axis.

Furthermore, by exploiting the θ-function representation of the solution, we are able to use DMD to approximately recover the necessary parameters.

We note in passing that since we expect purely imaginary eigenvalues in our system, it is a potential use-case for the physics-informed DMD method32 of finding a unitary matrix A to fit the data. However, we found that this obfuscates the results, as it prevents the use of the real part of the eigenvalue as a measure for how well-resolved a given mode is.

The analytic results of this paper could be extended to other integrable PDEs, which admit Lax pairs, such as the nonlinear Schrödinger equation, the sine-Gordon equation, or the Kadomtsev–Petviashvili equation. The latter could be particularly insightful since it describes two-dimensional wave fields, a significant increase in complexity over the one-dimensional PDEs, which have been studied heretofore.

See the supplementary material that consists of Mathematica notebooks reproducing the results of Secs. IV and V.

This work started life at the Geophysical Fluid Dynamics (GFD) summer school at the Woods Hole Oceanographic Institution, and Sec. VI represents a small part of the fellow’s project of C.V. The authors wish to thank Peter Schmid for his help with this project and everyone at GFD for many fruitful discussions. J.P.P. would like to thank Al Osborne for some useful pointers. J.P.P. was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant No. 865677). C.V. was supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1839302.

The authors have no conflicts to disclose.

Jeremy P. Parker: Conceptualization (lead); Investigation (equal); Methodology (lead); Software (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal). Claire Valva: Conceptualization (supporting); Investigation (equal); Methodology (supporting); Software (equal); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available within the supplementary material.

Here, we give a brief overview of the procedure to recover the parameters for the Riemann theta function, i.e., U⁠, W⁠, and B⁠, from a time series of ϑ⁠.

Again, assuming g = 2⁠, we determine W j from ϑ as we did with u⁠. The wavenumber vector U j is recovered using the same idea; as each DMD frequency ω k is an integer linear combination of i W j⁠, we expect that each DMD spatial mode u k will be a pure sinusoid with wavenumber ℓ j satisfying i ℓ j = n 1 k U 1 + n 2 k U 2⁠. To solve for entries of the period matrix B⁠, we will need three DMD modes, where we can infer n 1 k , n 2 k for each mode k⁠. We construct an invertible matrix M where each row M k has entries ( n 1 k 2 , 2 n 1 k n 2 k , n 2 k 2 )⁠. Then, letting c be a vector such that c k = 2 log ⁡ ( a k )⁠, we solve
M b = c , b = ( B 11 , B 21 , B 22 )
(A1)
for B⁠. We note that for g > 2⁠, we can still determine all parameters of the Riemann theta function give enough well-resolved DMD nodes. However, given that the symmetric g × g matrix B will have ( g + 1 ) g / 2 unique entries, we will need to identify the same number of well-resolved DMD modes, which can be a nontrivial task even for small g⁠.

We apply our DMD analysis to KdV data with initial condition u 0 ( x ) = sin ⁡ x⁠, where we analyze the value of ϑ⁠, rather than u⁠. To three decimal places, we recover the same values for the frequencies (29). Figure 7 shows the eigenvalues, along with the three modes corresponding to frequencies W ~ 1 , W ~ 2⁠, and W ~ 1 + W ~ 2⁠, which were used to determine B ~⁠. The DMD spectrum for ϑ is much cleaner than for u⁠, which shows that the Hirota transform has in some sense simplified the dynamics.

FIG. 7.
Left: As for (6) except that data analyzed were of ϑ rather than u. Horizontal lines show ω j for the well-resolved modes: W ~ 1 = 0.0356 i, W ~ 2 = 1.9312 i, W ~ 1 + W ~ 2 = 1.9668 i, and W ~ 3 = 6.7394 i. Right: DMD spatial modes corresponding to the given frequencies. The real part of the mode is shown in blue and the imaginary part in orange. For each pattern, we also compute the wavenumber ℓ.
View largeDownload slide

Left: As for (6) except that data analyzed were of ϑ rather than u⁠. Horizontal lines show ω j for the well-resolved modes: W ~ 1 = 0.0356 i⁠, W ~ 2 = 1.9312 i⁠, W ~ 1 + W ~ 2 = 1.9668 i⁠, and W ~ 3 = 6.7394 i⁠. Right: DMD spatial modes corresponding to the given frequencies. The real part of the mode is shown in blue and the imaginary part in orange. For each pattern, we also compute the wavenumber ℓ⁠.

FIG. 7.
Left: As for (6) except that data analyzed were of ϑ rather than u. Horizontal lines show ω j for the well-resolved modes: W ~ 1 = 0.0356 i, W ~ 2 = 1.9312 i, W ~ 1 + W ~ 2 = 1.9668 i, and W ~ 3 = 6.7394 i. Right: DMD spatial modes corresponding to the given frequencies. The real part of the mode is shown in blue and the imaginary part in orange. For each pattern, we also compute the wavenumber ℓ.
View largeDownload slide

Left: As for (6) except that data analyzed were of ϑ rather than u⁠. Horizontal lines show ω j for the well-resolved modes: W ~ 1 = 0.0356 i⁠, W ~ 2 = 1.9312 i⁠, W ~ 1 + W ~ 2 = 1.9668 i⁠, and W ~ 3 = 6.7394 i⁠. Right: DMD spatial modes corresponding to the given frequencies. The real part of the mode is shown in blue and the imaginary part in orange. For each pattern, we also compute the wavenumber ℓ⁠.

Close modal
1.
B. O.
Koopman
, “
Hamiltonian systems and transformation in Hilbert space
,”
Proc. Natl. Acad. Sci. U.S.A.
17
(
5
),
315
–
318
(
1931
).
https://doi.org/10.1073/pnas.17.5.315
Google Scholar
Crossref
Search ADS
PubMed
 
2.
I.
Mezić
, “
Spectral properties of dynamical systems, model reduction and decompositions
,”
Nonlinear Dyn.
41
,
309
–
325
(
2005
).
https://doi.org/10.1007/s11071-005-2824-x
Google Scholar
Crossref
Search ADS
 
3.
I.
Mezić
, “
Analysis of fluid flows via spectral properties of the Koopman operator
,”
Annu. Rev. Fluid Mech.
45
,
357
–
378
(
2013
).
https://doi.org/10.1146/annurev-fluid-011212-140652
Google Scholar
Crossref
Search ADS
 
4.
P. J.
Schmid
, “
Dynamic mode decomposition of numerical and experimental data
,”
J. Fluid Mech.
656
,
5
–
28
(
2010
).
https://doi.org/10.1017/S0022112010001217
Google Scholar
Crossref
Search ADS
 
5.
P. J.
Schmid
,
L.
Li
,
M. P.
Juniper
, and
O.
Pust
, “
Applications of the dynamic mode decomposition
,”
Theor. Comput. Fluid Dyn.
25
(
1
),
249
–
259
(
2011
).
https://doi.org/10.1007/s00162-010-0203-9
Google Scholar
Crossref
Search ADS
 
6.
J.
Nathan Kutz
,
S. L.
Brunton
,
B. W.
Brunton
, and
J. L.
Proctor
,
Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems
(
SIAM
,
2016
).
Google Scholar
Crossref
Search ADS
 
7.
P. J.
Schmid
, “
Dynamic mode decomposition and its variants
,”
Annu. Rev. Fluid Mech.
54
(
1
),
225
–
254
(
2022
).
https://doi.org/10.1146/annurev-fluid-030121-015835
Google Scholar
Crossref
Search ADS
 
8.
C. W.
Rowley
,
I.
Mezić
,
S.
Bagheri
,
P.
Schlatter
, and
D. S.
Henningson
, “
Spectral analysis of nonlinear flows
,”
J. Fluid Mech.
641
,
115
–
127
(
2009
).
https://doi.org/10.1017/S0022112009992059
Google Scholar
Crossref
Search ADS
 
9.
J.
Nathan Kutz
,
J. L.
Proctor
, and
S. L.
Brunton
, “
Applied Koopman theory for partial differential equations and data-driven modeling of spatio-temporal systems
,”
Complexity
2018
,
6010634
.
https://doi.org/10.1155/2018/6010634
Crossref
Search ADS
 
10.
J.
Page
 and
R. R.
Kerswell
, “
Koopman analysis of Burgers equation
,”
Phys. Rev. Fluids
3
(
7
),
071901
(
2018
).
https://doi.org/10.1103/PhysRevFluids.3.071901
Google Scholar
Crossref
Search ADS
 
11.
M.
Balabane
,
M.
Alfonso Mendez
, and
S.
Najem
, “
Koopman operator for Burgers’s equation
,”
Phys. Rev. Fluids
6
(
6
),
064401
(
2021
).
https://doi.org/10.1103/PhysRevFluids.6.064401
Google Scholar
Crossref
Search ADS
 
12.
H.
Nakao
 and
I.
Mezić
, “
Spectral analysis of the Koopman operator for partial differential equations
,”
Chaos
30
(
11
),
113131
(
2020
).
https://doi.org/10.1063/5.0011470
Google Scholar
Crossref
Search ADS
PubMed
 
13.
J. P.
Parker
 and
J.
Page
, “
Koopman analysis of isolated fronts and solitons
,”
SIAM J. Appl. Dyn. Syst.
19
(
4
),
2803
–
2828
(
2020
).
https://doi.org/10.1137/19M1305033
Google Scholar
Crossref
Search ADS
 
14.
N. J.
Zabusky
 and
M. D.
Kruskal
, “
Interaction of “solitons” in a collisionless plasma and the recurrence of initial states
,”
Phys. Rev. Lett.
15
(
6
),
240
(
1965
).
https://doi.org/10.1103/PhysRevLett.15.240
Google Scholar
Crossref
Search ADS
 
15.
E.
Fermi
,
P.
Pasta
,
S.
Ulam
, and
M.
Tsingou
, “Studies of the nonlinear problems,” Technical Report, Los Alamos National Lab (LANL), Los Alamos, NM, 1955.
16.
P. D.
Lax
, “
Almost periodic solutions of the KdV equation
,”
SIAM Rev.
18
(
3
),
351
–
375
(
1976
).
https://doi.org/10.1137/1018074
Google Scholar
Crossref
Search ADS
 
17.
E. D.
Belokolos
,
A. I.
Bobenko
,
V. Z.
Enolskii
,
A. R.
Its
, and
V. B.
Matveev
,
Algebro-Geometric Approach to Nonlinear Integrable Equations
(
Springer
,
1994
), Vol. 550.
Google Scholar
 
18.
D.
Johannes Korteweg
 and
G.
de Vries
, “
XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves
,”
Lond. Edinb. Dublin Phil. Mag. J. Sci.
39
(
240
),
422
–
443
(
1895
).
https://doi.org/10.1080/14786449508620739
Google Scholar
Crossref
Search ADS
 
19.
C. S.
Gardner
,
J. M.
Greene
,
M. D.
Kruskal
, and
R. M.
Miura
, “
Method for solving the Korteweg-deVries equation
,”
Phys. Rev. Lett.
19
(
19
),
1095
(
1967
).
https://doi.org/10.1103/PhysRevLett.19.1095
Google Scholar
Crossref
Search ADS
 
20.
D. J.
Benney
, “
Long non-linear waves in fluid flows
,”
J. Math. Phys.
45
(
1–4
),
52
–
63
(
1966
).
https://doi.org/10.1002/sapm196645152
Google Scholar
Crossref
Search ADS
 
21.
D.
Howell Peregrine
, “
Calculations of the development of an undular bore
,”
J. Fluid Mech.
25
(
2
),
321
–
330
(
1966
).
https://doi.org/10.1017/S0022112066001678
Google Scholar
Crossref
Search ADS
 
22.
V.
Iosifovich Karpman
,
Non-Linear Waves in Dispersive Media
, International Series of Monographs in Natural Philosophy (
Elsevier
,
1975
), Vol. 71.
Google Scholar
 
23.
T.
Kappeler
 and
P.
Topalov
, “
Global wellposedness of KdV in H − 1 ( T , R )
,”
Duke Math. J.
135
(
2
),
327
–
360
(
2006
).
https://doi.org/10.1215/S0012-7094-06-13524-X
Google Scholar
Crossref
Search ADS
 
24.
I. C.
Christov
, “
Hidden solitons in the Zabusky–Kruskal experiment: Analysis using the periodic, inverse scattering transform
,”
Math. Comput. Simul.
82
(
6
),
1069
–
1078
(
2012
).
https://doi.org/10.1016/j.matcom.2010.05.021
Google Scholar
Crossref
Search ADS
 
25.
R.
Hirota
,
The Direct Method in Soliton Theory
, Series No. 155 (
Cambridge University Press
,
2004
).
Google Scholar
Crossref
Search ADS
 
26.
S.
Novikov
,
S. V.
Manakov
,
L.
Petrovich Pitaevskii
, and
V.
Evgenevič Zakharov
,
Theory of Solitons: The Inverse Scattering Method
(
Springer Science & Business Media
,
1984
).
Google Scholar
 
27.
A.
Osborne
,
Nonlinear Ocean Waves and the Inverse Scattering Transform
(
Academic Press
,
2010
).
Google Scholar
 
28.
P. D.
Lax
, “
Integrals of nonlinear equations of evolution and solitary waves
,”
Commun. Pure Appl. Math.
21
(
5
),
467
–
490
(
1968
).
https://doi.org/10.1002/cpa.3160210503
Google Scholar
Crossref
Search ADS
 
29.
W.
Magnus
 and
S.
Winkler
,
Hill’s Equation
(
Courier Corporation
,
2013
).
Google Scholar
 
30.
I.
Mezić
, “
Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry
,”
J. Nonlinear Sci.
30
(
5
),
2091
–
2145
(
2020
).
https://doi.org/10.1007/s00332-019-09598-5
Google Scholar
Crossref
Search ADS
 
31.
H.
Arbabi
 and
I.
Mezic
, “
Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the Koopman operator
,”
SIAM J. Appl. Dyn. Syst.
16
(
4
),
2096
–
2126
(
2017
).
https://doi.org/10.1137/17M1125236
Google Scholar
Crossref
Search ADS
 
32.
P. J.
Baddoo
,
B.
Herrmann
,
B. J.
McKeon
,
J.
Nathan Kutz
, and
S. L.
Brunton
, “Physics-informed dynamic mode decomposition (piDMD),” arXiv:2112.04307 (2021).
© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Supplementary Material

Supplementary materials- zip file