dynamicslab / pysindy

A package for the sparse identification of nonlinear dynamical systems from data
https://pysindy.readthedocs.io/en/latest/
Other
1.36k stars 304 forks source link

polynomial library to include second order derivative #503

Open xiongxin9000 opened 1 month ago

xiongxin9000 commented 1 month ago

Hi,

I'm trying to reproduce a second order ode with polynomial library. I am struggling adding the second order derivative in the second equation. How can I achieve this by polynomial library. My code is as follows,

import pysindy as ps
import numpy as np
import matplotlib.pyplot as plt
model=ps.SINDy()
print(model)
########################
#pysindy.differentiation: X'
#ps.optimizers:coffecient
#ps.feature_library:basic function 
#########################
#import data from data.npy
#change to current directory
import os
os.chdir("xxxxxxxxxx")
data = np.loadtxt("xxxxxx")
data_lbm=data[:250,2]
R=1/data_lbm

R1 = R
R2 = model.differentiate(R1)
R3 = model.differentiate(R2)
t = np.arange(0, 250, 1)
# print(len(t))
x = np.column_stack((R1, R2))
# print(x.shape)
opt = ps.STLSQ(threshold=0, fit_intercept=False)
# fourier_library=ps.FourierLibrary(n_frequencies=2)
# model=ps.SINDy(feature_library=fourier_library,feature_names=["R1","R2"],optimizer=opt)

model=ps.SINDy(feature_names=["R1","R2","R3"],optimizer=opt,feature_library=ps.PolynomialLibrary(degree=3))
model.fit(x=x,t=t)
print('polinomial')
model.print()
R_model=model.simulate(x[0,:],t)
fig,ax= plt.subplots(1,2,figsize=(12,6))
ax[0].plot(R,label="lbm")
ax[0].plot(R_model[:,0],'--',label="SINdy")
ax[0].set_title("polinomial library")
ax[0].set(xlabel="t", ylabel="R")
ax[0].legend()

I want to include R3 into the second equation,

but what I got is,
polinomial (R1)' = 1.000 R2 (R2)' = -0.014 1 + 0.001 R1 + 0.293 R2 + -0.077 R1 R2 + 7.336 R2^2 + -0.004 R1^2 R2 + 1.531 R1 R2^2 + -163.617 R2^3

Thank you in advance, Best Regards, Xin

Jacob-Stevens-Haas commented 1 month ago

It is in the second equation, on the LHS. R2' is R3.

xiongxin9000 commented 1 month ago

Thank you for your reply Jacob I mean add R3 in RHS term because the terms on the right hand side only include R1 R2