Open
Description
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!
Metadata
Metadata
Assignees
Labels
No labels