SciML / SteadyStateDiffEq.jl

Solvers for steady states in scientific machine learning (SciML)
Other
30 stars 22 forks source link

Throw error only in case of singular mass matrix #17

Closed MatFi closed 4 years ago

MatFi commented 4 years ago

If the mass matrix is non-singular the the only solution to M*dn/dt = 0 is the trivial one. In this case solving f(u) = 0 leads to the correct steady state solution u.

MatFi commented 4 years ago

I'm thinking about extending the argumentation to to DAE's. In such a case the mass matrix will come with some rows of zeros and for the corresponding row index f(u) has to be zero as well. So we can ignore a row and column of the mass matrix with that index and only have to test if the rest matrix has full rank. If this is the case f(u)=0 is still the SS-solution for the DAE even if the mass matrix is singular.

ChrisRackauckas commented 4 years ago

Don't you need to make sure the algebraic variables satisfy the condition?

MatFi commented 4 years ago

In deed I missed that point (this happens when a experimental physicist tries to use math). So but can we say that if no time derivative of an algebraic variable exists ( saying that the column of the mass matrix with that index is all zero ) its the case? . There need to be literature on this.

ChrisRackauckas commented 4 years ago

Use ArrayInterface.issingular which specializes a little bit (rank is slower since it uses SVD)

luap-pik commented 4 years ago

I don't see how a (constant) singular mass matrix is a problem for nlsolve when I have on ODE in mass matrix form. This case should not error.

Take for instance an M=Diagonal([1, 1, 0]). As mentioned above, nlsolve solves for a point u0 that fulfills f(u0) = 0. But then, f(u0) automatically lies on the constraint manifold or not? I used nlsove previously on ODEFunctions with mass matrix, returning valid fixed points.

ChrisRackauckas commented 4 years ago

When there is a mass matrix, and then when we find f(u,p,t) = 0, then M*u' = f(u,p,t) is then always satisfied with u' = 0, so it is a steady state. I agree, so we could just completely remove the restriction. That's not necessarily the unique solution there, but I think it's better than erroring out. Given this, I'd accept a PR to remove the restriction.

luap-pik commented 4 years ago

You're right, I created #22.