dynamicslab / pysindy

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

Some problems about the PDE-FIND accuracy. #563

Open Camelllia-L opened 1 week ago

Camelllia-L commented 1 week ago

Well, thanks a lot for your work. It's fantastic.

I'm now trying to identify the 1D Terzaghi Consolidation Equations (u_t = a*u_xx), but I met some problems. The accuracy is not very good, the PDE can only be identified with clean data or data with about 0.15% noise, and only with the weak form method(even clean data).

I don't know whether it is the problem of the data I input into the model. The data was obtained through Finite Difference Method.

Actually, I noticed your paper about PDE-FIND first, there is a poly- difference method in the STRidge algorithm. Can I use this difference method in pySINDy? Or in which ways can I obtain the same effect?

I tried the Smoothed-FD here but it has no effect, just like the normal FD. Does this problem mean that it has nothing to do with the differential method? If yes, then maybe it is the problem of input data itself?

Thank you!

Jacob-Stevens-Haas commented 1 week ago

Hey, can you start with the version of pysindy, and at the very least, the code you used to create the ps.SINDy object?

Camelllia-L commented 1 week ago

The version is 1.7.5. And the code is as follows,

Import

import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from IPython.display import Image import numpy as np import pysindy as ps from scipy.io import loadmat from sklearn.metrics import mean_squared_error

np.random.seed(100)

Data Prep

data = loadmat('data/VacuumFDM.mat') u = data['usol'] x = data['x'][:,0] t = data['t'][:,0] dt = t[1]-t[0] dx = x[2]-x[1]

u_dot = ps.FiniteDifference(axis=1)._differentiate(u, t=dt)

u = u.reshape(len(x), len(t), 1) u_dot = u_dot.reshape(len(x), len(t), 1)

Library Define

library_functions = [lambda x: x, lambda x: x x, lambda x: x x * x] library_function_names = [lambda x: x, lambda x: x + x, lambda x: x + x + x]

pde_lib = ps.PDELibrary( library_functions=library_functions, function_names=library_function_names, derivative_order=3, spatial_grid=x, diff_kwargs={"is_uniform": True, "periodic": True}, )

Fit

optimizer = ps.STLSQ(threshold=1, alpha = 1e-5, normalize_columns=False)
model = ps.SINDy(feature_library=pde_lib, feature_names=['u'], optimizer=optimizer) model.fit(u, t=dt) model.print()

Weak form SINDy

u = np.real(data['usol']) x = np.real(data['x'][:,0]) t = np.real(data['t'][:,0]) dt = t[1]-t[0] dx = x[2]-x[1]

u_dot = ps.SmoothedFiniteDifference(axis=1)._differentiate(u, t=dt) u = u.reshape(len(x), len(t), 1) u_dot = u_dot.reshape(len(x), len(t), 1)

library_functions = [lambda x: x, lambda x: x x, lambda x: x x * x] library_function_names = [lambda x: x, lambda x: x + x, lambda x: x + x + x]

X, T = np.meshgrid(x, t) XT = np.asarray([X, T]).T

pde_lib = ps.WeakPDELibrary( library_functions=library_functions, function_names=library_function_names, derivative_order=3, spatiotemporal_grid=XT, is_uniform=True, K=1000, include_bias=False )

optimizer = ps.STLSQ(threshold=1e-3, alpha = 1e-5, normalize_columns=True)
model = ps.SINDy(feature_library=pde_lib, feature_names=['u'], optimizer=optimizer) model.fit(u, t=dt) model.print()

Camelllia-L commented 1 week ago

I tried STRidge algorithm (from Data-driven discovery of partial differential equations) with the same dataset, the PDE can be identified with clean data in a strong form (grid size is 101*201, the same as the default size in pySINDy model).

However the parameter input into the may be wierd? w = TrainSTRidge(R,Ut,0.1,3000), means the dtol is 3000. Whatever, can this result rule out issues with the quality of the dataset? And how can I achieve the same effect in pySINDy? I noticed that there seems no iteration about the threshold (TrainSTRidge function in STRidge) in pySINDy.

makeabhishek commented 3 days ago

Can you please provide the link to the data?

Camelllia-L commented 3 days ago

@makeabhishek I upload the data in my repositories, the link is https://github.com/Camelllia-L/1DConsolidationVacuumAnalyticalData Thank you! A problem I found with this PDE is that tit will be discontinuous near point (0,0). The initial boundary is u(0,x) = 0, here 0<x<H. And the top boundary is u(t,0) = 1, here 0<t<T. I suspect this is what caused pysindy to fail.