arXiv CSopen access

An Introduction to Sparse Identification of Nonlinear Dynamics for Engineering Applications

An Introduction to Sparse Identification of Nonlinear Dynamics for Engineering Applications Yao Cheng Li∗1 , Ana Larrañaga2,3 , Steven L. Brunton2,3 , and Urban Fasel1 1

Department of Aeronautics, Imperial College London, SW7 2AZ, United Kingdom Department of Mechanical Engineering, University of Washington, Seattle, WA 98195, USA 3 NSF AI Institute in Dynamic Systems, University of Washington, Seattle, WA 98195, USA

2

∗

corresponding author: [email protected]

arXiv:2607.15077v1 [cs.LG] 16 Jul 2026

Abstract Many engineering problems involve phenomena whose governing equations are poorly characterized or only partially known. Surrogate modeling techniques such as neural networks can capture the behavior of these systems, but they typically demand large training datasets that are difficult to obtain in engineering contexts and yield models with limited physical interpretability. The Sparse Identification of Nonlinear Dynamics (SINDy) method addresses both limitations by performing sparse regression over libraries of candidate nonlinear terms, recovering interpretable governing equations from comparatively small datasets. Although SINDy has been demonstrated extensively on canonical benchmark systems, its application to practical engineering problems is less widely documented. This tutorial introduces the SINDy method and progressively builds toward its main extensions, from noise-robust weak-form and ensembling-based variants to constrained and parametrizable formulations. The paper and the accompanying tutorial 1 is organized in three parts: the first introduces the standard SINDy algorithm and progressively extends it, inviting readers without prior knowledge to follow each step and adapt the methods to their own problems; the remaining two parts present detailed case studies on (1) the system identification of an unmanned aerial vehicle and (2) a chaotic thermosyphon heat exchanger. Through these examples, we aim to demonstrate that SINDy is simple to implement yet flexible enough to serve as a valuable identification tool for advanced engineering applications.

1

Introduction

Consider the task of modeling the reaction dynamics of a continuous stirred tank reactor using live data, in order to develop a controller that adjusts the heat transfer rate into the tank to control temperature and reactant concentration [1]. Given the strong nonlinearities involved, one could approach this challenge with a neural network, but such models capture input–output correlations rather than the governing equations themselves, offering little insight into the underlying mechanisms [2]. An alternative is to recover the governing dynamics in interpretable form, so that the resulting model itself informs the design decision. This is the motivation behind the Sparse Identification of Nonlinear Dynamics (SINDy) framework [3]. High-stakes engineering decisions require an understanding of a system’s governing physics, not merely a prediction of its behavior. SINDy addresses this need through sparse regression on time-series data, identifying the few nonlinear terms that best characterize the system’s dynamics. Despite its simplicity, the framework has proven effective across a wide range of applications, such as fluid dynamics [4], neuroscience [5], and epidemiology [6]. Researchers have repeatedly shown that SINDy can identify the dynamics of both ordinary differential equations (ODEs) [3] and partial differential equations (PDEs) [7]. Common examples include the Lorenz system and the Van der Pol oscillator for ODEs, and the inviscid Burgers 1

Tutorial available at: https://github.com/paullililili/SINDy4Engineers

1

2 SINDY

Tutorial Structure

DATA

START

ODE PDE

Ensemble Active Learning

1.1

SINDy PDE-Find

correct dynamics?

general recommendations: 1.2

Library Normalization

Embed Prior 1.3 Knowledge

1.8

END

Weak-SINDy

1.7

noisy?

Smoothing

1.5

enough data?

Augment

1.2

right coord?

Manual

1.4

UAV

heat transfer

POD/DMD/AE

1.4

S2

S3

Applied cases

Figure 1: SINDy Tutorial Structure. The repository (hosted in GitHub) is organized into three distinct sections. Section 1 provides a self-contained implementation of all SINDy methods without class dependencies. Sections 2 and 3 present applied examples that use the official SINDy package, PySINDy, to learn the dynamics of an unmanned aerial vehicle and a chaotic flow, respectively. Numbers in the figure indicate the corresponding tutorial notebook in which each component is explained.

equation or a 3D reaction-diffusion model for PDEs. These benchmarks, however, rarely match the conditions found in practice. Measured engineering data are hardly ever clean: sensors introduce noise, sampling is often sparse or irregular, data can come from different sources and vary in quality, acquisition is limited by cost or operational constraints, and the systems of interest are frequently high-dimensional. These limitations are rarely treated from an engineering perspective in the literature, yet, each has an applicable SINDy extension: the weak or integral formulation for noisy and sparsely sampled data [8–11], active or ensemble SINDy for limited or noisy measurements [12, 13], multi-fidelity SINDy for heterogeneous data [14], constrained scenarios [15–17], and its combination with autoencoders for high-dimensional systems [18–21]. This tutorial builds an understanding of SINDy and its main extensions from the ground up (see Fig 1), separating the method itself from the software choices used to implement it [22, 23]. The accompanying notebooks present the mathematics alongside the corresponding code, keeping the two distinct before applying it to several ODE and PDE examples. It then applies these techniques to two engineering problems: (1) an unmanned aerial vehicle and (2) a chaotic thermosyphon flow.

2

Sparse Identification of Nonlinear Dynamics

SINDy is a data-driven method for identifying the governing equations of a system directly from time-series measurements of its state. It applies sparse regression to select a small set of terms from a library of candidate functions, identifying a set of differential equations that describe how the system evolves [3]. Because the resulting model is an explicit equation rather than a blackbox model, it is physically interpretable and tends to extrapolate better to unseen conditions than over-parametrized models such as neural networks. Building from this motivation, suppose that we have a dynamical system defined by a state vector with n variables, x(t) ∈ Rn , such that d x(t) = f (x(t)), dt

(1)

where f is potentially a nonlinear function. In an RC circuit (resistor-capacitor) example, the state can be the voltage decay whose dynamics are governed by dV V (t) = f (V (t)) = − . dt RC

(2)

SINDy assumes that f can be approximated as a linear combination of nonlinear candidate functions θ collected in a library matrix Θ(X) ∈ Rm×l , where X ∈ Rm×n is the collated timeseries data matrix over m measurements to be used in l number of candidate functions (see 2

2.1 Coordinate Basis Selection

2 SINDY

Tutorial 1.1). The library is chosen from prior knowledge of the system, typically starting with simpler terms (such as polynomials of the state variables) and adding complexity as needed (for instance, trigonometric terms). The objective is to find a sparse coefficient matrix Ξ ∈ Rl×n such that Ẋ ≈ Θ(X) Ξ. (3) Higher-degree polynomial terms can make Θ(X) ill-conditioned, meaning that small perturbations in the data (e.g., measurement noise) can produce large changes in the computed coefficients. This is quantified by the condition number κ(Θ): a large value indicates strong amplification of such errors, which in turn increases sensitivity to noise and to the choice of the sparsification threshold (see Tutorial 1.2). It is also good practice to normalize the columns of Θ(X) when the candidate functions have very different magnitudes, as this affects the sparse regression step to solve Ẋ = Θ(X)Ξ, which is the core idea of SINDy. Sparsity is the key assumption: most physical systems are governed by only a few terms, so we seek a solution in which most coefficients are exactly zero. Originally, Brunton et al. [3] proposed sequentially thresholded least squares (STLSQ), a simple approach that first solves an ordinary least squares regression and then zeros out any coefficient falling below a chosen threshold λ (interpreting small coefficients as spurious). The ordinary least squares regression and thresholding steps are then repeated on the remaining terms until a sparse solution is found. The problem is solved for each state variable xi , identifying the sparse coefficient vector ξi (the i-th column of Ξ) by minimizing min ∥Ẋi − Θ(X) ξi ∥2 , ξi

(4)

where Ẋi is the corresponding column of Ẋ. In practice, the sparsification step is applied after the least squares regression by setting all coefficients with a magnitude below λ to zero. If λ is too high, we threshold out the key dynamics, resulting in under-fitting. If λ is too low, we include excess terms that overfit to the dataset. Hence, a right balance must be found. To correctly identify the governing equations, SINDy requires data that sufficiently captures the behavior of the system. Trajectories that are too short, sampled too infrequently, or initialized from a limited range of initial conditions may fail to resolve the full dynamics and lead to an incorrect identification. When data is limited, incorporating prior physical knowledge can help guide the sparse regression toward the correct equations, which can be done for example by constraining the candidate function library [15, 16, 24], or by selecting a model based on coefficient stability [25] (see Tutorial 1.3).

2.1

Coordinate Basis Selection

In many practical applications, the observed data collected from experimentation or numerical simulations may not directly measure the state(s) of the system. For instance, if we track an undamped pendulum using video cameras, the motion-tracking software may return the xy position of the mass rather than its angular displacement. Learning the system dynamics using xy coordinates would require learning the rotation transformation simultaneously, whereas directly learning from the angular states would yield the critical and yet sparse governing equation of θ̈ = −g sin(θ) (see Tutorial 1.4). If a system’s intrinsic coordinate basis is known, coordinate transformation should always be performed as part of data pre-processing. Whilst SINDy excels at discovering dynamics in low-dimensional problems, its candidate function library scales combinatorially with the number of states, making it computationally intractable for high-dimensional systems. In such cases, SINDy can be combined with dimensionalityreduction methods such as Proper Orthogonal Decomposition (POD) or Dynamic Mode Decomposition (DMD), allowing it to discover dynamics in a lower-dimensional latent space. For 3

2.2 Noise Robustness

2 SINDY

instance, the flow field around a cylinder can be decomposed into several wake shedding modes using POD [3], after which SINDy can identify the governing equations of the reduced system. In POD, measurements at each snapshot in time x(tj ) ∈ Rm are stacked as columns of a data matrix X ∈ Rm×n , where each column is one of n snapshots of the m-dimensional state such that X = [x(t1 ) x(t2 ) . . . x(tn )]. It extracts a set of orthogonal spatial modes from this matrix using the singular value decomposition (SVD) [26] X = U ΣV T ,

(5)

where the columns of the orthonormal matrix U ∈ Rm×m are the spatial modes, the diagonal matrix Σ ∈ Rm×n holds the singular values σi that rank each mode by its amplitude, and the orthonormal rows of V T ∈ Rn×n describe the temporal evolution of the modes 2 . The POD modes are arranged in descending amplitude, from most to least energetic, in the column space of U . The first r columns Ur can then be used to project the data onto a lower-dimensional space via Z = UrT X. Here, Z ∈ Rr×n contains the temporal evolution of the r POD modes, which can then be used with SINDy to discover the dynamics in the latent space (see Tutorial 1.4). Unlike POD, DMD computes a linear mapping matrix A that relates the snapshot matrix X = [x(t1 ) x(t2 ) . . . x(tn−1 )] to its time-shifted counterpart Y = [x(t2 ) x(t3 ) . . . x(tn )] via Y = AX. Using reduced-rank regression to compute A [27], the eigenvectors Ψ ∈ Rm×r of A can be used to project the data to an r-dimensional space Z = Ψ† X, which can again be used with SINDy to discover dynamics in a low-dimensional latent space. POD and DMD are both linear dimensionality reduction methods that often require many modes to capture the dominant physics. Autoencoders can be considered a nonlinear extension of these methods, compressing the same amount of variance (or energy in the fluid flow example) into significantly fewer latent variables [28]. Specifically, an autoencoder is a type of neural network that learns a nonlinear mapping through an encoder φ : Rm → Rr and a decoder ψ : Rr → Rm where z(tj ) = φ(x(tj )) and x(tj ) = ψ(z(tj )), enabling it to capture more complex structures with fewer latent variables than linear methods. Combined with SINDy, as first demonstrated by Champion et al. [18], this approach enables the simultaneous learning of the dimensionality reduction mapping and the sparse governing dynamics within the latent space (see Tutorial 1.4).

2.2

Noise Robustness

Obtaining accurate time derivatives is one of the main challenges when using SINDy, as they are rarely measured directly and must instead be computed numerically using techniques such as finite differences (see Tutorial 1.1). When the data is noisy, however, direct numerical differentiation becomes unreliable, as illustrated by the following example. Let us define y = f (x) + ϵ(σ)

(6)

where ϵ is an additive Gaussian noise component with a standard deviation of σ. Specifically, let us define ϵ as a single sine wave with frequency ωnoise and amplitude A such that y = f (x) + A sin(ωnoise x), ( n n (−1) 2 Aωnoise sin(ωnoise x) for n even dn y (n) = f (x) + n−1 n dxn 2 Aωnoise cos(ωnoise x) for n odd (−1)

(7a) (7b)

Therefore, the higher the frequency, the larger the amplitude of the noise component in the derivative, which is further amplified by using higher-order derivatives. Some methods can 2

In this tutorial, we assume that typical engineering data is real-valued. A more general SVD formulation uses the Hermitian transpose V H for complex-valued data, or adjoint V ∗ for non-standard inner products.

4

2.3 Controller Dynamics Identification

2 SINDY

mitigate this, such as data smoothing algorithms and differentiation with total variation regularization [29] (see Tutorial 1.5), as well as weak-form SINDy formulations [8, 9] that bypass point-wise differentiation altogether (see Tutorials 1.7.1 and 1.7.2).

2.3

Controller Dynamics Identification

A natural extension of SINDy is SINDy with Controls (SINDyc) [30], which enables the discovery of parameterized systems or systems subjected to forcing or control inputs u(t) ∈ Rq such that d x(t) = f (x(t), u(t)). dt

(8)

In SINDyc, the candidate function library Θ(X, U ) is extended to include controls U . For a polynomial based library, the candidate functions are x, y, z, u, x2 , xy, . . . , zu2 , u3 . Once constructed, the exact same sparse regression problem is repeated to solve Ẋ ≈ Θ(X, U )Ξ.

(9)

The learned model can for example be used with a receding horizon model predictive controller (MPC), which is an optimal control strategy that uses a model of the system to predict the evolution of the state over a finite horizon [31] (see Tutorial 1.6). Alternatively, SINDy can be combined with reinforcement learning (RL), a framework in which an agent learns a control policy through trial-and-error interaction with its environment by maximizing a cumulative reward. Standard deep RL is notoriously data-hungry and yields black-box policies that are not interpretable. SINDy-RL instead learns sparse, interpretable surrogate models of the system dynamics, the reward function, and the control policy [32]. By replacing expensive environment interactions with cheap evaluations of the learned dynamics model, it improves sample efficiency whilst obtaining compact, interpretable policies.

2.4

Ensembling & Active Learning

Finally, uncertainty quantification is a desirable capability, especially when operating in the low data regime. Ensembling methods can tackle this, combining multiple models trained on different subsets of the data, most commonly through bootstrap aggregation (bagging). In the SINDy setting, this is known as Ensemble-SINDy (E-SINDy) [13], which bootstraps from the data samples and, optionally, from the candidate library terms. Each of the m ensemble members is then fitted independently,
Ẋ (k) ≈ Θ X (k) Ξ(k) , k = 1, . . . , m, (10) and the resulting coefficients are aggregated by their mean to give a bagged estimate, while the median yields a more outlier-robust variant (bragging), m

1 X (k) Ξ̄ = Ξ m

Ξ̃ = mediank Ξ(k)

(bagging),

(bragging).

(11)

k=1

The model ensemble can provide an uncertainty estimate for both the predictions and the coefficient identification, making it possible to track which coefficients appear most consistently across the ensemble (see Tutorial 1.8.1). In fact, such uncertainty estimates can be used as a sampling criterion in active learning scenarios, where gathering new data is expensive and each sampled time-series must be used efficiently [12]. Rather than sampling long periods of time, short trajectories are distributed across different regions of the state space, enabling new dynamics to be discovered at a lower data cost (see Tutorial 1.8.2). 5

3 UAV SYSTEM IDENTIFICATION

3

ODE Case Study: UAV System Identification

In recent years, unmanned aerial systems have seen increased deployment across a wide range of applications, from cargo and personnel transportation to monitoring, surveillance, mapping, and inspection, establishing themselves as a cornerstone enabling technology [33]. This calls for rapid innovation and iterations to seek out unique design solutions that address various requirements. For instance, a hybrid fixed-wing vertical takeoff and landing (VTOL)-capable unmanned aerial vehicle (UAV) exploits both the long range and endurance of a fixed wing, while offering the flexible operational capabilities of a rotorcraft vehicle [34]. However, the dynamics of such vehicles are challenging to model due to effects such as wing-propeller interactions or vertical to forward flight transitions. Furthermore, UAVs have become widely accessible and can be rapidly prototyped and tested, making flight test data increasingly available compared to high-fidelity computational simulations or wind tunnel experiments. This accessibility motivates the use of data-driven methods to identify UAV dynamics purely from time-series measurements, or to correct existing low to medium-fidelity models using system identification methods. The common approach to UAV system identification involves parametric estimation using a known analytical model [35–37]. Recent developments in deep learning have also spurred the application of black-box neural networks in the identification process [38, 39] to circumvent the need for known models. In our tutorial, we present a constrained SINDy approach that enforces known equation structures and coefficients, whilst still allowing SINDy to learn additional terms that describe physics without known priors (see Tutorial 2). The use of constrained regression, noisy data smoothing, and custom function library design methods introduced in the previous tutorials will be applied here.

3.1

Quadcopter Model and Simulation

In this tutorial, we use a quadcopter UAV configuration, as shown in Figure 2, propelled by four motors and propellers in a cross layout. All four motors contribute equally to the vertical thrust of the UAV, and differential thrust contributes to the rolling moments Mx or pitching moments My . Opposing pairs of motors rotate in the same direction to generate equal motor torque τm , and differential control pairs contribute to yawing moments Mz . The full six degrees of freedom system defining a quadcopter’s kinematics, defined by its translational velocities  T  T and angular velocities B ω = p q r . Whilst most applied forces and B v = vx vy vz moments, such as those generated by thrust inputs, are oriented in the body frame of reference (denoted by B (·)), gravitational force is oriented in the inertial frame of reference. As a result, the governing equation contains rotational matrices and cross product terms from rigid body kinematics shown in Figure 2. The full expanded dynamics of the UAV system is not easily discoverable with the standard SINDy implementation and are defined as ϕ̇ = p + tan(θ) [q sin(θ) + r cos(ϕ)] ,

(12a)

θ̇ = q cos(ϕ) − r sin(ϕ),

(12b)

ψ̇ = sec(θ) [q sin(θ) + r cos(ϕ)] ,

(12c)

v̇x = (rvy − qvz ) − g sin(θ) + Fx /m,

(12d)

v̇y = (pvz − rvx ) + g cos(θ) sin(ϕ) + Fy /m,

(12e)

v̇z = (qvx − pvy ) + g cos(ϕ) cos(θ) + Fz /m,

(12f)

ṗ = [Mx + qr(Iyy − Izz )] /Ixx ,

(12g)

q̇ = [My + pr(Izz − Ixx )] /Iyy ,

(12h)

ṙ = [Mz + pq(Ixx − Iyy )] /Izz .

(12i)

We can omit the position states of the UAV for this tutorial, since they do not affect the 6

3.1 Quadcopter Model and Simulation

3 UAV SYSTEM IDENTIFICATION

Quadcopter

Partially Known Dynamics

𝜔!!

𝜏! !

𝑚 !𝒗̇ + !𝝎 × !𝒗 = !𝑭 𝑱 !𝝎̇ + !𝝎 × 𝑱 !𝝎 = !𝑴

𝜏!" 𝜔!" "𝑥

SINDy System Identification Process Data smoothing

𝑦 "𝑧 " 𝜔!# 𝜏!#

𝜏! $

𝜔!$ Constrained sparse regression

Noisy Sensor Measurements

| 𝜙̇ |

| 𝜃̇ |

| 𝑟𝑣 = … " |

𝑋̇

subject to CΞ ∶ = d

| 𝑞𝑣# | 𝛩(𝑋)

…

…

Ξ

Model Predictions

Figure 2: Schematic of the UAV system identification process. We first extract sensor measurements from simulated flight test data before applying smoothing techniques. A constrained sparse regression approach is then taken using partial knowledge of the model to identify the governing equations and predict future trajectories with it.

dynamics when ground effect is not present and wind is negligible. Whilst complex, these functions are known from rigid body kinematics, and only the applied forces F and moments M remain unknown. We can therefore use a constrained regression approach with constrained SR3 algorithm [16] to allow SINDy to discover the unknown forces and moments, while retaining known components of the dynamics. In the tutorial, we explore two separate approaches: the construction of a tailored candidate function library to include all non-externally applied force and moment functions in Equations 12, and second, applying linear equality constraints on known coefficients in Equations 12. We generate the training and validation data using Mathwork’s Parrot Minidrone simulation in Matlab 2025b, which is a medium-fidelity simulation that models certain physics such as blade flapping, pitch and roll rotor damping, and body drag, but ignores higher-fidelity physics such as unsteady aerodynamics with time-scales too small to be significant here [36, 40, 41]. Furthermore, it simulates onboard Kalman filter-based state estimators, which use sensor data from a combination of simulated IMU sensors and a mimicked visual odometry sensor for position and velocity drift corrections. Furthermore, an onboard cascaded-PID controller controls the motor RPM to match reference position states and yaw angle. To obtain data that is sufficiently rich in dynamics that SINDy may learn from, we generate a series of decoupled and coupled step and chirp reference signals that the controller must follow. Finally, we collate all time-series

7

3.2 Constrained SINDy Learning

3 UAV SYSTEM IDENTIFICATION

state measurements into a data matrix X ∈ Rm×n and control matrix U ∈ Rm×q .

3.2

Constrained SINDy Learning

Since many of the known functions, specifically ones that originate from the application of the rotation matrices, within Equations 12 cannot be built from simple basis function libraries such as polynomial or Fourier libraries, we define custom functions using the appropriate terms such that  T p + tan(θ) (q sin(ϕ) + r cos(ϕ))   q cos(ϕ) − r sin(ϕ)     sec(θ) (q sin(ϕ) + r cos(ϕ))     − sin(θ) Θcustom (X) =  (13)  .   cos(θ) sin(ϕ)     cos(ϕ) cos(θ)   .. . In addition to the known terms present in the Equations 12, we include additional linear velocity squared terms, and pitch and roll rates to allow SINDy to learn the effects of body drag and rotor damping, respectively. We can further introduce a 3rd order polynomial basis library   for control inputs only such that Θpolynomial (U ) = u1 (t) u2 (t) u3 (t) u4 (t) . . . (u4 (t))3 where ui (t) corresponds to the time-series input for the ith motor. This library allows SINDy to discover the control terms present in the system. Finally, we concatenate the libraries together such that   Θ(X, U ) = Θcustom (X) Θpolynomial (U ) , (14) where the library matrix Θ(X, U ) ∈ Rm×l contains l candidate functions. Next, to enforce sparsity or enforce known coefficients as shown in Section 2, we can constrain known functions and equation structure by creating linear constraints CΞ(:) = d. In the tutorial, we programmatically create the constraints. For instance, to enforce unity on p + tan(θ) (q sin(ϕ) + r cos(ϕ)) for ϕ̇, we set our first constraint as ( 1 if j = 1, C1,j = (15) 0 if j ̸= 1, with d1 = 1. We must also enforce sparsity on this function for all other states such that ( 1 if j = (i − 1)l + 1, Ci,j = (16) 0 if j ̸= (i − 1)l + 1, and di = 0 for i ∈ [2, n]. The same approach can be applied to other custom functions, such as functions with gravitational acceleration constant g, to enforce both coefficient and sparsity. In some cases, however, we omit the application of Equation 15 to allow SINDy to learn the coefficients. For instance, for the cross-coupling terms of angular rates qr, pr, and pq, we allow SINDy to learn the coefficients which are a function of the UAV’s moments of inertia Ixx , Iyy , and Izz . Further constraints can also be enforced, such as control input symmetry. For instance, the forward and rear sets of motors should have equal and opposite pitching moments. In this tutorial, the combination of defining custom candidate functions and constraining known coefficients and sparsity yields an interpretable model that matches the expected equation structure and produces accurate predictions despite the noisy data, as shown in Figure 2. While demonstrated here on a UAV dataset, the approach generalizes naturally to other system identification problems with partially known coefficients and equation structure, and we encourage the reader to test it on their own engineering datasets. For readers interested specifically in UAVs, there are also other recent approaches to applying SINDy in UAV system identification. For example, Osman et al. [42] who presented an application of an adaptive SINDy framework 8

4 CHAOTIC THERMOSYPHON FLOW

Discover governing PDEs

𝑇 𝜋/2, 𝑅" = 𝑇!

| 𝜔̇ |

r

| | 𝑇̇ = 𝒖 + 𝛻𝜔 | | 𝑈̇

𝜃

| 𝛻"𝜔 |

…

𝛩(𝑈)

… Ξ

Discover low-dimensional ODEs

g 𝑅!

𝑅"

𝑇 −𝜋/2, 𝑅" = 𝑇"

Figure 3: A thermosyphon contains convection-driven flow within its concentric inner and outer walls, subjected to a Dirichlet temperature boundary condition on its outer wall. Its inner wall is adiabatic. The governing PDEs may be discovered using PDE-FIND with a differential operator library. We may also discover a lower-dimensional latent space attractor (shown in normalized coordinates) using a fully data-driven approach that combines SINDy with DMD or using an "informed" approach where we know the variables to measure or reduce to.

integrated with a Lyapunov-based Model Predictive Controller to produce a real-time controller for a VTOL-capable flying wing. Manaa et al. [43] have also presented an application of SINDy to learn the governing dynamics of a quadrotor system constrained to three degrees of freedom only, using a combination of polynomial and Fourier libraries.

4

PDE Case Study: Chaotic Thermosyphon Flow

In this section of the tutorial, we apply SINDy to identify the governing PDEs of a flow within a thermosyphon shown in Figure 3, and learn a lower-dimensional representation that can be used to analyze and predict the flow’s evolution. PDEs present unique challenges to the identification process. For example, directly adding partial derivatives into the candidate function library like in PDE-FIND [7] results in high degrees of collinearity between functions. Furthermore, SINDy is primarily designed to discover equations where the time derivatives of its states can be explicitly defined, and its terms are spatially invariant. If such criteria are not met, such as in this tutorial (see Equations 17), we have to develop alternative workarounds. Additionally, we are often interested in identifying a nonlinear dynamical system that lies within a lowerdimension latent space from the spatiotemporal solution, which is useful to analyze and predict future evolution of the solution fields. A thermosyphon is a passive heat exchanger that relies on convection-driven flow to transport heat around the geometry. Early studies found that self-excited oscillations can occur in closedloop thermosyphon geometries due to the interplay between buoyant and viscous forces [44, 45]. This was later shown to have similar bifurcation behaviors as the chaotic Lorenz system [46]. More recently, Huang et al. [47] derived a low-dimensional model from Naiver-Stokes Boussinesq equations, which treats the flow as incompressible but assumes that its temperature T varies linearly with its density, using a truncated Fourier-Laurent expansion. Loiseau [27] instead implemented a rank-constrained DMD to reduce the solution’s dimensionality before applying constrained SINDy to identify a Lorenz like system. The low-dimensional time-series trajectories from both approaches are shown in Figure 3. In this tutorial, we follow the flow setup in [47]

9

4 CHAOTIC THERMOSYPHON FLOW

to demonstrate that the derived analytical system can be discovered purely from data using SINDy, as well as demonstrating an alternative approach to the fully data-driven DMD method presented in [27]. The annulus geometry of the thermosyphon is presented in Figure 3. The heat exchanger is subjected to a temperature Dirichlet boundary condition on its outer walls that varies linearly with the vertical height, and the inner wall is adiabatic. The Navier-Stokes Boussinesq equations can be reformulated into the divergence-free vorticity-streamfunction equations   ∂ω ∂T 1 ∂T 2 + u · ∇ω = Pr∇ ω + Pr Ra cos θ − sin θ , (17a) ∂t ∂r r ∂θ ∂T + u · ∇T = ∇2 T, (17b) ∂t −∇2 ψ = ω, (17c) which are parameterized with the Rayleigh number Ra, which defines the ratio between buoyant and viscous forces, and Prandtl number Pr, which defines the ratio between momentum and viscous diffusivity. At low Rayleigh numbers, any transients in the flow quickly decay away and the flow returns to a zero-circulation static state. Past a critical Rayleigh number, the flow undergoes a pitchfork bifurcation and settles into either a clockwise or counter-clockwise constant circulation. At even higher Rayleigh numbers, the flow transitions into a chaotic regime and exhibits aperiodic flow reversal behavior. The ability to predict such bifurcation points is critical to any engineer designing such a heat exchanger (see Tutorial 3.1). In an "informed" approach, we can extract low-dimensional states from the spatiotemporal solution that capture the key dynamics. Huang et al. [47] suggested extracting Z 2π Z R1 1 X= r2 udrdθ, A0 0 R2 Z 1 xT dA, Y =− A0 Ω Z 1 Z=− yT dA. A0 Ω

(18a) (18b) (18c)

which are the average angular momentum of the flow X, flow horizontal center of mass Y , and flow vertical center of mass Z, respectively, all normalized by the area of the annulus A0 . An example trajectory extracted from a flow with Ra = 1.0 × 107 and P r = 4.0 is shown in Figure 3, which shares a similar attractor structure with the Lorenz system. Using these measured variables, a set of parameterized equations can be derived in the form of Ẋ = −RaPrY − αPrX,

(19a)

Ẏ = −kX(Z − z1 ) − βY,

(19b)

Ż = kXY − β(Z − z0 ),

(19c)

where the constants α, β, k, z0 , and z1 are all functions of the radii of the thermosyphon [47]. In the tutorial, we show that SINDyc can be used to learn this parameterized set of equations directly from data, using a cubic polynomial library and STLSQ. The predicted derivatives learned by SINDyc produce an even more accurate model than the "informed" approached, shown in Figure 4. Furthermore, we also show how the nonlinear system learned by SINDy can be ported to a symbolic Maths package such as SymPy [48], where a nonlinear system analysis can be carried out to accurately identify the critical Rayleigh numbers where bifurcation occurs (see Tutorial 3.4).

10

4 CHAOTIC THERMOSYPHON FLOW

4e6

1e7

2e6

0

0

0

-1e7 200 0

50

200

50

50

0

0

Z

50

Z

Z

Y

0

Y

0

X

X

2e6

Y

X

Increasing Rayleigh number

0 0.00

0.01

Time (s)

0.02

200 0.00

Learned system from data

0.01

Time (s)

0.02

Derived system from (Huang et al, 2023)

0.00

0.01

Time (s)

0.02

Data from simulated system

Figure 4: A comparison of the state derivative predictions between the learned SINDy model and derived model against simulated data for increasing Rayleigh number.

The tutorial also demonstrates how the rank-constrained DMD from [27] can be applied to discover the underlying chaotic attractor dynamics. In this approach, we can extract one DMD mode from the azimuthal velocity field and two DMD modes from the temperature field to obtain the low-dimensional system shown in Figure 4. Whilst this coordinate transformation with DMD results in a different time-series dataset compared to the "informed" approach shown in 3, they both share the same chaotic attractor dynamics (see Tutorial 3.5). Lastly, the tutorial explores how SINDy can be applied to a challenging PDE identification problem, learning Equation 17 directly from data. This problem is particularly difficult as the time derivatives cannot be explicitly defined due to the Poisson equation in Equation 17c), the system is spatially variant, and the equations are expressed in polar coordinates. Similar to [49], we can define a candidate function library with differential operators, such as the Laplacian operator ∇2 (·). With this library, we can learn the governing PDE in the form

∂ω 1 ∂ψ ∂ω 1 ∂ψ ∂ω + − = Pr ∂t |r ∂θ ∂r {z r ∂r ∂θ} u·∇ω

   ∂T ∂ 2 ω 1 ∂ω 1 ∂2ω 1 ∂T + + 2 2 +Pr Ra cos θ − sin θ , ∂r2 r ∂r r ∂θ ∂r r ∂θ {z } | 

(20a)

∇2 ω

1 ∂ψ ∂T 1 ∂ψ ∂T ∂2T 1 ∂T 1 ∂2T ∂T + − = + + 2 2, 2 ∂t r ∂r |r ∂θ ∂r {z r ∂r ∂θ} |∂r {z r ∂θ }

(20b)

∇2 T

u·∇T

 ∂ 2 ψ 1 ∂ψ 1 ∂2ψ ω=− + + . ∂r2 r ∂r r2 ∂θ2 | {z } 

(20c)

∇2 ψ

We show that using such a library eliminates the collinearity that arises from learning individual partial derivatives, and we show that the true equation using a non-weak form based method can be accurately recovered (see Tutorial 3.2).

11

REFERENCES

5

Conclusion

In this tutorial paper and its accompanying Github repository, we introduce the SINDy method and show how SINDy can be implemented, extended, and applied to relevant engineering applications. This ranges from learning parametrizable systems with actuation that can be integrated with model-based controllers [30, 31], to formulating noise robust extensions such as weak form SINDy [8,9]. The tutorial also applies SINDy and its extensions to two engineering applications: a UAV system identification and a chaotic heat exchanger modeling problem. We showcase how commonly faced challenges can be mitigated, for example, how noisy data can initially be smoothed, or how sparsity threshold hyperparameters may be selected using statistical stability criteria [25]. Whilst we have covered several key SINDy extensions, many more exist that tackle other challenges. For example, Kaptanoglu et al. [17] developed the Trapping SINDy optimizer, which promotes the discovery of systems with bounded trajectories. An alternative line of work adopts a Bayesian approach to the regression problem, performing thresholding based on coefficient likelihood rather than magnitude [50–52]. Separately, Boninsegna et al. [53] adapted SINDy to discover stochastic differential equations. We recommend readers to check out the PySINDy repository [22, 23], which contains an extensive list of readily available SINDy implementations in Python.

Acknowledgments This work was jointly funded by UK Research & Innovation’s Industrial Doctoral Landscape Awards (reference number UKRI1983) with Airbus UK, and the National Science Foundation AI Institute in Dynamic Systems (grant number 2112085).

References [1] Bhavana Bhadriraju, Abhinav Narasingam, and Joseph Sang-Il Kwon. Machine learningbased adaptive model identification of systems: Application to a chemical process. Chemical Engineering Research and Design, 2019. [2] Cynthia Rudin. Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. Nature Machine Intelligence, 2019. [3] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 2016. [4] Kai Fukami, Takaaki Murata, Kai Zhang, and Koji Fukagata. Sparse identification of nonlinear dynamics with low-dimensionalized flow representations. Journal of Fluid Mechanics, 2021. [5] Robin Delabays, Giulia De Pasquale, Florian Dörfler, and Yuanzhao Zhang. Hypergraph reconstruction from dynamics. Nature Communications, 2025. [6] Jonathan Horrocks and Chris T. Bauch. Algorithmic discovery of dynamic models from infectious disease data. Scientific Reports, 2020. [7] Samuel H. Rudy, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 2017. [8] Daniel A. Messenger and David M. Bortz. Weak SINDy: Galerkin-Based Data-Driven Model Selection. Multiscale Modeling & Simulation, 2021. 12

REFERENCES

REFERENCES

[9] Daniel A. Messenger and David M. Bortz. Weak SINDy for partial differential equations. Journal of Computational Physics, 2021. [10] Hayden Schaeffer and Scott G. McCalla. Sparse model selection via integral terms. Physical Review E, 2017. [11] Patrick A. K. Reinbold, Daniel R. Gurevich, and Roman O. Grigoriev. Using noisy or incomplete data to discover models of spatiotemporal dynamics. Physical Review E, 2020. [12] Ana Larrañaga, Urban Fasel, and Steven L. Brunton. How low can you go? active learning for sparse model discovery in the ultra-low-data limit, 2026. [13] Urban Fasel, J. Nathan Kutz, Bingni W. Brunton, and Steven L. Brunton. EnsembleSINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2022. [14] Filippo Zacchei, Ana Larrañaga, Attilio Frangi, Andrea Manzoni, and Steven L. Brunton. Multi-fidelity sindy: Sparse discovery of nonlinear dynamical systems with fidelity-weighted measurements, 2026. [15] J.-C. Loiseau and S. L. Brunton. Constrained sparse Galerkin regression. Journal of Fluid Mechanics, 2018. [16] Kathleen Champion, Peng Zheng, Aleksandr Y. Aravkin, Steven L. Brunton, and J. Nathan Kutz. A unified sparse optimization framework to learn parsimonious physics-informed models from data. IEEE Access, 2020. [17] Alan A Kaptanoglu, Jared L Callaham, Christopher J Hansen, Aleksandr Aravkin, and Steven L Brunton. Promoting global stability in data-driven models of quadratic nonlinear dynamics. Physical Review Fluids, 2021. [18] Kathleen Champion, Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 2019. [19] Joseph Bakarji, Kathleen Champion, J. Nathan Kutz, and Steven L. Brunton. Discovering governing equations from partial measurements with deep delay autoencoders. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2023. [20] L Mars Gao and J Nathan Kutz. Bayesian autoencoders for data-driven discovery of coordinates, governing equations and fundamental constants. Proceedings of the Royal Society A, 2024. [21] Paolo Conti, Jonas Kneifl, Andrea Manzoni, Attilio Frangi, Jörg Fehr, Steven L Brunton, and J Nathan Kutz. VENI, VINDy, VICI: a generative reduced-order modeling framework with uncertainty quantification. Neural Networks, 2026. [22] Brian M de Silva, Kathleen Champion, Markus Quade, Jean-Christophe Loiseau, J Nathan Kutz, and Steven L Brunton. Pysindy: A python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 2020. [23] Alan A Kaptanoglu, Brian M de Silva, Urban Fasel, Kadierdan Kaheman, Jared L Callaham, Charles B Delahunt, Kathleen Champion, Zachary G. Nicolaou, Kathleen Champion, Jean-Christophe Loiseau, J Nathan Kutz, and Steven L Brunton. Pysindy: A comprehensive python package for robust sparse system identification. Journal of Open Source Software, 2022.

13

REFERENCES

REFERENCES

[24] Peng Zheng, Travis Askham, Steven L Brunton, J Nathan Kutz, and Aleksandr Y Aravkin. Sparse relaxed regularized regression: SR3. IEEE Access, 2019. [25] Suryanarayana Maddu, Bevan L. Cheeseman, Ivo F. Sbalzarini, and Christian L. Müller. Stability selection enables robust learning of differential equations from limited noisy data. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2022. [26] Julien Weiss. A tutorial on the proper orthogonal decomposition. In AIAA aviation 2019 forum, 2019. [27] Jean-Christophe Loiseau. Data-driven modeling of the chaotic thermal convection in an annular thermosyphon. Theoretical and Computational Fluid Dynamics, 2020. [28] Kookjin Lee and Kevin T Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. Journal of Computational Physics, 2020. [29] Rick Chartrand. Numerical Differentiation of Noisy, Nonsmooth Data. ISRN Applied Mathematics, 2011. [30] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Sparse identification of nonlinear dynamics with control (SINDYc). IFAC-PapersOnLine, 2016. [31] E. Kaiser, J. N. Kutz, and S. L. Brunton. Sparse identification of nonlinear dynamics for model predictive control in the low-data limit. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2018. [32] Nicholas Zolman, Christian Lagemann, Urban Fasel, J Nathan Kutz, and Steven L Brunton. SINDy-RL for interpretable and efficient model-based reinforcement learning. Nature Communications, 2025. [33] Anuradha M. Annaswamy, Karl H. Johansson, and George Pappas. Control for societalscale challenges: Road map 2030. IEEE Control Systems Magazine, 2024. [34] Carl Olsson, Sebastian L. Verling, Thomas Stastny, and Roland Siegwart. Full Envelope System Identification of a VTOL Tailsitter UAV. In AIAA Scitech 2021 Forum, VIRTUAL EVENT, 2021. [35] Mehmet Baskın and Kemal Leblebicioğlu. Frequency-domain estimation of a transfer matrix of an uncommon quadrotor in hover. IEEE Transactions on Control Systems Technology, 2022. [36] Paul Pounds, Robert Mahony, and Peter Corke. System identification and control of an aerobot drive system. In 2007 Information, Decision and Control. IEEE, 2007. [37] Benjamin M. Simmons, James L. Gresham, and Craig A. Woolsey. Nonlinear Dynamic Model Identification for Aircraft with Unknown Mass Properties Using Flight Data. In AIAA SCITECH 2023 Forum, National Harbor, MD & Online, 2023. [38] Joshua Harris, Frank Arthurs, James V. Henrickson, and John Valasek. Aircraft system identification using artificial neural networks with flight test data. In 2016 international conference on unmanned aircraft systems (ICUAS). IEEE, 2016. [39] Mohammad Fahmi Pairan, Syariful Syafiq Shamsudin, and Mohd Fadhli Zulkafli. Neural network based system identification for quadcopter dynamic modelling: A review. Journal of Advanced Mechanical Engineering Applications, 2020. [40] Paul Pounds, Robert Mahony, and Peter Ian Corke. Modelling and control of a quad-rotor robot. In Proceedings of the 2006 Australasian Conference on Robotics and Automation, 2006.

14

REFERENCES

REFERENCES

[41] Fabian Riether. Agile quadrotor maneuvering using tensor-decomposition-based globally optimal control and onboard visual-inertial estimation. Master’s thesis, Massachusetts Institute of Technology, 2016. [42] Mohammed Osman, Yuanqing Xia, Mohammed Mahdi, Tayyab Manzoor, Abdulrahman H. Bajodah, Asif Ali, Abid Ali, and Azzam Ahmed. An Adaptive SINDy-Lyapunov Model Predictive Control Framework for Dual-System VTOL UAVs. International Journal of Robust and Nonlinear Control, 2025. [43] Zeyad M. Manaa, Mohammed R. Elbalshy, and Ayman M. Abdallah. Data-driven Discovery of The Quadrotor Equations of Motion Via Sparse Identification of Nonlinear Dynamics. In AIAA SCITECH 2024 Forum, Orlando, FL, 2024. [44] Joseph B. Keller. Periodic oscillations in a model of thermal convection. Journal of Fluid Mechanics, 1966. [45] H. F. Creveling, J. F. De Paz, J. Y. Baladi, and R. J. Schoenhals. Stability characteristics of a single-phase free convection loop. Journal of Fluid Mechanics, 1975. [46] James A. Yorke, Ellen D. Yorke, and John Mallet-Paret. Lorenz-like chaos in a partial differential equation for a heated fluid loop. Physica D: Nonlinear Phenomena, 1987. [47] Jinzi Mac Huang and Nicholas J. Moore. A convective fluid pendulum revealing states of order and chaos. arXiv preprint arXiv:2307.13146, 2023. [48] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, and Sartaj Singh. SymPy: symbolic computing in Python. PeerJ Computer Science, 2017. [49] Rohit Supekar, Boya Song, Alasdair Hastewell, Gary PT Choi, Alexander Mietke, and Jörn Dunkel. Learning hydrodynamic equations for active matter from particle simulations and experiments. Proceedings of the National Academy of Sciences, 2023. [50] Robert K Niven, Ali Mohammad-Djafari, Laurent Cordier, Markus Abel, and Markus Quade. Bayesian identification of dynamical systems. Multidisciplinary Digital Publishing Institute Proceedings, 2020. [51] Seth M. Hirsh, David A. Barajas-Solano, and J. Nathan Kutz. Sparsifying priors for Bayesian uncertainty quantification in model discovery. Royal Society Open Science, 2022. [52] Lloyd Fung, Urban Fasel, and Matthew Juniper. Rapid Bayesian identification of sparse nonlinear dynamics from scarce and noisy data. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2025. [53] Lorenzo Boninsegna, Feliks Nüske, and Cecilia Clementi. Sparse learning of stochastic dynamical equations. The Journal of Chemical Physics, 2018.

15