Can parameterized_pattern_formation work with time?

[ziyinyuan]

I am trying a spring mass equation with x''+w^2*x=0 (w=omega)

My goal is to have the SINDy model would print something like (x)' = -omega sin(omega t). I can get this right now with one set of data with one omega value. I am wondering if I can generate multiple data with multiple omega value, then use the data to get the model using sindy_parametrized_pattern?

Here's my code for now:

def mass_spring_equation(y, t, w):
    x, dx_dt = y
    dy_dt = [dx_dt, -omega**2 * x]
    return dy_dt

# Parameters
m = 1.0  # Mass (kg)
k = 9.0  # Spring constant (N/m)
omega = np.sqrt(k / m)  # Angular frequency (rad/s)

# Initial conditions
x0 = 1.0  # Initial displacement (meters)
v0 = 0.0  # Initial velocity (m/s)

# Time points
t = np.linspace(0, 10, 1000)  # 0 to 10 seconds, 1000 time points for simulation

# Solve the ODE using scipy's odeint function
sol = odeint(mass_spring_equation, [x0, v0], t, args=(omega,))

# Extract displacement and velocity from the solution
displacement = sol[:, 0]
velocity = sol[:, 1]

t_2d=t.reshape(-1,1) # change row to column
omega_2d=omega*np.ones_like(t_2d)
u1=np.concatenate((omega_2d,t_2d),axis=1)

f1 = lambda omega, t: np.cos(omega * t)
f1_name = lambda omega, t: f"(cos({omega}*{t}))"
f2 = lambda omega, t: np.sin(omega * t)
f2_name = lambda omega, t: f"(sin({omega}*{t}))"
f3 = lambda omega, t: omega * np.cos(omega * t)
f3_name = lambda omega, t: f"{omega}*(cos({omega}*{t}))"
f4 = lambda omega, t: omega * np.sin(omega * t)
f4_name = lambda omega, t: f"{omega}*(sin({omega}*{t}))"

lib=ps.CustomLibrary(
    library_functions=[f1,f2,f3,f4],
    function_names=[f1_name,f2_name,f3_name,f4_name]
)
generalized_library=ps.GeneralizedLibrary(
    libraries=[ps.PolynomialLibrary(degree=1),lib],
    inputs_per_library=np.array([[0,0],[1,2]])
)

model = ps.SINDy(
    feature_library=generalized_library,
    optimizer=ps.STLSQ(threshold=0.5),
    feature_names=["x","omega","t"],
)
model.fit(displacement, t=t,u=u1)
model.print()

sim=model.simulate([displacement[0]],t,u=u1)

[Jacob-Stevens-Haas]

I'm less familiar with sindy_parametrized_pattern or the example notebook on parametrized pattern formation - that's @znicolaou 's bailiwick.

However, there's another issue, and that's that the two ODEs you wrote aren't equivalent. The standard way to do $\ddot x =- \omega x$ is to not split up sol into displacement and velocity, but to run model.fit(sol, t=t) without control variables.

[znicolaou]

Sure, should be possible to fit multiple trajectories for the first order model x'=p and p'=-w^2 x with differing w using the ParameterizedLibrary. I think I can write up a short example quickly in the next day or two.

[ziyinyuan]

@Jacob-Stevens-Haas Thank you for pointing this out. I'll definitely try this approach.

@znicolaou That would be great! Then I assume it also works with couple spring mass system?

[ziyinyuan]

I tried model.fit with 1 alpha value and it worked. So I am trying a simple multiple trajectory model. But I can't output the right model. Am I doing anything wrong here?

def mass_spring_equation(y, t, w):
    x, dx_dt = y
    dy_dt = [dx_dt, -omega**2 * x]
    return dy_dt

omega_values = [1,2,3]
u_trains = [[omega_values[i]] for i in range(len(omega_values))]

x0 = 1.0 
v0 = 0.0
t = np.linspace(0, 5, 1000)

results = []
for omega in omega_values:
    sol = odeint(mass_spring_equation, [x0, v0], t, args=(omega,))
    results.append(sol)

model = ps.SINDy(
    feature_library=ps.PolynomialLibrary(),
    optimizer=ps.STLSQ(threshold=0.1,normalize_columns=False),
    feature_names=["x","v","alpha"],
)
model.fit(results, u=u_trains, t=t, multiple_trajectories=True)
model.print()

[znicolaou]

Hey @ziyinyuan , you were almost there! You just need to tell ps.SINDy to use ps.ParameterizedLibrary for the feature_library option, where you can us instances of ps.PolynomialLibrary for the parameterized libraries feature_library and the parameter_library options, as in

def mass_spring_equation(y, t, w):
    x, dx_dt = y
    dy_dt = [dx_dt, -omega**2 * x]
    return dy_dt

omega_values = [1,2,3]
u_trains = [[omega_values[i]] for i in range(len(omega_values))]

x0 = 1.0 
v0 = 0.0
t = np.linspace(0, 5, 1000)

results = []
for omega in omega_values:
    sol = odeint(mass_spring_equation, [x0, v0], t, args=(omega,))
    results.append(sol)

lib=ps.ParameterizedLibrary(
    feature_library=ps.PolynomialLibrary(),
    parameter_library=ps.PolynomialLibrary(),
    num_features=2,
    num_parameters=1,
)
model = ps.SINDy(
    feature_library=lib,
    optimizer=ps.STLSQ(threshold=0.1,normalize_columns=False),
    feature_names=["x","v","alpha"],
)
model.fit(results, u=u_trains, t=t, multiple_trajectories=True)
model.print()

The num_features=2 corresponds to the fact that each trajectory has features x and v, and the num_parameters=1 corresponds to the single parameter alpha. That gives the right model for me:

(x)' = 1.000 1 v
(v)' = -1.000 alpha^2 x

[ziyinyuan]

I couldn't get this to work. Can someone check if I made any mistakes here?

def coupled_equations(x, t, zeta, w, beta):
    x1, x1_dot, x2, x2_dot = x

    dx1_dot_dt = -2 * zeta * x1_dot - w**2 * x1 + beta * w**2 * x2
    dx2_dot_dt = -2 * zeta * x2_dot - w**2 * x2 + beta * w**2 * x1

    return [x1_dot, dx1_dot_dt, x2_dot, dx2_dot_dt]

# Parameters
zeta_values = [0.1, 0.3, 0.5, 0.7, 0.9]
w_values = [1, 2, 3, 4, 5]
beta = 0
t = np.linspace(0, 60, 1000)  # Time values

# Generate the result list of arrays
result_list = []

# Loop through parameter combinations
for zeta_val, w_val in zip(zeta_values, w_values):
    x_initial = [0.5, 0.0, -1.0, 0.0] # Initial conditions for positions and velocities
    x_solution = odeint(coupled_equations, x_initial, t, args=(zeta_val, w_val, beta))
    result_list.append(np.array(x_solution))

# Convert the result list to an array
result_array = np.array(result_list)

u_trains = [[zeta_values[i],w_values[i]] for i in range(len(w_values))]
lib = ps.ParameterizedLibrary(
    feature_library=ps.PolynomialLibrary(degree=1),
    parameter_library=ps.PolynomialLibrary(degree=2),
    num_features=4, # each trajectories have x1,x1dot,x2,x2dot
    num_parameters=2, # 2 parameters (zeta and w)
)

model = ps.SINDy(
    feature_library=lib,
    optimizer=ps.STLSQ(threshold=0.2, normalize_columns=False),
    feature_names=["x1", "x1dot", "x2", "x2dot", "zeta","w"],
)
model.fit(result_list, u=u_trains, t=t, multiple_trajectories=True)
model.print()