search

dynamicslab / pysindy

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

PDE-finding when control variables and data have same space-time grid (compressible fluid continuity equation) #580

Open wzy7178 opened 3 weeks ago

wzy7178 commented 3 weeks ago

Dear Pysindy developers: I hope to use Pysindy to discover the continuity equation of a compressible fluid ρ_t + (ρ v)_x = 0 from the time and spatial-dependent data ρ(x,t) and v(x,t).

Since the time derivative only acts on ρ_t, I tried to use SINDyCP to input v(x,t) as the control variables. However, it is not clear to me how to do so. I hope to input ρ(x,t) and v(x,t) as two-dimensional arrays with x and t grids so I can use the finite difference method to compute the derivative terms.

Here is what I tried:

import numpy as np
import pysindy as ps

#note that, this data is just for testing the program, as its scale is probably too small to find the correct pde. But I don't know how to attach a larger data set here.

rho=np.array([[0.091552  , 0.09155687, 0.09157203, 0.09160017, 0.09164898,
        0.0917346 , 0.09188373, 0.09213237, 0.09251997, 0.0930796 ,
        0.093827  , 0.09475279, 0.09582145, 0.0969782 , 0.09816153,
        0.09931652, 0.10040419, 0.10140393, 0.10230954, 0.10312231],
       [0.09397726, 0.09400235, 0.09408249, 0.0942304 , 0.09446122,
        0.09478432, 0.09519488, 0.09566922, 0.09616765, 0.09664576,
        0.09707048, 0.09743374, 0.09775648, 0.09808004, 0.09844847,
        0.09889003, 0.09940686, 0.09997703, 0.10056701, 0.10114723],
       [0.09879489, 0.0988553 , 0.0990323 , 0.09931141, 0.09966473,
        0.10005034, 0.10041737, 0.10071607, 0.10090855, 0.10097523,
        0.10091475, 0.10073915, 0.10046871, 0.10012991, 0.0997565 ,
        0.0993903 , 0.09907804, 0.09886317, 0.09877574, 0.09882529],
       [0.10524154, 0.10526484, 0.10531927, 0.10536524, 0.10535576,
        0.1052549 , 0.10504934, 0.10474882, 0.10437702, 0.10395887,
        0.10351054, 0.10303556, 0.10252728, 0.10197501, 0.10137148,
        0.10071921, 0.10003473, 0.09934867, 0.09870055, 0.09812839],
       [0.11043432, 0.11032064, 0.10999391, 0.10949278, 0.10886931,
        0.10817585, 0.10745468, 0.10673352, 0.10602682, 0.10534054,
        0.10467723, 0.10403876, 0.10342609, 0.10283683, 0.10226203,
        0.10168394, 0.10107618, 0.10040721, 0.09964715, 0.09877678],
       [0.11043432, 0.11032064, 0.10999391, 0.10949278, 0.10886931,
        0.10817585, 0.10745468, 0.10673352, 0.10602682, 0.10534054,
        0.10467723, 0.10403876, 0.10342609, 0.10283683, 0.10226203,
        0.10168394, 0.10107618, 0.10040721, 0.09964715, 0.09877678],
       [0.10524154, 0.10526484, 0.10531927, 0.10536524, 0.10535576,
        0.1052549 , 0.10504934, 0.10474882, 0.10437702, 0.10395887,
        0.10351054, 0.10303556, 0.10252728, 0.10197501, 0.10137148,
        0.10071921, 0.10003473, 0.09934867, 0.09870055, 0.09812839],
       [0.09879489, 0.0988553 , 0.0990323 , 0.09931141, 0.09966473,
        0.10005034, 0.10041737, 0.10071607, 0.10090855, 0.10097523,
        0.10091475, 0.10073915, 0.10046871, 0.10012991, 0.0997565 ,
        0.0993903 , 0.09907804, 0.09886317, 0.09877574, 0.09882529],
       [0.09397726, 0.09400235, 0.09408249, 0.0942304 , 0.09446122,
        0.09478432, 0.09519488, 0.09566922, 0.09616765, 0.09664576,
        0.09707048, 0.09743374, 0.09775648, 0.09808004, 0.09844847,
        0.09889003, 0.09940686, 0.09997703, 0.10056701, 0.10114723],
       [0.091552  , 0.09155687, 0.09157203, 0.09160017, 0.09164898,
        0.0917346 , 0.09188373, 0.09213237, 0.09251997, 0.0930796 ,
        0.093827  , 0.09475279, 0.09582145, 0.0969782 , 0.09816153,
        0.09931652, 0.10040419, 0.10140393, 0.10230954, 0.10312231]])

v=np.array([[ 0.00000000e+00, -2.11598576e-04, -4.50226199e-04,
        -7.87760998e-04, -1.36860698e-03, -2.40264144e-03,
        -4.11204891e-03, -6.63899900e-03, -9.94673051e-03,
        -1.37634371e-02, -1.76102266e-02, -2.09192513e-02,
        -2.32024527e-02, -2.42009981e-02, -2.39495898e-02,
        -2.27289995e-02, -2.09344183e-02, -1.89260997e-02,
        -1.69304712e-02, -1.50259830e-02],
       [ 0.00000000e+00, -1.27809508e-03, -2.77769560e-03,
        -4.68527373e-03, -7.10891342e-03, -1.00360184e-02,
        -1.33186090e-02, -1.67094281e-02, -1.99450647e-02,
        -2.28388374e-02, -2.53315724e-02, -2.74675955e-02,
        -2.93082917e-02, -3.08375964e-02, -3.19222826e-02,
        -3.23569079e-02, -3.19703358e-02, -3.07322922e-02,
        -2.87981397e-02, -2.64672430e-02],
       [ 0.00000000e+00, -3.63335027e-03, -7.29883472e-03,
        -1.09364165e-02, -1.43626796e-02, -1.73222011e-02,
        -1.95917433e-02, -2.10775616e-02, -2.18544681e-02,
        -2.21332751e-02, -2.21838623e-02, -2.22578806e-02,
        -2.25410907e-02, -2.31368893e-02, -2.40640538e-02,
        -2.52553632e-02, -2.65613557e-02, -2.77751634e-02,
        -2.86860778e-02, -2.91456653e-02],
       [ 0.00000000e+00, -4.30056357e-03, -8.07412635e-03,
        -1.09465908e-02, -1.27866896e-02, -1.37002917e-02,
        -1.39421835e-02, -1.37951572e-02, -1.34739517e-02,
        -1.30889744e-02, -1.26694449e-02, -1.22193511e-02,
        -1.17731151e-02, -1.14273947e-02, -1.13398058e-02,
        -1.16957008e-02, -1.26496500e-02, -1.42540576e-02,
        -1.63964653e-02, -1.87746294e-02],
       [ 0.00000000e+00, -8.50854122e-17,  5.89760099e-16,
         1.48936070e-15,  1.97322433e-15,  1.67179102e-15,
         1.05379315e-15,  6.15001390e-16,  9.61462630e-16,
         1.65657344e-15,  2.28278528e-15,  2.52559538e-15,
         2.35247843e-15,  2.01566286e-15,  1.66257714e-15,
         1.35840071e-15,  9.23688589e-16,  4.19149945e-16,
         4.72884828e-17,  1.45806526e-16],
       [ 0.00000000e+00,  4.17270431e-03,  7.83836849e-03,
         1.06381601e-02,  1.24470027e-02,  1.33662069e-02,
         1.36382845e-02,  1.35320980e-02,  1.32518360e-02,
         1.29023859e-02,  1.25114762e-02,  1.20838412e-02,
         1.16554394e-02,  1.13244480e-02,  1.12498722e-02,
         1.16189110e-02,  1.25894750e-02,  1.42194136e-02,
         1.64026804e-02,  1.88410780e-02],
       [ 0.00000000e+00,  3.46123685e-03,  6.96276533e-03,
         1.04541054e-02,  1.37622962e-02,  1.66411322e-02,
         1.88728270e-02,  2.03631321e-02,  2.11795458e-02,
         2.15205919e-02,  2.16431667e-02,  2.17895659e-02,
         2.21418806e-02,  2.28049759e-02,  2.38021624e-02,
         2.50713564e-02,  2.64667875e-02,  2.77828746e-02,
         2.88081930e-02,  2.93925074e-02],
       [ 0.00000000e+00,  1.22129926e-03,  2.65514617e-03,
         4.48142553e-03,  6.80639160e-03,  9.62189952e-03,
         1.27902977e-02,  1.60783289e-02,  1.92373318e-02,
         2.20931425e-02,  2.45941262e-02,  2.67856155e-02,
         2.87238938e-02,  3.03833390e-02,  3.16184417e-02,
         3.22106237e-02,  3.19768506e-02,  3.08774387e-02,
         2.90608350e-02,  2.68221282e-02],
       [ 0.00000000e+00,  2.06205653e-04,  4.38534249e-04,
         7.66795138e-04,  1.33138243e-03,  2.33666982e-03,
         4.00037432e-03,  6.46525609e-03,  9.70293796e-03,
         1.34563673e-02,  1.72625942e-02,  2.05655446e-02,
         2.28818534e-02,  2.39493143e-02,  2.37913079e-02,
         2.26717762e-02,  2.09689257e-02,  1.90313645e-02,
         1.70812886e-02,  1.51987475e-02],
       [ 0.00000000e+00, -1.21929100e-15, -2.33477241e-15,
        -2.28218825e-15, -1.01918245e-15,  6.76215627e-16,
         1.65753004e-15,  1.38582347e-15,  3.47212153e-16,
        -5.49764465e-16, -4.65121426e-16,  4.80556770e-16,
         1.35514519e-15,  1.26832931e-15,  1.53881905e-16,
        -9.49645641e-16, -1.13465964e-15, -3.23798612e-16,
         5.67327347e-16,  6.81221338e-16]])

x=np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])  #spatial grid
t=np.array([0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6. ,
       6.5, 7. , 7.5, 8. , 8.5, 9. , 9.5]) #time grid

dt = t[1]-t[0]

rhot = ps.FiniteDifference(axis=1)._differentiate(rho, dt)

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

feature_lib = ps.PDELibrary(
    library_functions=library_functions,
    function_names=library_function_names,
    derivative_order=1,
    spatial_grid=x,
    periodic=True
)

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

parameter_lib = ps.PDELibrary(
    library_functions=library_functions,
    derivative_order=1,
    function_names=library_function_names,
    #include_interaction=False,
    spatial_grid=x,
    #include_bias=True,
    periodic=True
)

lib = ps.ParameterizedLibrary(
    parameter_library=parameter_lib,
    feature_library=feature_lib,
    num_parameters=1,
    num_features=1,
)  #I input both 

opt = ps.STLSQ(threshold=1e-5, alpha=1e-4, normalize_columns=False)
model = ps.SINDy(feature_library=lib, optimizer=opt, feature_names=["rho_feature","v_para"])

model.fit(rho, u=v, x_dot=rhot)  
model.print()

When I run this, I get following output:

(rho_feature)' = -0.001 v_para rho_feature + 0.630 v_para_1 rho_feature_1
(v_para)' = -0.011 v_para rho_feature + 39.841 v_para_1 rho_feature_1 + -210639.142 v_para_1 rho_featurerho_feature_1 + 210248.267 v_parav_para_1 rho_feature_1
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[3], line 170
    167 model = ps.SINDy(feature_library=lib, optimizer=opt, feature_names=["rho_feature","v_para"])
    169 model.fit(rho, u=v, x_dot=rhot)  
--> 170 model.print()

File [/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pysindy/pysindy.py:538](http://localhost:8888/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pysindy/pysindy.py#line=537), in SINDy.print(self, lhs, precision)
    536 elif lhs is None:
    537     if not sindy_pi_flag or not isinstance(self.optimizer, SINDyPI):
--> 538         names = "(" + feature_names[i] + ")"
    539         print(names + "' = " + eqn)
    540     else:

IndexError: list index out of range

I find it confusing. Especially, it outputs the fitted equation for both "rho_feature" and "v_para", but I only want to fit "rho_feature".

As I am new to PySINDy, I would appreciate it if you could explain to me a bit about how to treat such PDE-finding cases where your control parameters [in this case v(x,t)] have the same space-time grid as the data [in this case rho(x,t)]. Thank you!