Need help with PDE

[ziyinyuan]

I am trying apply SINDy on PDE using the kdv example you have. So I first start with loading the kdv file from matlab and just to plot the first set to see the result:

kdV = loadmat('data/kdv.mat')
t = np.ravel(kdV['t'])
x = np.ravel(kdV['x'])
u = np.real(kdV['usol'])
dt = t[1] - t[0]
dx = x[1] - x[0]
plt.figure()
for i in range(1):
  plt.plot(t,u[i,:])
plt.grid(True)

And then just follow along with the pde_lib, optimizer, and create the SINDy model:

u1 = u.reshape(len(x), len(t), 1)
library_functions = [lambda x: x,
                     lambda x: x * x]
library_function_names = [lambda x: x,
                          lambda x: x + x]
pde_lib = ps.PDELibrary(library_functions=library_functions,
                        function_names=library_function_names,
                        derivative_order=2, spatial_grid=x,
                        include_bias=True, is_uniform=True)

# Fit the model with different optimizers.
# Using normalize_columns = True to improve performance.
print('STLSQ model: ')
optimizer = ps.STLSQ(threshold=0.01, alpha=1e-5, normalize_columns=True)
model1 = ps.SINDy(feature_library=pde_lib, optimizer=optimizer)
model1.fit(u1, t=dt)
model1.print()

But the result is not matching when I plot it out:

sol1=model1.predict(u1)
plt.figure()
for i in range(1):
  plt.plot(t,u[i,:])
  plt.plot(t,sol1[i,:],'k--')
plt.grid(True)

Can someone help me with this? I'll be greatly appreciated.

[znicolaou]

Hello @ziyinyuan,

The issue is that the SINDy.predict function returns the right hand side of the SINDy model, which corresponds to the time derivatives of the input data rather than the input itself. So your sol1[i,:] should be compared to ps.FiniteDifference(d=1,axis=1)._differentiate(u,dt)[i,:]. Note also that the time derivatives for i=1 are all almost zero, so they are a little noisy. You'll see a better comparison if you use, say i=200, as in

sol1=model1.predict(u1)
plt.figure()
for i in range(200,201):
  plt.plot(t,ps.FiniteDifference(d=1,axis=1)._differentiate(u,dt)[i,:])
  plt.plot(t,sol1[i,:],'k--')
plt.grid(True)

[znicolaou]

PS: if you'd like to use the SINDy model to predict the values of u from an initial conditions rather than the time derivative, you'd need to numerically integrate the SINDy model. This is only implemented in PySINDy for ODE models (with the model.simulate function), since numerically integrating PDE models requires many specific choices, including boundary conditions and discretization methods. I'd recommend looking into either PDE integration with spectral or finite different methods and implementing something outside of PySINDy by retrieving the model coefficients with the model.optimizer.coef_ attribute, for example.