RasmahRasmah

Sparse Identification of Nonlinear Dynamics (SINDy)

Overview

SINDy (Brunton, Proctor & Kutz 2016) discovers governing equations directly from time-series data. It builds a library Θ(x) of candidate monomials and solves ẋ = Θ(x) Ξ for a sparse coefficient matrix Ξ with sequentially thresholded least squares — the small entries are zeroed, so only the true terms survive. The result is a human-readable equation, not a black box.

using Rasmah
using LinearAlgebra

Discover a damped oscillator

Generate snapshots of ẋ = y, ẏ = −x − 0.2y and let SINDy recover the sparse dynamics from a degree-2 polynomial library.

m = 400
X = zeros(2, m); X[:, 1] = [1.0, 0.0]
dt = 0.01
for k in 1:(m - 1)
    x, y = X[:, k]
    X[:, k + 1] = [x + dt * y, y + dt * (-x - 0.2y)]
end
Xdot = [X[2, :]'; (-X[1, :] .- 0.2 .* X[2, :])']   # ẋ = y, ẏ = −x − 0.2y

Ξ, exps = sindy(X, Xdot, 2)
println("monomials: ", exps)
monomials: [[0, 0], [0, 1], [1, 0], [0, 2], [1, 1], [2, 0]]

Verification — the coefficients match

The true dynamics ẋ = y, ẏ = −x − 0.2y correspond to the monomials x₂, x₁ and x₂ with coefficients 1, −1 and −0.2:

# monomial index of x₁ (=exps == [1,0]) and x₂ (=exps == [0,1])
i1 = findfirst(e -> e == [1, 0], exps)
i2 = findfirst(e -> e == [0, 1], exps)
println("ẋ = ", round(Ξ[1, i2]; digits = 3), "·x₂  (true 1)")
println("ẏ = ", round(Ξ[2, i1]; digits = 3), "·x₁ + ", round(Ξ[2, i2]; digits = 3), "·x₂  (true −1, −0.2)")
ẋ = 1.0·x₂  (true 1)
ẏ = -1.0·x₁ + -0.2·x₂  (true −1, −0.2)

Next steps