question about second-order ode identification

[xiongxin9000]

Hi,
I would like to identify second-order ode below:
$$(R 1)^{\prime}=1.000 R 2 $$

$$ (R2)' = \frac{3.6521}{R_1 \cdot \ln\left(\frac{1}{R_1}\right)} + \frac{0.0051}{R_2 \cdot \ln\left(\frac{1}{R_2}\right)} - \frac{34.8451}{R_1^{2} \cdot \ln\left(\frac{1}{R_1}\right)} + \frac{339.376 \cdot R_2}{R_1^{2} \cdot \ln\left(\frac{1}{R_1}\right)} - \frac{2287.1912^{2}}{R_1 \cdot \ln\left(\frac{1}{R_1}\right)} - \frac{0.242 \cdot R_1^{2} \cdot R_2^{2}}{\ln\left(\frac{1}{R_1}\right)} - \frac{732.746 \cdot R_2^{2}}{R_1} $$

but I whatever I use customised or generalized library, it always comes with R1 and R2 term in second equation which I don't want or that two terms disappear but the first equation is not what I want with $(R 1)^{\prime}=1.000 R 2 $

[Jacob-Stevens-Haas]

Can you share the code to create the generalized or customized library? I'm not sure if this is a coding issue or a data/regression issue.

Also, FWIW, you called this a PDE, but I think you mean ODE since there are no spatial derivatives in your equation.

[xiongxin9000]

Thank you for your reply, please see my code below.

https://github.com/xiongxin9000/sindy.git

[Jacob-Stevens-Haas]

Typically, GeneralizedLibrary creates a composite library with different inputs for each library. However, SINDy will use the complete library for every equation. There is the option to set some library coefficients to zero so that you don't discover those terms in different equations. Look at the ConstrainedSR3 optimizer.

[xiongxin9000]

Hi Jacob,

Thank you for you reply, please let me clarify my question.
what I want is some thing similar as follows,
$$(R 1)^{\prime}=1.000 R 2 $$

$$ (R2)' = \frac{3.6521}{R_1 \cdot \ln\left(\frac{1}{R_1}\right)} + \frac{0.0051}{R_2 \cdot \ln\left(\frac{1}{R_2}\right)} - \frac{34.8451}{R_1^{2} \cdot \ln\left(\frac{1}{R_1}\right)} + \frac{339.376 \cdot R_2}{R_1^{2} \cdot \ln\left(\frac{1}{R_1}\right)} - \frac{2287.1912^{2}}{R_1 \cdot \ln\left(\frac{1}{R_1}\right)} - \frac{0.242 \cdot R_1^{2} \cdot R_2^{2}}{\ln\left(\frac{1}{R_1}\right)} - \frac{732.746 \cdot R_2^{2}}{R_1} $$

but for the first equation it should be $$(R 1)^{\prime}=1.000 R 2 $$ no other term. however, with customsied equation, the other terms ,like second equation, shows up in some condition like below

$$ (R 1)^{\prime}=0.002 R 1+0.578 R 2+-0.5781 /(R 1 * \ln (1 / R 1))+-0.0031 /(R 2 * \ln (1 / R 2))+-0.0011 /\left(R 1^{2} * \ln (1 / R 1)\right)+-0.074 R 1^* R 2^{2} /(\ln(1 / R 1))+0.001 R 2^ 2 / R 1 $$

In terms of the second equation, $R_1,R_2$ term shows up in some condition like below

$$ (R2)' = R_1 + R_2+\frac{3.6521}{R_1 \cdot \ln\left(\frac{1}{R_1}\right)} + \frac{0.0051}{R_2 \cdot \ln\left(\frac{1}{R_2}\right)} - \frac{34.8451}{R_1^{2} \cdot \ln\left(\frac{1}{R_1}\right)} + \frac{339.376 \cdot R_2}{R_1^{2} \cdot \ln\left(\frac{1}{R_1}\right)} - \frac{2287.1912^{2}}{R_1 \cdot \ln\left(\frac{1}{R_1}\right)} - \frac{0.242 \cdot R_1^{2} \cdot R_2^{2}}{\ln\left(\frac{1}{R_1}\right)} - \frac{732.746 \cdot R_2^{2}}{R_1} $$

I want to cancel out the $R_1, R_2$ terms in the second equation. I play around also with ConstrainedSR3 optimizer, but don't quite understand how exactly set the coefficient and the threshold. when I set the threshold=0. The first equation is not correct. when I set threshold to be 0.5, the first equation is zero, but the second equation is what i want.

I add the new code related to constrainedSR3 in same repository named as R-P-SINDY-GEN-constrained.py as below.
https://github.com/xiongxin9000/sindy.git
named as R-P-SINDY-GEN-constrained.py
Thank you for your help.

[Jacob-Stevens-Haas]

ConstrainedSR3 uses the constraint_lhs and constraint_rhs arguments to set terms equal to zero. For example, let's say you're fitting the equations:

$$ \begin{aligned} \dot x_1 &= \theta_{11} f_1(x_1) + \theta_{12} f_2(x_1) +\theta_{13} f_1(x2) \\ \dot x_2 &= \theta_{21} f_1(x_1) +\theta_{22} f_2(x_1) + \theta_{23} f_1(x2). \end{aligned} $$

But If you want to set certain elements to zero, you need to understand the order that feature_library creates the functions. E.g. to set $\theta_{11}, \theta_{22}=0$, you would need the following constraint_lhs:

constraint_lhs = np.array([[1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0])
constraint_rhs = np.array([0, 0])

This ordering is $\theta_{11} , \theta_{12}, \theta_{13} , \theta_{21} , \theta_{22} , \theta_{23} $.

[xiongxin9000]

Thanks for reply. And if I want to make both $\theta_{11}, \theta_{12}$ in equation 1 to be zero. How can achieve that?

[xiongxin9000]

And I managed to define what I wanted about the customized library but it gave me the "ConvergenceWarning: SR3._reduce did not converge after 30 iterations." what would be the issue?

[Jacob-Stevens-Haas]

Thanks for reply. And if I want to make both θ11,θ12 in equation 1 to be zero. How can achieve that?

We don't have a great way to identify the names of equations before fitting the data. The best way is to fit the data once without constraints, then check model.feature_library.get_feature_names(input_features=model.feature_names) to identify the order of feature names, and then fit the model again with a ConstrainedSR3 optimizer with the correct constraints.

[Jacob-Stevens-Haas]

And I managed to define what I wanted about the customized library but it gave me the "ConvergenceWarning: SR3._reduce did not converge after 30 iterations." what would be the issue?

Try increasing the max_iter parameter.

[xiongxin9000]

Hi Jacob,

My code get stuck when using simulate function and after I manually stop it, it give me "capi_return is NULL
Call-back cb_f_in_lsoda__user__routines failed". Also I meet with error that "Input X contains infinity value" which I can catch that exception but for the first error it just stuck and cannot be caught by exception, My question is how can I catch the first error and contine the execution of the code and why I experience the second error since my input is confirmed.

Thanks in advance
Xin

[Jacob-Stevens-Haas]

The infinity value is because you have discovered an unstable model and the integration blows up. You have to catch it by using a warnings filter to turn Warning into Exception.

Also, see scipy/scipy#15940.