Sindy Identification

SINDy — Sparse Identification of Nonlinear Dynamics

Overview

Discover governing ODEs from time-series data using SINDy (Brunton et al., PNAS 2016). Given measurements of state variables x(t), SINDy identifies the sparse set of terms in a library of candidate functions that best describe dx/dt.

When to Use

• You have trajectory data and want to find the governing ODE
• System is expected to have a sparse representation (few active terms)
• Input variables are known (you measured the right state variables)
• Data is relatively clean (or you can denoise it)

Do NOT Use When

• You want arbitrary symbolic expressions (use symbolic-regression with PySR)
• You only have steady-state data (SINDy needs time derivatives)
• System is stochastic (SINDy assumes deterministic dynamics)
• You have PDE data (use PDE-FIND variant, not standard SINDy)

Installation

pip install pysindy

Core Workflows

1. Basic SINDy (Discover ODE from Data)

import numpy as np
import pysindy as ps
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Generate training data from Lorenz system (ground truth)
def lorenz(t, y, sigma=10, rho=28, beta=8/3):
    return [sigma*(y[1]-y[0]), y[0]*(rho-y[2])-y[1], y[0]*y[1]-beta*y[2]]

dt = 0.001
t_train = np.arange(0, 10, dt)
sol = solve_ivp(lorenz, (0, 10), [1, 1, 1], t_eval=t_train, rtol=1e-10)
x_train = sol.y.T  # shape (N, 3)

# Fit SINDy model
model = ps.SINDy(
    feature_names=["x", "y", "z"],
    optimizer=ps.STLSQ(threshold=0.1),  # sparsity threshold
    feature_library=ps.PolynomialLibrary(degree=2),
)
model.fit(x_train, t=dt)
model.print()

# Expected output:
# x' = -10.000 x + 10.000 y
# y' = 28.000 x + -1.000 y + -1.000 x z
# z' = -2.667 z + 1.000 x y

2. Choosing the Sparsity Threshold

# Sweep thresholds and check model complexity vs error
thresholds = [0.001, 0.01, 0.05, 0.1, 0.5, 1.0]
for thresh in thresholds:
    model_t = ps.SINDy(
        optimizer=ps.STLSQ(threshold=thresh),
        feature_library=ps.PolynomialLibrary(degree=2),
    )
    model_t.fit(x_train, t=dt)
    complexity = model_t.complexity  # number of nonzero terms
    # Simulate and compute error
    x_sim = model_t.simulate(x_train[0], t_train)
    rmse = np.sqrt(np.mean((x_sim - x_train)**2))
    print(f"threshold={thresh:.3f}: {complexity} terms, RMSE={rmse:.4f}")

3. Custom Function Libraries

# Include trigonometric functions for oscillatory systems
library = ps.PolynomialLibrary(degree=2) + ps.FourierLibrary(n_frequencies=3)

# Or build a custom library
import pysindy as ps
custom_library = ps.CustomLibrary(
    library_functions=[
        lambda x: x,           # linear
        lambda x: x**2,        # quadratic
        lambda x: np.sin(x),   # sinusoidal
        lambda x: np.cos(x),
    ],
    function_names=[
        lambda x: x,
        lambda x: f"{x}^2",
        lambda x: f"sin({x})",
        lambda x: f"cos({x})",
    ]
)

model = ps.SINDy(
    feature_library=custom_library,
    optimizer=ps.STLSQ(threshold=0.1),
)
model.fit(x_train, t=dt)
model.print()

4. Noisy Data (Smoothed Differentiation)

# Add noise
x_noisy = x_train + 0.1 * np.random.randn(*x_train.shape)

# Use smoothed finite differences for derivatives
model_noisy = ps.SINDy(
    differentiation_method=ps.SmoothedFiniteDifference(),
    optimizer=ps.STLSQ(threshold=0.2),  # higher threshold for noisy data
    feature_library=ps.PolynomialLibrary(degree=2),
)
model_noisy.fit(x_noisy, t=dt)
model_noisy.print()

5. Validate: Simulate and Compare

# Simulate discovered model from same IC
x_sim = model.simulate(x_train[0], t_train)

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
labels = ['x', 'y', 'z']
for i, ax in enumerate(axes):
    ax.plot(t_train, x_train[:, i], 'b-', linewidth=0.5, label='Ground truth')
    ax.plot(t_train, x_sim[:, i], 'r--', linewidth=0.5, label='SINDy model')
    ax.set_ylabel(labels[i], fontsize=13)
    ax.legend(loc='upper right')
    ax.grid(True, alpha=0.3)
axes[-1].set_xlabel('Time [s]')
axes[0].set_title('SINDy Model vs Ground Truth')
plt.tight_layout()
plt.savefig('sindy_validation.png', dpi=150, bbox_inches='tight')

# Quantitative error
rmse = np.sqrt(np.mean((x_sim - x_train)**2, axis=0))
print(f"RMSE per variable: x={rmse[0]:.4f}, y={rmse[1]:.4f}, z={rmse[2]:.4f}")

6. Coefficient Extraction

# Get coefficient matrix (rows = equations, cols = library functions)
coefficients = model.coefficients()
feature_names = model.get_feature_names()

print("Discovered equations:")
for i, eq_name in enumerate(['dx/dt', 'dy/dt', 'dz/dt']):
    terms = []
    for j, name in enumerate(feature_names):
        if abs(coefficients[i, j]) > 1e-10:
            terms.append(f"{coefficients[i,j]:+.3f}·{name}")
    print(f"  {eq_name} = {' '.join(terms)}")

Key Parameters

Tips

1. Start with polynomial degree 2 — most physics ODEs are low-order polynomial
2. Threshold tuning is critical — too low = overfitting, too high = missing terms
3. Clean data first — SINDy is sensitive to noise in derivatives
4. Check the library — if the true terms aren't in the library, SINDy can't find them
5. Validate by simulation — a good SINDy model should reproduce the trajectory

Troubleshooting