Incorrect identified system equations for a coupled system with symmetric coefficients

[wanikhawar]

How can I resolve the issue of SINDy identifying incorrect system equations for the following system:

$\ddot{x_1}-\lambda_1 (1-x_1^2) \dot{x_1} + \alpha_1 x_1 = \beta_1 (x_1-x_2)$
$\ddot{x_2}-\lambda_2 (1-x_2^2) \dot{x_2} + \alpha_2 x_2 = \beta_2 (x_2-x_1)$

Here's my code

# import all the necessary libraries
import pysindy as ps
import numpy as np
from scipy.integrate import solve_ivp

# The coupled van der pol oscillator system 
def coupled_vdp(t, x, lambda1=3, lambda2=3, alpha1=1.5, alpha=1.5, beta1=0.5, beta2=0.5):
    return [x[1],
           beta1*(x[1]-x[3]) + lambda1*(1-x[0]**2) * x[1] - alpha1**2 * x[0],
           x[3],
           beta1*(x[3]-x[1]) + lambda1*(1-x[2]**2) * x[3] - alpha2**2 * x[2]]

# generate the data for testing
start = 0
end = 100
step = 0.01
t_span = [start, end]
x0 = [0.1, 0, 0.1, 0]
sol = solve_ivp(coupled_vdp, t_span, x0, dense_output=True)
t = np.arange(start, end+step/10, step)
x = sol.sol(t).T

# The feature names
feature_names = ['x1', 'x1dot', 'x2', 'x2dot']

feature_library = ps.PolynomialLibrary(degree=3)

# Differentiation method
differentiation_method = ps.FiniteDifference(order=2)


# Optimizer
opt = ps.STLSQ(threshold=0.1)
# opt = ps.SR3(threshold=0.1,thresholder="L0")
# opt = ps.SSR()

# Building the model
model=ps.SINDy(differentiation_method=differentiation_method,
               feature_library = feature_library,
               feature_names=feature_names, 
              optimizer=opt)


# Fitting the model to the training data
model.fit(x,
          t=t)

# Print the model
model.print()
# model.get_feature_names()

If I take $\lambda_1 = \lambda_2$, $\alpha_1 = \alpha_2$ and $\beta_1 = \beta_2$, SINDy identifies incorrect system equations.
For example for $\lambda_1 = \lambda_2 = 3$, $\alpha_1 = \alpha_2 = 1.5$ and $\beta_1 = \beta_2 = 0.5$, SINDy identified:

(x1)' = 0.500 x1dot + 0.500 x2dot

(x1dot)' = -1.126 x1 + 1.500 x1dot + -1.126 x2 + 1.500 x2dot + -0.500 x1^2 x1dot + -0.500 x1^2 x2dot + -0.500 x1 x1dot x2 + -0.500 x1 x2 x2dot + -0.500 x1dot x2^2 + -0.500 x2^2 x2dot

(x2)' = 0.500 x1dot + 0.500 x2dot

(x2dot)' = -1.126 x1 + 1.500 x1dot + -1.126 x2 + 1.500 x2dot + -0.500 x1^2 x1dot + -0.500 x1^2 x2dot + -0.500 x1 x1dot x2 + -0.500 x1 x2 x2dot + -0.500 x1dot x2^2 + -0.500 x2^2 x2dot

I've tried using STLSQ, SR3 and SSR but all have the same issue.

How do I resolve this? Am I supposed to impose the constraints:
(x1)' = x1dot
and
(x2)' = x2dot

I tried doing this but I didn't understand the documentation.

Additionally the documentation has broken links that don't lead anywhere.

Most links on the following page are broken.

The page with broken links

[akaptano]

We are fixing the links, but I cannot debug your issues without seeing the code you are running

[wanikhawar]

@akaptano I have added my code to the question.

[akaptano]

There is nothing wrong with your code -- it is the physics you need to check. You have chosen an initial condition that guarantees x1(t) = x2(t) for all time (plot the trajectories and you will see) and therefore x1' = 0.5 * x1dot + 0.5 * x2dot is correct.

Try with the initial condition [1, 0, 3, 0] and your problem will fix itself.

Best,
Alan

[wanikhawar]

Thanks.
And sorry for wasting your time with this silly question.

[akaptano]

no worries, didn't take long to find since I have fiddled with these methods extensively. Best wishes.