Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

polynomial library to include second order derivative #503

Open
xiongxin9000 opened this issue May 14, 2024 · 2 comments
Open

polynomial library to include second order derivative #503

xiongxin9000 opened this issue May 14, 2024 · 2 comments

Comments

@xiongxin9000
Copy link

xiongxin9000 commented May 14, 2024

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
Copy link
Collaborator

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

@xiongxin9000
Copy link
Author

xiongxin9000 commented May 14, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants