To demonstrate how SINDy works in practice, we begin with a classic benchmark from the study of chaos: the Lorenz system. Originally derived as a simplified model for atmospheric convection, the Lorenz equations are now a textbook example of deterministic chaos in continuous dynamical systems.

The Lorenz system consists of three coupled nonlinear differential equations:

Using standard parameters \( \sigma = 10 \), \( \rho = 28 \), and \( \beta = \frac{8}{3} \), the system exhibits chaotic dynamics: sensitive dependence on initial conditions, bounded trajectories, and a fractal-like attractor in phase space. This attractor forms a twisted double-lobed structure where solutions appear to jump unpredictably between lobes — despite being fully deterministic.

We numerically integrate the Lorenz system using SciPy's solve_ivp, starting from an initial condition \( (-8, 7, 27) \) and sampling over 25 seconds with 10,000 time points. The resulting trajectory is shown below:

With this dataset in hand, we now apply SINDy to learn the underlying equations directly from the simulated data. Thanks to the PySINDy library, the process is remarkably straightforward:

The learned equations are:

        (x0)' = -10.003 x0 + 10.003 x1
        (x1)' = 27.804 x0 + -0.953 x1 + -0.994 x0 x2
        (x2)' = -2.667 x2 + 0.999 x0 x1
            

These results closely match the original Lorenz parameters. Despite being learned purely from data, the identified model captures both the structure and coefficients of the true system to high precision.

To assess the model’s predictive performance, we simulate both the original system and the learned SINDy model from the same initial condition. The trajectories for \( x(t) \) are plotted below:

To understand how SINDy works internally, let us step through the main components of the algorithm using the Lorenz system. This helps demystify the magic behind model.fit(X, t) and shows how the method discovers equations from data.

We start with our data — a matrix \( X \in \mathbb{R}^{m \times n} \) containing \( m \) time samples of the system’s state \( \mathbf{x}(t) \in \mathbb{R}^n \). But SINDy doesn’t learn from position alone — it needs the rate of change, the time derivatives \( \dot{\mathbf{x}}(t) \).

Since we don’t have analytical derivatives, we compute them numerically. The simplest method is finite differences: for small time steps \( \Delta t \),

This gives us \( \dot{X} \), a matrix of the same shape as \( X \). Libraries like pysindy provide smoother estimators, but the principle is the same.

The next step is to build a function library \( \Theta(X) \). This is a transformation of our data into many nonlinear combinations — the candidate right-hand-side terms in the differential equations. We’ll use all polynomial combinations up to degree 3.

As a small example: suppose \( x = 2 \) and \( y = -1 \). Then the second-degree polynomial library might include:

The full library \( \Theta(X) \in \mathbb{R}^{m \times p} \) evaluates such terms for all time points. In practice, this step is handled by PySINDy, but we can inspect its shape and content directly.

Now comes the key part: solving the sparse regression problem \( \dot{X} = \Theta(X) \cdot \Xi \). The matrix \( \Xi \) contains the coefficients of each function in the library. Most entries should be zero — only a few terms are active in each equation.

We use the STLSQ algorithm (Sequentially Thresholded Least Squares), which starts by fitting a standard least-squares solution, then iteratively sets small coefficients to zero and refits only the remaining ones. Here is a simplified version:

The threshold parameter controls sparsity — if set too high, important terms may be removed. If too low, the model becomes dense and overfits. With threshold = 0.1, the recovered coefficient matrix looks like:

        [[  0.           0.           0.        ]
        [-10.003        27.804        0.        ]
        [ 10.003        -0.953        0.        ]
        [  0.            0.          -2.667      ]
        [  0.            0.           0.999      ]
        [  0.           -0.994        0.        ]]
            

This matrix corresponds to the learned equations:

These match the original Lorenz system, and we obtained them using our own code. For convenience, we still use PySINDy’s feature library and differentiation tools, but the core regression is fully transparent.

A key question in data-driven modeling of chaotic systems is: how much data is needed to recover the correct dynamics? To answer this, we compare how well the Sparse Identification of Nonlinear Dynamical Systems (SINDy) framework performs on four different chaotic systems under varying amounts of available data. These systems differ in complexity and chaotic behavior, measured by their largest Lyapunov exponent.

We generate clean simulation data for each system and then subsample the time series to simulate lower data availability. Subsampling here refers to reducing the number of timepoints used for fitting, for instance by taking every k-th point or slicing the first N values from the total trajectory. We define a successful identification as finding exactly the correct terms in the governing equations — no more, no less.

Real-world measurements are rarely noise-free. Even the most precise experimental data contain fluctuations due to sensor error, digitization, or environmental interference. To test the robustness of sparse identification, we revisit our four chaotic systems and now add Gaussian noise to the simulated trajectories.

Noise is added to the state vector \( \mathbf{x}(t) \) such that: \[ \mathbf{x}_{\text{noisy}}(t) = \mathbf{x}(t) + \mathcal{N}(0, \sigma^2) \] where \( \sigma \) is a percentage of the standard deviation of the signal (here set to 1% of the signal amplitude). We then perform identification exactly as before using pysindy with polynomial libraries.

Due to the destabilizing nature of noise — especially in numerical differentiation — even small levels of perturbation drastically degrade performance. This test therefore acts as a stress test for both the numerical methods and the sparsity assumption in SINDy.