Weak SINDY for higher order ODEs

[bthakur21]

Hi,

Can Weak SINDY be used on higher order ordinary differential equations, for example the van der Pol equation (a second order differential equation):
x''(t) - mu (1-x^2) x' + x = 0
(where x' is the first order time derivative, x'' is the second order time derivative, and mu is a constant),
such that we only give the time series of x in the training data and Weak SINDy gives back either a second order differential equation or two first order differential equations for x and x'?

There is an option "derivative_order" in the weak SINDy library but I guess it is only meant for the order of the spatial derivatives in feature library.

If I have a time series of x given to me and I calculate x' (lets name it y) using some numerical methods and give the time series of both x and x' in the training data, Weak SINDy is correctly able to identify the van der Pol equation in the form of coupled set of two first order differential equations:
x' = y
y'= mu (1-x^2) x' - x

So if I have a time series x as a solution of a fourth order differential equation, say,
x''''(t) = f(x,x',x'',x''')
do I need to calculate the numerical derivatives of x up to third order and give x,x',x'',and x''' in the training data to weak SINDy so that it can identify the coupled set of 4 first order differential equations in x,x',x'', and x'''?

Thanks very much in advance.

Best regards,
Bhumika

[akaptano]

Hi Bhumika,

The quick answer is that yes, you can do this, and the easiest way is to convert everything into a coupled system of 1st-order ODEs. In that case, you can repeat exactly what you did with the Van Der Pol equation, and indeed you need to have x, x', x'', and x''' available. Then you could integrate each of these (integrating by parts the derivative terms) and stack them together as [weak_form(x), weak_form(x'), ...] and fit a system of 1st-order equations.

This is rather laborious so instead I have made a few tiny changes to the code so that you can fit a fourth order equation like this directly from data in x, including allowing for time derivatives in the library. I have pushed these changes to a new branch of the code weak_form_temporal_derivatives, and let me know if you have any issues with obtaining the code.

Here is an example that should work correctly with the new code:

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import trapezoid
from scipy.interpolate import RegularGridInterpolator

import pysindy as ps

# integration keywords for solve_ivp, typically needed for chaotic systems
integrator_keywords = {}
integrator_keywords['rtol'] = 1e-12
integrator_keywords['method'] = 'LSODA'
integrator_keywords['atol'] = 1e-12

def convert_u_nth_order_integral(u, weak_pde_library, n):
    K = weak_pde_library.K
    gdim = weak_pde_library.grid_ndim
    u_dot_integral = np.zeros((K, u.shape[-1]))
    deriv_orders = np.zeros(gdim)
    deriv_orders[-1] = n
    w_diff = weak_pde_library._smooth_ppoly(deriv_orders)
    for j in range(u.shape[-1]):
        u_interp = RegularGridInterpolator(
            tuple(weak_pde_library.grid_pts), np.take(u, j, axis=-1)
        )
        for k in range(K):
            u_new = u_interp(np.take(weak_pde_library.XT, k, axis=0))
            u_dot_integral_temp = trapezoid(
                w_diff[k] * u_new,
                x=weak_pde_library.xtgrid_k[k, :, 0],
                axis=0,
            )
            for i in range(1, gdim):
                u_dot_integral_temp = trapezoid(
                    u_dot_integral_temp, 
                    x=weak_pde_library.xtgrid_k[k, :, i], axis=0
                )
            u_dot_integral[k, j] = u_dot_integral_temp
    return u_dot_integral

# Imagine system is x'''' = - 2 * x'' - 1 * x
# Set x' = y, x'' = z, x''' = w
# System of 1st order ODEs becomes:
# x' = y
# y' = z
# z' = w
# w' = - 2 * z - x
def higher_order_system(t, x):
    return [x[1], x[2], x[3], -2 * x[2] - x[0]]

    
# Generate measurement data
dt = 0.01
t_train = np.arange(0, 200, dt)
t_train_span = (t_train[0], t_train[-1])
u0_train = np.random.rand(4)
u_train = solve_ivp(higher_order_system, t_train_span, u0_train,
                    t_eval=t_train, **integrator_keywords).y.T
plt.plot(u_train[:, -1])

# Define weak form ODE library
library_functions = [lambda x: x, lambda x: x * x, lambda x, y: x * y]
library_function_names = [lambda x: x, lambda x: x + x, lambda x, y: x + y]
ode_lib = ps.WeakPDELibrary(
    library_functions=library_functions,
    function_names=library_function_names,
    spatiotemporal_grid=t_train,
    derivative_order=3,
    is_uniform=True,
    num_pts_per_domain=100,
    K=200,
)

# Instantiate and fit the SINDy model with the integral of u_dot
optimizer = ps.SR3(
    threshold=0.1, 
    thresholder="l0", 
    max_iter=1000, 
    normalize_columns=True, 
    tol=1e-10
)
model = ps.SINDy(feature_library=ode_lib, 
                 optimizer=optimizer,
                 feature_names=['x'])

# fit on the non-weak form of x and pass the weak form of the x''''
x_train = u_train[:, 0:1]
u_4th_order_derivative = convert_u_nth_order_integral(x_train, ode_lib, 4)
model.fit(x_train, x_dot=u_4th_order_derivative)
model.print()
model.get_feature_names()

Best,
Alan

[akaptano]

@znicolaou I think we should include a flag that decides whether or not to include the temporal derivatives in the PDELibrary and WeakPDELibrary going forward, since Bhumika's case is actually not a particularly uncommon one. I'll add this to your branch whenever I get around to the pull request review.

[bthakur21]

Many thanks Alan @akaptano for the modified code and the example. It is working well. However, I have another question regarding this.

What if the differential equation in the example above had higher powers of derivatives? For example, consider the equation
x'''' = - 2 * (x'')^2 - 1 * x.
In the above example, the library functions are applied only to the variable x and not to its derivatives. How can one apply the library functions to the derivatives also? Thanks very much again.

Best regards,
Bhumika

[akaptano]

To generate those kind of terms, you can use the GeneralizedLibrary class. With this class you can specify two (or more) separate libraries for x, x', x'', etc. so that for instance, your total library has cubic terms with x and up to quadratic terms in x''. You can also tensor all the sub-libraries together to generate mixed-library terms. Check out https://www.youtube.com/watch?v=UA8C8sccU54 and examples are also included in the Example 1 and Example 15 Jupyter notebooks.

[bthakur21]

@akaptano Thanks very much! I looked into the GeneralizedLibrary class but I am facing a few issues in implementing it.
(1) Since the derivatives are generated within the WeakPDE library, I am not sure how to give them as input variables to other libraries, say the polynomial library.

(2) So, I tensored the WeakPDE library with itself (see the example code below) which generated a lot of duplicate terms in the tensored library. To remove the duplicate terms, I used the constrained optimizer but I haven't been able to recover the "correct" model.

(3) Also, when I use the tensored library, I have to take the value of "K" equal to the length of the training data. Otherwise I get the following error (where the length of my training data is 10000 and K=100)

IndexError: boolean index did not match indexed array along dimension 0; dimension is 10000 but corresponding boolean dimension is 100

See below a code for solving the differential equation
x'' = (1- (x')^2)x'-x
where I have tensored the weak PDE library with itself and applied constraints to remove duplicate terms :

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import trapezoid
from scipy.interpolate import RegularGridInterpolator
import pysindy as ps

# integration keywords for solve_ivp, typically needed for chaotic systems
integrator_keywords = {}
integrator_keywords['rtol'] = 1e-12
integrator_keywords['method'] = 'LSODA'
integrator_keywords['atol'] = 1e-12

def convert_u_nth_order_integral(u, weak_pde_library, n):
    K = weak_pde_library.K
    gdim = weak_pde_library.grid_ndim
    u_dot_integral = np.zeros((K, u.shape[-1]))
    deriv_orders = np.zeros(gdim)
    deriv_orders[-1] = n
    w_diff = weak_pde_library._smooth_ppoly(deriv_orders)
    for j in range(u.shape[-1]):
        u_interp = RegularGridInterpolator(
            tuple(weak_pde_library.grid_pts), np.take(u, j, axis=-1)
        )
        for k in range(K):
            u_new = u_interp(np.take(weak_pde_library.XT, k, axis=0))
            u_dot_integral_temp = trapezoid(
                w_diff[k] * u_new,
                x=weak_pde_library.xtgrid_k[k, :, 0],
                axis=0,
            )
            for i in range(1, gdim):
                u_dot_integral_temp = trapezoid(
                    u_dot_integral_temp, 
                    x=weak_pde_library.xtgrid_k[k, :, i], axis=0
                )
            u_dot_integral[k, j] = u_dot_integral_temp
    return u_dot_integral

# Consider the equation: x'' = (1- (x')^2)*x'-x

def Rayleigh_osc(t, x):
    return [x[1], (1-(x[1])**2)*x[1]-x[0]]

    
# Generate measurement data
dt = 0.01
t_train = np.arange(0, 100, dt)
t_train_span = (t_train[0], t_train[-1])
u0_train = np.random.rand(2)
u_train = solve_ivp(Rayleigh_osc, t_train_span, u0_train,
                    t_eval=t_train, **integrator_keywords).y.T
plt.plot(u_train[:, -1])

# Define weak form ODE library
library_functions = [lambda x: x, lambda x, y: x * y]
library_function_names = [lambda x: x, lambda x, y: x + y]
ode_lib = ps.WeakPDELibrary(
    library_functions=library_functions,
    include_bias=True,
    include_interaction = False,
    function_names=library_function_names,
    spatiotemporal_grid=t_train,
    derivative_order=1,
    is_uniform=True,
    num_pts_per_domain=100,
    K=10000,
)


# Tensored library:
concat_lib =  ode_lib * ode_lib *ode_lib


# Constraints to remove the duplicate terms
constraint_zeros = np.zeros((17))
constraint_matrix = np.zeros((17, 27))
indices_list = [3,6,7,9,10,11,12,15,16,18,19,20,21,22,23,24,25]
for i in range(17):
    constraint_matrix[i,indices_list[i]] = 1
    

# Instantiate and fit the SINDy model with the integral of u_dot
optimizer = ps.ConstrainedSR3(
    threshold=0.5,
    nu=1,
    max_iter=1000,
    thresholder='l0',
    constraint_lhs=constraint_matrix,
    constraint_rhs=constraint_zeros,
    constraint_order="target",
    tol=1e-8
)

model = ps.SINDy(feature_library=concat_lib, 
                 optimizer=optimizer,
                 feature_names=['x'])

# fit on the non-weak form of x and pass the weak form of the x''
x_train = u_train[:, 0:1]
u_2nd_order_derivative = convert_u_nth_order_integral(x_train, ode_lib, 2)
model.fit(x_train, x_dot=u_2nd_order_derivative)
model.print()
model.get_feature_names()

[akaptano]

Looks like there are more bugs in the Tensor/GeneralizedLibrary than I anticipated. Given that, I think your easiest route forward is to convert this into a system of 1st-order ODEs and use that with the weak form. I have push a few changes to get things working on the weak_form_temporal_derivatives and you should be able to correctly run the following:

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import trapezoid
from scipy.interpolate import RegularGridInterpolator
import pysindy as ps

# integration keywords for solve_ivp, typically needed for chaotic systems
integrator_keywords = {}
integrator_keywords['rtol'] = 1e-12
integrator_keywords['method'] = 'LSODA'
integrator_keywords['atol'] = 1e-12


# Consider the equation: x'' = (1- (x')^2)*x'-x
def Rayleigh_osc(t, x):
    return [x[1], (1-(x[1])**2)*x[1]-x[0]]

    
# Generate measurement data
dt = 0.01
t_train = np.arange(0, 100, dt)
t_train_span = (t_train[0], t_train[-1])
u0_train = np.random.rand(2)
u_train = solve_ivp(Rayleigh_osc, t_train_span, u0_train,
                    t_eval=t_train, **integrator_keywords).y.T

# Define custom functions up to cubic terms
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]

# Library for x and x_dot, note omitting the bias term
ode_lib = ps.WeakPDELibrary(
    library_functions=library_functions,
    include_bias=False,
    include_interaction = False,
    function_names=library_function_names,
    spatiotemporal_grid=t_train,
    derivative_order=0,
    is_uniform=True,
    num_pts_per_domain=200,
    K=1000,
)

# Initialize the default inputs, i.e. each library
# uses all the input variables
inputs_temp = np.tile([0, 1], 2)
inputs_per_library = np.reshape(inputs_temp, (2, 2))

# Only use one of the inputs for each library
inputs_per_library[0, 1] = 0
inputs_per_library[1, 0] = 1
print(inputs_per_library)
full_lib = ps.GeneralizedLibrary([ode_lib, ode_lib], 
                                 tensor_array=[[1, 1]],
                                 inputs_per_library=inputs_per_library)

# Instantiate and fit the SINDy model with the integral of u_dot
optimizer = ps.SR3(
    threshold=0.05,
    nu=1,
    max_iter=1000,
    thresholder='l0',
    tol=1e-8
)

model = ps.SINDy(feature_library=full_lib, 
                 optimizer=optimizer,
                 feature_names=['x', 'x_dot'])

# fit on the non-weak form of x and pass the weak form of the x''
model.fit(u_train)
model.print()
model.get_feature_names()

[akaptano]

@znicolaou Just a heads up for the weak formulation update... currently the TensoredLibrary does not allow for weak terms (and therefore neither does the GeneralizedLibrary if any of the libraries are tensored together) and this would be a good fix. The fix is straightforward... basically just need to make sure "n_samples" is correct, since this number changes in the weak form. I made a hack-ish fix to this branch so that it works for Bhumika's application.

[bthakur21]

Thanks very much Alan @akaptano for your help and the fix. I really appreciate it!

[akaptano]

Very happy to help, and please let me know if you find any other issues. This portion of the code is fairly buggy but hopefully you can make do with the type of code I showed above :)

[bthakur21]

Sure, I will let you know. Thanks again.

[akaptano]

This is now implemented in the most recent branch!