Discovering governing equations from data by sparse identification of nonlinear dynamical systems

Chaotic Lorenz System.

As a first example, consider a canonical model for chaotic dynamics, the Lorenz system (40):x˙=σ(y−x),[7a]y˙=x(ρ−z)−y,[7b]z˙=xy−βz.[7c]Although these equations give rise to rich and chaotic dynamics that evolve on an attractor, there are only a few terms in the right-hand side of the equations. Fig. 1 shows a schematic of how data are collected for this example, and how sparse dynamics are identified in a space of possible right-hand-side functions using convex ℓ1 minimization.

For this example, data are collected for the Lorenz system, and stacked into two large data matrices X and X˙, where each row of X is a snapshot of the state x in time, and each row of X˙ is a snapshot of the time derivative of the state x˙ in time. Here, the right-hand-side dynamics are identified in the space of polynomials Θ(X) in (x,y,z) up to fifth order, although other functions such as sin,cos,exp, or higher-order polynomials may be included:Θ(X)=[x(t)||y(t)||z(t)||x(t)2||x(t)y(t)||⋯z(t)5||].Each column of Θ(X) represents a candidate function for the right-hand side of Eq. 1. Because only a few of these terms are active in each row of f, we solve the sparse regression problem in Eq. 3 to determine the sparse vectors of coefficients Ξ=[ξ1ξ2⋯ξn] that determine which terms are active in the dynamics. This is illustrated schematically in Fig. 1, where sparse vectors ξk are found to represent the derivative x˙k as a linear combination of the fewest terms in Θ(X).

In the Lorenz example, the ability to capture dynamics on the attractor is more important than the ability to predict an individual trajectory, because chaos will quickly cause any small variations in initial conditions or model coefficients to diverge exponentially. As shown in Fig. 1, the sparse model accurately reproduces the attractor dynamics from chaotic trajectory measurements. The algorithm not only identifies the correct terms in the dynamics, but it accurately determines the coefficients to within .03% of the true values. We also explore the identification of the dynamics when only noisy state measurements are available (SI Appendix, Fig. S7). The correct dynamics are identified, and the attractor is preserved for surprisingly large noise values. In SI Appendix, section 4.5, we reconstruct the attractor dynamics in the Lorenz system using time-delay coordinates from a single measurement x(t).

PDE for Vortex Shedding Behind an Obstacle.

The Lorenz system is a low-dimensional model of more realistic high-dimensional PDE models for fluid convection in the atmosphere. Many systems of interest are governed by PDEs (24), such as weather and climate, epidemiology, and the power grid, to name a few. Each of these examples is characterized by big data, consisting of large spatially resolved measurements consisting of millions or billions of states and spanning orders of magnitude of scale in both space and time. However, many high-dimensional, real-world systems evolve on a low-dimensional attractor, making the effective dimension much smaller (35).

Here we generalize the SINDy method to an example in fluid dynamics that typifies many of the challenges outlined above. In the context of data from a PDE, our algorithm shares some connections to the dynamic mode decomposition, which is a linear dynamic regression (41⇓–43). Data are collected for the fluid flow past a cylinder at Reynolds number 100 using direct numerical simulations of the 2D Navier–Stokes equations (44, 45). The nonlinear dynamic relationship between the dominant coherent structures is identified from these flow-field measurements with no knowledge of the governing equations.

The flow past a cylinder is a particularly interesting example because of its rich history in fluid mechanics and dynamical systems. It has long been theorized that turbulence is the result of a series of Hopf bifurcations that occur as the flow velocity increases (46), giving rise to more rich and intricate structures in the fluid. After 15 years, the first Hopf bifurcation was discovered in a fluid system, in the transition from a steady laminar wake to laminar periodic vortex shedding at Reynolds number 47 (47, 48). This discovery led to a long-standing debate about how a Hopf bifurcation, with cubic nonlinearity, can be exhibited in a Navier–Stokes fluid with quadratic nonlinearities. After 15 more years, this was resolved using a separation of timescales and a mean-field model (49), shown in Eq. 8. It was shown that coupling between oscillatory modes and the base flow gives rise to a slow manifold (Fig. 2, Left), which results in algebraic terms that approximate cubic nonlinearities on slow timescales.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Example of high-dimensional dynamical system from fluid dynamics. The vortex shedding past a cylinder is a prototypical example that is used for flow control, with relevance to many applications, including drag reduction behind vehicles. The vortex shedding is the result of a Hopf bifurcation. However, because the Navier–Stokes equations have quadratic nonlinearity, it is necessary to use a mean-field model with a separation of timescales, where a fast mean-field deformation is slave to the slow vortex shedding dynamics. The parabolic slow manifold is shown (Left), with the unstable fixed point (C), mean flow (B), and vortex shedding (A). A POD basis and shift mode are used to reduce the dimension of the problem (Middle Right). The identified dynamics closely match the true trajectory in POD coordinates, and most importantly, they capture the quadratic nonlinearity and timescales associated with the mean-field model.

This example provides a compelling test case for the proposed algorithm, because the underlying form of the dynamics took nearly three decades for experts in the community to uncover. Because the state dimension is large, consisting of the fluid state at hundreds of thousands of grid points, it is first necessary to reduce the dimension of the system. The POD (35, 37), provides a low-rank basis resulting in a hierarchy of orthonormal modes that, when truncated, capture the most energy of the original system for the given rank truncation. The first two most energetic POD modes capture a significant portion of the energy, and steady-state vortex shedding is a limit cycle in these coordinates. An additional mode, called the shift mode, is included to capture the transient dynamics connecting the unstable steady state (“C” in Fig. 2) with the mean of the limit cycle (49) (“B” in Fig. 2). These modes define the x,y,z coordinates in Fig. 2.

In the coordinate system described above, the mean-field model for the cylinder dynamics is given by (49)x˙=μx−ωy+Axz,[8a]y˙=ωx+μy+Ayz,[8b]z˙=−λ(z−x2−y2).[8c]If λ is large, so that the z dynamics are fast, then the mean flow rapidly corrects to be on the (slow) manifold z=x2+y2 given by the amplitude of vortex shedding. When substituting this algebraic relationship into Eqs. 8a and 8b, we recover the Hopf normal form on the slow manifold.

With a time history of these three coordinates, the proposed algorithm correctly identifies quadratic nonlinearities and reproduces a parabolic slow manifold. Note that derivative measurements are not available, but are computed from the state variables. Interestingly, when the training data do not include trajectories that originate off of the slow manifold, the algorithm incorrectly identifies cubic nonlinearities, and fails to identify the slow manifold.

Normal Forms, Bifurcations, and Parameterized Systems.

In practice, many real-world systems depend on parameters, and dramatic changes, or bifurcations, may occur when the parameter is varied. The SINDy algorithm is readily extended to encompass these important parameterized systems, allowing for the discovery of normal forms (31, 50) associated with a bifurcation parameter μ. First, we append μ to the dynamics:x˙=f(x;μ),[9a]μ˙=0.[9b]It is then possible to identify f(x;μ) as a sparse combination of functions of x as well as the bifurcation parameter μ.

Identifying parameterized dynamics is shown in two examples: the 1D logistic map with stochastic forcing,xk+1=μxk(1−xk)+ηk,and the 2D Hopf normal form (51),x˙=μx+ωy−Ax(x2+y2)y˙=−ωx+μy−Ay(x2+y2).The logistic map is a classical model for population dynamics, and the Hopf normal form models spontaneous oscillations in chemical reactions, electrical circuits, and fluid instability.

The noisy measurements and the sparse dynamic reconstructions for both examples are shown in Fig. 3. In the logistic map example, the stochastically forced trajectory is sampled at 10 discrete parameter values, shown in red. From these measurements, the correct parameterized dynamics are identified. The parameterization is accurate enough to capture the cascade of bifurcations as μ is increased, resulting in the detailed bifurcation diagram in Fig. 3. Parameters are identified to within .1% of true values (SI Appendix, Appendix C).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

SINDy algorithm is able to identify normal forms and capture bifurcations, as demonstrated on the logistic map (Left) and the Hopf normal form (Right). Noisy data from both systems are used to train models. For the logistic map, a handful of parameter values μ (red lines), are used for the training data, and the correct normal form and bifurcation sequence is captured (below). Noisy data for the Hopf normal form are collected at a few values of μ, and the total variation derivative (33) is used to compute time derivatives. The accurate Hopf normal form is reproduced (below). The nonlinear terms identified by the algorithm are in SI Appendix, section 4.4 and Appendix C.

In the Hopf normal-form example, noisy state measurements from eight parameter values are sampled, with data collected on the blue and red trajectories in Fig. 3 (Top Right). Noise is added to the position measurements to simulate sensor noise, and the total variation regularized derivative (33) is used. In this example, the normal form is correctly identified, resulting in accurate limit cycle amplitudes and growth rates (Bottom Right). The correct identification of a normal form depends critically on the choice of variables and the nonlinear basis functions used for Θ(x). In practice, these choices may be informed by machine learning and data mining, by partial knowledge of the physics, and by expert human intuition.

Similarly, time dependence and external forcing or feedback control u(t) may be added to the vector field:x˙=f(x,u(t),t),t˙=1.Generalizing the SINDy algorithm makes it possible to analyze systems that are externally forced or controlled. For example, the climate is both parameterized (50) and has external forcing, including carbon dioxide and solar radiation. The financial market is another important example with forcing and active feedback control.