Application to a different system of ODEs

[Ricvalp]

Hi,

I was trying to apply this method to a different system od ODEs, namely the 2D Hopf normal form defined by:

function Hopf2D(du, u, p, t)
μ, ω, A1, A2 = p
du[1] = μ*u[1] - ω*u[2] + A1*u[1]*(u[1]^(2) + u[2]^(2))
du[2] = μ*u[2] + ω*u[1] + A2*u[2]*(u[1]^(2) + u[2]^(2))
end

where I approximate the third term of the right-hand side of the equations with a UA L(u, p).

Even though it seems that I get good convergence (at least visually) between the (unavailable) ideal data L̄ and L̂ = L(Xₙ, res2.minimizer), the SiNDy algorithm does not output the correct linear combination of u_1 and u_2 (namely A1u_1*(u_1^(2) + u_2^(2) and A2u_2*(u_1^(2) + u_2^(2)).

I was wondering whether the part of the code

opt = SR3()
λ = exp10.(-7:0.1:3)
g(x) = x[1] < 1 ? Inf : norm(x, 2)

is in some way "specific" for the Lotka-Volterra problem.

Or are there any other parts of the code that cannot be applied to some similar systems of ODEs? Am I missing something?

Here's the trainined UA L̂, together with the correct L̄:

The output of the SiNDy is:

2 dimensional basis in ["u₁", "u₂"]
du₁ = u₁ ^ 3 * p₁
du₂ = p₂ * u₂

[AlCap23]

Hey there!

Generally speaking, the methodology is universal :) . However, given the structure of the equations to recover and the magnitude and curvature of the training trajectory the result may vary. The snippet

opt = SR3()
λ = exp10.(-7:0.1:3)
g(x) = x[1] < 1 ? Inf : norm(x, 2)

Can indeed be the reason why this is failing on your specific example. You can try out different λs, which is the threshold of the sparse regression ( meaning it removes candidates under a specific λ ). g(x) is the selection criterion for the pareto front optimization, so it selects the pareto optimal solution based on this specific function.

So here are some rough guidelines I can give you without seeing the parameters etc

• Try different initial conditions, simulations times and savetimes
• Try different thresholds and switch off normalization
• Try different g(x), so maybe use something like g(x) = x[1] < 1 ? Inf : norm([0.1; 10.0] .* x, 2) which results in a higher penalties for the reconstruction error

[ChrisRackauckas]

How well does it extrapolate? I'd be curious to see.

I don't think that's really the issue. The problem is that it both found a sparse solution and it seems to fit well. My guess is that, given the data it has seen, it's unidentifiable between those two equations and you'd need more data. We're still working to try and find out how to establish when that could be the case, but this is probably a good equation to look into it (so thanks!)

[Ricvalp]

Here's how well the SiNDy linear combination fits the trained UA:

• Training tspan = (0.0f0, 4.0f0), solve's saveat = 0.1

• Training tspan = (0.0f0, 4.5f0), solve's saveat = 0.05

It indeed worked better with a larger training set... still one of them doesn't fit well and I get an Instability Warning when I solve the resulting system of ODEs.

PS: I tried with @AlCap23 suggestion for g(x) but SiNDy did not converge to a sparse solution...

PPS: The Jupyter notebook is in this folder if you want to take a look. Thank you both for the prompt reply!

[ChrisRackauckas]

Oh, I thought you were showing the sparse regression fit. This shows that the sparse regression fit isn't good at all, and that's a SINDY problem.