In racdale/sindyr: Sparse Identification of Nonlinear Dynamics

knitr::opts_chunk$set(echo = TRUE)

From:

Dale, R. & Bhat, H. S. (in press). Equations of mind: Data science for inferring nonlinear dynamics of socio-cognitive systems. Cognitive Systems Research.

This brief code obtains a simulation of the Lorenz system using the resources in the crqa library (Coco and Dale, 2014). For illustration, we plot the system, variables, and the features that are used in SINDy. The code then runs the sindyr function to recover the system of equations.

library(sindyr)
library(crqa) # for Lorenz

set.seed(666)
dt = .001
numsteps = 50000; dt = dt; sigma = 10; r = 28; b = 2.6;
xs = data.frame(lorenzattractor(numsteps, dt, sigma, r, b, plots=F))
colnames(xs) = list('x','y','z')

xs = xs[2000:nrow(xs),] # cut out initialization

points3D(xs$x,xs$y,xs$z,type='l',col='black')

# plot the individual variables
par(mfrow=c(3,1),oma = c(0,0,0,0) + 0.1,mar = c(4.1,4.1,1,1) + 0.1)
plot(xs$x,type='l',xlab='t',ylab='x')
plot(xs$y,type='l',xlab='t',ylab='y')
plot(xs$z,type='l',xlab='t',ylab='z')

# plot the features
Theta = features(xs,3) # grid of features
par(mfrow=c(7,3),oma = c(2,0,0,0) + 0.1,mar = c(1,1,1,1) + 0.1)
for (i in 2:ncol(Theta)) {
  plot(Theta[,i],xlab='t',main=gsub(':','',colnames(Theta)[i]),type='l',xaxt='n',yaxt='n')
}

Here is the reconstruction of this system.

sindy.obj = sindy(xs=xs,dt=dt,lambda=.5) # let's reconstruct
sindy.obj$B

Let's initialize the parameters, the first system state, then simulated using our reconstructed B.

#
# let's model under particular values
#
par(mfrow=c(1,2))
x.now = xs[1,]
xs_modeled = xs[1,]
for (i in 1:10000) {
  x.now = x.now + (10^(-3))*features(x.now,3) %*% sindy.obj$B
  xs_modeled = rbind(xs_modeled,x.now)
}

points3D(xs_modeled$x,xs_modeled$y,xs_modeled$z,type='l',col='black')