simulating error(overflow) when expanded to the third order

[BK201-kkk]

I am currently utilizing the SINDy (Sparse Identification of Nonlinear Dynamical Systems) method to capture the dynamics of fluid Proper Orthogonal Decomposition (POD) modes. While I have been successful in modeling the system up to the second order, I have encountered a simulation error when attempting to extend the model to the third order.

here is my code:
optimizer = ps.STLSQ(threshold=0.00001, alpha=0, max_iter=1000)
library = ps.PolynomialLibrary(degree=3, include_bias=False)
model = ps.SINDy(optimizer=optimizer, feature_library=library)
model.fit(normalized_coefs_i.T, t=dt2)
tsim =np.linspace(0,1000,15000)
x_sim2 = model.simulate(normalized_coefs_i.T[200, :], tsim)

I have a set of 8-dimensional data with 4000 time steps. I used interpolation method to make it into 8-dimensional data with 12000 time points (the time range is still 0-4000). A third-order SINDy model was trained. However, in the simulate step, it is difficult to train data with a time range exceeding 650. If you increase the simulate time, overflow will occur(this problem does not occur in the second-order model). I suspect that my wrong choice of step size is causing the numerical problem. Do you know the integration method used by psindy so that I can calculate the appropriate step size based on the relevant formula? Or do you have another opinion on this issue? Thank you so much!

This is the equation obtained from one of the modes:
(x0)' = 0.047702 x0 + 0.698551 x1 + -0.114620 x2 + 0.045224 x3 + 0.036848 x4 + 0.016167 x5 + -0.011996 x6 + 0.002499 x7 + -0.003969 x0^2 + -0.004245 x0 x1 + -0.002821 x0 x2 + -0.000790 x0 x3 + 0.006424 x0 x4 + 0.006871 x0 x5 + -0.006000 x0 x6 + 0.001644 x0 x7 + 0.014452 x1^2 + 0.012945 x1 x2 + 0.038963 x1 x3 + 0.002168 x1 x4 + 0.025544 x1 x5 + -0.003518 x1 x6 + 0.001677 x1 x7 + -0.018828 x2^2 + 0.015936 x2 x3 + -0.000574 x2 x4 + 0.023369 x2 x5 + 0.006994 x2 x6 + 0.004906 x2 x7 + 0.014055 x3^2 + -0.014039 x3 x4 + 0.024345 x3 x5 + -0.000743 x3 x6 + 0.010239 x3 x7 + -0.006994 x4^2 + -0.007396 x4 x5 + 0.002928 x4 x6 + 0.001220 x4 x7 + -0.001067 x5^2 + 0.000672 x5 x6 + -0.003117 x5 x7 + -0.005111 x6^2 + -0.002606 x6 x7 + -0.005483 x7^2 + -0.085164 x0^3 + -0.008811 x0^2 x1 + -0.048474 x0^2 x2 + -0.058962 x0^2 x3 + -0.056181 x0^2 x4 + -0.051650 x0^2 x5 + 0.021425 x0^2 x6 + -0.009091 x0^2 x7 + -0.077157 x0 x1^2 + 0.046257 x0 x1 x2 + -0.021555 x0 x1 x3 + -0.021816 x0 x1 x4 + 0.017303 x0 x1 x5 + 0.006793 x0 x1 x6 + 0.006328 x0 x1 x7 + -0.084383 x0 x2^2 + 0.015495 x0 x2 x3 + -0.030385 x0 x2 x4 + -0.089059 x0 x2 x5 + 0.036263 x0 x2 x6 + -0.015960 x0 x2 x7 + -0.027957 x0 x3^2 + -0.027038 x0 x3 x4 + -0.011440 x0 x3 x5 + -0.001476 x0 x3 x6 + 0.006846 x0 x3 x7 + -0.046979 x0 x4^2 + -0.027632 x0 x4 x5 + 0.011633 x0 x4 x6 + 0.013083 x0 x4 x7 + -0.030085 x0 x5^2 + 0.004959 x0 x5 x6 + 0.009029 x0 x5 x7 + -0.007982 x0 x6^2 + -0.002972 x0 x6 x7 + -0.007946 x0 x7^2 + -0.008613 x1^3 + -0.009264 x1^2 x2 + -0.035641 x1^2 x3 + -0.040932 x1^2 x4 + -0.038617 x1^2 x5 + 0.013465 x1^2 x6 + -0.003863 x1^2 x7 + -0.002351 x1 x2^2 + -0.006123 x1 x2 x3 + -0.062009 x1 x2 x4 + 0.023593 x1 x2 x5 + 0.007480 x1 x2 x6 + 0.018165 x1 x2 x7 + -0.016489 x1 x3^2 + -0.046186 x1 x3 x4 + -0.044545 x1 x3 x5 + 0.012919 x1 x3 x6 + -0.013692 x1 x3 x7 + -0.009534 x1 x4^2 + -0.010912 x1 x4 x6 + 0.007654 x1 x4 x7 + 0.002451 x1 x5^2 + -0.014082 x1 x5 x6 + -0.001499 x1 x5 x7 + 0.005488 x1 x6^2 + -0.001060 x1 x6 x7 + 0.002456 x1 x7^2 + -0.006246 x2^3 + -0.009150 x2^2 x3 + 0.048150 x2^2 x4 + -0.025777 x2^2 x5 + -0.006008 x2^2 x6 + -0.008126 x2^2 x7 + 0.014103 x2 x3^2 + -0.036111 x2 x3 x4 + -0.015379 x2 x3 x5 + 0.012985 x2 x3 x6 + 0.001658 x2 x3 x7 + 0.002437 x2 x4^2 + -0.003909 x2 x4 x5 + 0.005555 x2 x4 x6 + -0.037399 x2 x4 x7 + -0.008545 x2 x5^2 + -0.019145 x2 x5 x6 + -0.008405 x2 x5 x7 + 0.008877 x2 x6^2 + -0.000243 x2 x6 x7 + 0.008131 x2 x7^2 + -0.019641 x3^3 + -0.002111 x3^2 x4 + -0.030087 x3^2 x5 + -0.003555 x3^2 x6 + -0.002981 x3^2 x7 + 0.016271 x3 x4^2 + -0.015623 x3 x4 x5 + 0.008905 x3 x4 x6 + -0.005665 x3 x4 x7 + -0.006595 x3 x5^2 + 0.000828 x3 x5 x6 + 0.004818 x3 x5 x7 + -0.015619 x3 x6^2 + 0.000014 x3 x6 x7 + -0.012090 x3 x7^2 + 0.002953 x4^3 + 0.003208 x4^2 x5 + 0.001527 x4^2 x6 + -0.020351 x4^2 x7 + -0.002172 x4 x5^2 + -0.001665 x4 x5 x6 + -0.007673 x4 x5 x7 + 0.004046 x4 x6^2 + -0.000794 x4 x6 x7 + 0.007569 x4 x7^2 + 0.005030 x5^3 + -0.002880 x5^2 x6 + -0.011006 x5^2 x7 + 0.007198 x5 x6^2 + 0.002669 x5 x6 x7 + 0.008865 x5 x7^2 + 0.000707 x6^3 + -0.003524 x6^2 x7 + -0.001119 x6 x7^2 + -0.004221 x7^3

[Jacob-Stevens-Haas]

The TL;DR is that the code uses scipy.integrate.odeint or solve_ivp depending on integration keywords.

Also, to clarify, by "third order", you mean $d^3 x/dt^3=f(x, x_t, x_{tt})$? Or do you mean 3rd degree polynomial?

[BK201-kkk]

@Jacob-Stevens-Haas
Thanks for your answer, I was referring to 3rd degree polynomial
library = ps.PolynomialLibrary(degree=3, include_bias=False)

As rigidity increases, if the time step of the numerical method is not sufficient to capture the rapid changes in the system, it may lead to numerical instability(I'm wondering if overflow is due to this reason).
ValueError: Input contains NaN, infinity or a value too large for dtype('float64')
If you have any further suggestions, could you please point them out?

[Jacob-Stevens-Haas]

In general, there's no guarantee that models discovered by pysindy are bounded. 3rd degree polynomial ODEs can experience blow up, reaching inifinity in finite time. Even if blow up doesn't occur, yes, models can get too rigid for numerical solvers.

The TrappingSR3 and StabilizedLinearSR3 optimizers provide a way to discover "guaranteed" stable models, but not of 3rd degree polynomials. These require additional math.