Discovery of Nonlinear Dynamics with Automated Basis Function Generation

Mohammad Amin Basiri  ·  Charles Nicholson
Data Science and Analytics Institute & School of Industrial and Systems Engineering
University of Oklahoma, Norman, OK, USA

---

Overview

AutoSINDy is a hybrid Discovery-then-Solve framework that automatically identifies sparse governing equations of nonlinear dynamical systems from noisy observational data without any prior specification of the system's mathematical form.

Standard SINDy methods require a researcher-specified candidate library. If the right basis functions are omitted, the method fails. Symbolic regression is flexible but noise-sensitive and often produces bloated, unstable equations. AutoSINDy bridges both worlds:

Noisy Data  →  [PySR: Symbolic Mining]  →  [Curation Pipeline]  →  [Ensemble SINDy]  →  Sparse Governing Equations

Across 540 total experimental trials on six canonical nonlinear systems and six noise levels, AutoSINDy achieves:

Metric	AutoSINDy	Standard SINDy	Standard PySR
Excellent Derivative Prediction (R² ≥ 0.99)	92.8%	58.9%	28.3%
Excellent Simulation Stability (R² ≥ 0.99)	68.3%	18.9%	8.3%
Hard Simulation Crashes	0 / 180	20 / 180	14 / 180
Median Equation Complexity	≈ Ground Truth	10–29× over	Under-complex

---

The Three-Stage Pipeline

Stage 1: Discover (Symbolic Mining via PySR)

PySR is applied to K short, randomly sampled data windows ("chunks"). Every expression on the Pareto-optimal front is harvested for each chunk and each state variable, building a diverse pool of candidate functional forms.

Stage 2: Curate (Library Construction)

The raw expression pool is processed through three sequential operations:

1. Symbolic Decomposition: compound expressions are split into additive atoms; numeric prefactors are stripped.
2. Algebraic Expansion: atoms are optionally expanded (configurable: gentle, severe, or hybrid) and sorted by SymPy operator count to enforce a simplicity bias.
3. Collinearity Pruning: a greedy forward-selection procedure rejects any atom whose Pearson correlation with already-accepted terms exceeds a threshold ρ_max, implementing Occam's razor at the library level. VIF-based pruning is also supported.

Stage 3: Identify (Ensemble SINDy)

The curated library is handed to an ensemble SINDy optimizer (STLSQ or SR3) running over B bootstrap replicates. Terms are retained only if their inclusion probability across the ensemble exceeds a hard cutoff κ = 0.80, providing robust noise rejection.

---

Benchmark Systems

System	Governing Equations	True Complexity
Harmonic Oscillator	ẋ₀ = x₁, ẋ₁ = −k₁x₀ − k₂x₁	3
Damped Pendulum	ẋ₀ = x₁, ẋ₁ = −bx₁ − c·sin(x₀)	4
Modulated Oscillator	ẋ₀ = x₁, ẋ₁ = −bx₁·cos(x₀) − kx₀	5
Van der Pol	ẋ₀ = x₁, ẋ₁ = µ(1 − x₀²)x₁ − x₀	6
Duffing Oscillator	ẋ₀ = x₁, ẋ₁ = −δx₁ − αx₀ − βx₀³	6
Complex Lorenz	ẋ₀ = σ(x₁−x₀), ẋ₁ = x₀(ρ−x₂)−x₁, ẋ₂ = x₀x₁ − βx₂ + γx₁·sin(x₀+x₂)	15

Experiments sweep 6 noise levels (σ ∈ {0, 0.01, 0.02, 0.03, 0.04, 0.05}) × 5 random seeds per system = 180 trials per method, 540 total.

---

Repository Structure

AutoSINDy/
│
├── AutoSINDy.py              # Core framework: all three pipeline stages + evaluation
├── systems.py                # Data generators for all benchmark dynamical systems
├── run_sweep_v2.py           # Full publication sweep with crash-recovery logic
├── Plot_raw_signal.py        # Visualization: states, derivatives, phase portrait
├── Paper_Figures.ipynb       # Visualization: visualizing the summary of the results
├── requirements.txt          # Python dependencies
│
├── figures/                  # Output figures (generated by sweep)
├── tables/                   # Result tables
│
├── autosindy_results_log_tidy - sweep on noise - V...csv  # Noise sweep results
├── autosindy_started_log.json               # Crash-recovery sidecar log
└── table_results_summary.csv               # Aggregated summary table

---

Installation

Prerequisites: Python 3.9+

git clone https://github.com/mabasiri95/AutoSINDy.git
cd AutoSINDy
pip install -r requirements.txt

Dependencies

matplotlib==3.10.7
scipy==1.16.3
numpy==1.26.4
pandas==2.3.3
pysindy==1.7.5
pysr==1.5.8
sympy==1.14.0
seaborn==0.13.2

Note on PySR: PySR requires Julia to be installed and will attempt to install it automatically on first run. See the PySR documentation for details.

---

Quick Start

Run a Single Experiment

Open AutoSINDy.py and locate the config dictionary near the top of the file. Set your desired system and parameters, then run:

python AutoSINDy.py

Key configuration options:

config = {
    "system_to_run": "damped_pendulum",   # Which system to identify
    "use_unified_library": False,          # False = separate library per state var (recommended)
    "discovery_chunks": 10,                # Number of PySR bootstrap chunks (K)
    
    "data_params": {
        "damped_pendulum": {
            "noise_level": 0.05,           # Noise-to-signal ratio σ
            "noise_seed": 32,              # For reproducibility
        }
    },
    
    "pysr_params": {
        "niterations": 40,
        "random_state": 32,
        "deterministic": True,
        "parallelism": "serial",
    },
    
    "curation_params": {
        "pruning_method": "correlation",   # 'correlation' or 'vif'
        "correlation_threshold": 0.95,     # ρ_max
        "expansion_strategy": "gentle",    # 'gentle', 'severe', or 'hybrid'
    },
    
    "optimizer_params": {
        "name": "STLSQ",                   # 'STLSQ' or 'SR3'
        "threshold": 0.21,                 # Sparsity threshold λ
        "use_ensemble": True,              # Ensemble SINDy (E-SINDy)
        "n_models": 20,                    # Bootstrap replicates B
        "ensemble_aggregator": "median",
    },
}

Visualize Raw Data

python Plot_raw_signal.py

This generates a three-panel dashboard (state trajectories, derivatives, phase portrait) and saves it as publication_figure.svg.

To switch systems, edit the generator call at the top of the script. Examples for damped pendulum, Duffing, and complex Lorenz are all included.

Run the Full Publication Sweep

python run_sweep_v2.py

Set SWEEP_GROUP at the top of run_sweep_v2.py to control what runs:

Value	Description
'main'	All 6 systems × 6 noise levels × 5 seeds (180 trials)
'ablation'	Expansion strategy / library / optimizer ablations
'chunks'	Chunk-count sensitivity analysis
'all'	Everything above

Crash Recovery: The sweep automatically handles interruptions. A sidecar JSON log (autosindy_started_log.json) records each run's status before it starts. On restart, completed runs are skipped; crashed runs are retried with a fallback seed (nominal_seed + 100). Results accumulate in autosindy_results_log_tidy.csv.

---

Available Dynamical Systems (systems.py)

Each generator returns (X_noisy, X_dot, t_eval) and accepts a noise_seed parameter for reproducibility.

import systems
import numpy as np

# Damped Pendulum
X, Xdot, t = systems.generate_damped_pendulum_data(
    b=0.25, c=5.0, x0=[np.pi - 0.1, 0], t_end=10, n_samples=5000, noise_level=0.05
)

# Van der Pol
X, Xdot, t = systems.generate_vanderpol_data(mu=2.0)

# Lorenz
X, Xdot, t = systems.generate_lorenz_data(sigma=10.0, rho=28.0, beta=8./3.)

# Complex Lorenz (with cross-term γx₁sin(x₀+x₂))
X, Xdot, t = systems.generate_complex_lorenz_data(gamma=1.5)

# Duffing Oscillator
X, Xdot, t = systems.generate_duffing_data(delta=0.3, alpha=-1.0, beta=1.0)

# Harmonic Oscillator
X, Xdot, t = systems.generate_harmonic_oscillator_data(k1=5.0, k2=1.0)

# Modulated Oscillator
X, Xdot, t = systems.generate_modulated_oscillator_data(b=0.25, k=5.0)

---

Evaluation Metrics

AutoSINDy reports five complementary metrics, each probing a different aspect of discovery quality:

Metric	What It Tests
Derivative R² / MSE (train/test)	Local structural fit on training data
Derivative R² / MSE (new trajectory)	Generalization to unseen initial conditions
Simulation R² / MSE	Long-horizon stability under forward integration
Canonical Complexity	sympy.count_ops(sympy.expand(eq)) a fair cross-method comparison
Recovery Rate	Fraction of trials with new-trajectory R² > 0.99 (proxy for exact identification)

Discovery and simulation wall-clock times are also logged separately.

---

Key Results

Overall Reliability (180 trials per method, all systems & noise levels)

AutoSINDy     Derivative:  92.8% Excellent │  6.1% Good │  0.6% Poor │  1.1% Failed
              Simulation:  68.3% Excellent │  8.3% Good │  6.1% Poor │ 21.1% Failed

Std. SINDy    Derivative:  58.9% Excellent │ 11.1% Good │ 12.8% Poor │ 17.2% Failed
              Simulation:  18.9% Excellent │  8.3% Good │ 21.1% Poor │ 51.7% Failed

Std. PySR     Derivative:  28.3% Excellent │ 32.8% Good │ 30.0% Poor │  8.9% Failed
              Simulation:   8.3% Excellent │ 16.1% Good │ 75.0% Poor

Model Complexity (Median canonical operator count)

System	Ground Truth	AutoSINDy	Std. SINDy	Std. PySR
Harmonic Osc.	3	≈3	≈68	≈5
Damped Pendulum	4	≈5	≈59	≈6
Modulated Osc.	5	≈4	≈71	≈6
Van der Pol	6	≈12	≈175	≈8
Duffing Osc.	6	≈7	≈11	≈7
Complex Lorenz	15	≈17	≈35	≈9

Simulation Failure Counts (out of 30 trials per system)

System	AutoSINDy	Std. SINDy	Std. PySR
Harmonic Osc.	0	4	5
Damped Pendulum	0	2	2
Modulated Osc.	0	2	0
Van der Pol	0	6	0
Duffing Osc.	0	8	7
Complex Lorenz	0	0	0

---

How AutoSINDy Outperforms Its Components

Why not just PySR? Without a sparsity-promoting mechanism, genetic programming optimizes fit on training chunks; even small coefficient errors accumulate into trajectory divergence. PySR achieves 0% recovery on Van der Pol because it consistently finds the right form but not the right coefficients.

Why not just SINDy with an enriched library? The polynomial-Fourier basis introduces severe multicollinearity. Under noisy derivatives, the optimizer distributes coefficient mass across hundreds of correlated spurious terms, producing equations with 10–29× ground-truth complexity that diverge immediately under integration.

AutoSINDy resolves both failure modes by construction. The Discover stage mines system-specific functional forms (including non-standard compound terms like x₁·sin(x₀+x₂)). The Curate stage eliminates collinearity before any regression is performed. The Identify stage concentrates coefficient mass on the small number of genuinely active terms.

---

Citation

If you use AutoSINDy in your research, please cite:

@misc{basiri2026discoverynonlineardynamicsautomated,
      title={Discovery of Nonlinear Dynamics with Automated Basis Function Generation}, 
      author={Mohammad Amin Basiri and Charles Nicholson},
      year={2026},
      eprint={2605.09696},
      archivePrefix={arXiv},
      primaryClass={cs.LG},
      url={https://arxiv.org/abs/2605.09696}, 
}

---

Limitations and Future Work

• Computational cost scales linearly with the number of state variables and chunks. High-dimensional systems (n >> 3) may require amortized symbolic regression or parallel PySR calls.
• Chaotic long-horizon simulation: Even with the correct governing equations, chaotic sensitivity means any finite-precision model will eventually diverge from a specific reference trajectory. Attractor-geometry metrics (e.g., Lyapunov exponents, attractor dimension) are more appropriate for assessing long-horizon performance on chaotic systems.
• PDEs and real experimental data: All benchmarks use synthetic ODE data with controlled Gaussian noise. Extension to partial differential equations and real experimental time series (irregular sampling, partial observability) is an active direction.

---

License

This project is licensed under the MIT License. See the LICENSE file for details.

---

Contact

Mohammad Amin Basiri	ma.basiri@ou.edu
Charles Nicholson	cnicholson@ou.edu
GitHub	https://github.com/mabasiri95/AutoSINDy