search

racdale / sindyr

sindyr: Tools for Sparse Identification of Nonlinear Dynamics with R
14 stars 1 forks source link

Finite differences problem #1

Closed barneyricca closed 5 months ago

barneyricca commented 1 year ago

I created some test data to use with sindy(): set.seed(42) matrix(c(1:6 + rnorm(6, 0, 0.1), 10:15 + rnorm(6, 0, 0.1)), nrow = 6, ncol = 2, byrow = FALSE) -> z c("X1", "X2") -> colnames(z)

The command: sindy(x = z, verbose = TRUE) returns the following for dx: X1 X2 [1,] 0.950 1.025 [2,] 1.060 1.002 [3,] 1.002 0.964 [4,] 0.000 0.000 [5,] 0.000 0.000 [6,] 0.949 1.098

I don't think the fourth and fifth rows of dx[,] can be correct.

The resulting computed coefficients in matrix B are all over the place, too. I haven't (yet) dug into the finite_difference() code to find the mistake, but it seems to me that there must be one.

The issue appears to be in finite_difference(). Lines 14 & 15: for (i in 3:(n-1)) { fdx[i-2] = (x[i] - x[i-2]) / (2*S) # the intermediate ones

should be: for (i in 3:n) { fdx[i-1] = (x[i] - x[i-2]) / (2*S) # the intermediate ones

I had a colleague run SINDy using the Python package. The results report something slightly different for the 1st data point in each series, but supports my proposed change.

Barney

racdale commented 1 year ago

Hi there. Thanks for checking this and suggesting the change!

I'm convinced by your suggested change and have now committed the update to the repo. I am going to correspond with co-author to consider a note on the main publication so that readers can be aware. This does not impact the illustrations in the main paper, however it does potentially impact how it will be used in the future (esp. on deploying sample code for the paper and other applications). The shift yields r = 1 by a single time sample, but this may be non-trivial for stochastic processes with low autocorrelation.

Here's the update.

Thank you again!

racdale commented 5 months ago

Thanks again for this input. Function updated and some further minor revisions to the repo. I've now updated CRAN with version 0.2.4 and can close the thread.