Solve ODEs in the form [M] {x}' = {F}

[matheusbrf]

Hello,

I'm trying to convert a code I wrote in MATLAB to C++. It is about the dynamics of a structural beam. Because of the formulation, I naturally get a set of coupled ODEs in the form (already written in a state space representation):

[M] {x}' = {F}

In which [M] is a constant matrix, {x}' is the state-vector and {F} is a nonlinear vector. A priori, I can invert the [M] matrix (which is non-singular) and get a system in the standard form

{x}' = {G}

Are there some way to integrate the equations without inverting the [M] matrix? This is, I want to solve the system directly in the form [M] {x}' = {F}. I'm new with C++. I've read all the examples available in the documentation, but due to my lack of experience with C ++, I couldn't get a response. Thank you in advance.

[mariomulansky]

Differential equations of the form M x' = F(x) are, in general, not ODEs - so odeint does not provide solvers for this. If M is constant and non-singular, the transformation to the ODE: x' = M^-1 F(x) seems the best, most efficient way to go I think. Otherwise you are dealing with an algebraic differential equation which is considerably harder than an ODE.

[matheusbrf]

Since [M] is always invertible, I'm dealing with a system of ODEs, not DAEs. Are there some way to deal with this problem in the form [M] {x}' = {F} using odeint? Thanks for answering.

[mariomulansky]

Well the only way I see is to invert M and treat it as an ODE. You can use Eigen to invert M, for example.

[ddemidov]

One could also solve for y in My = F(x) using an iterative solver. Depending on M this could be more memory efficient.

[mariomulansky]

One could also solve for y in My = F(x) using an iterative solver. Depending on M this could be more memory efficient.

Right, but then you would need to solve in every Runge-Kutta step (or whatever ode solver is used), wouldn't you? While otherwise you compute M^{-1} only once for the whole integration.

[headmyshoulder]

Why an iterative solver? If M is small enough, LU might be sufficient and the decomposition can be computed only once.

I also think that the precision is better you it is solved by LU instead of the using the inverse, but maybe this is not important, since a solution of an ODE is approximated and the precision is below the error of the approximation. But performance might suffer.

On 22. Sep 2018, at 19:49, Mario Mulansky notifications@github.com wrote:

One could also solve for y in My = F(x) using an iterative solver. Depending on M this could be more memory efficient.

Right, but then you would need to solve in every Runge-Kutta step (or whatever ode solver is used), wouldn't you? While otherwise you compute M^{-1} only once for the whole integration.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/231#issuecomment-423761468, or mute the thread https://github.com/notifications/unsubscribe-auth/AAzSlPRTanMQ9sMMoGGCw2qmpS1gpvwwks5udng0gaJpZM4WxFVC.

[ddemidov]

That was just a suggestion in case the system is too large to be solved by a direct solver (as compute cost in general grows as O(n^3) for a matrix inverse, and an inverse of a sparse matrix in general may be not sparse). I agree that if M is small enough then LU decomposition is the way to go.

[matheusbrf]

Thank you all. What does 'small enough' mean? The [M] matrix is about 10,000x10,000. Is this too large for LU decomposition? I was thinking of evaluating [M]^-1 using LU decomposition itself. Could it be a good idea? Thank you all again!!! A lot of learning for me

[ddemidov]

10000 could be too big for a dense LU decomposition, but you could use a sparse LU solver, from SuperLU or PaStiX. You could also borrow an implementation of SkylineLU solver here.

[matheusbrf]

I'm currently using LAPACKE for LU decomposition. I'll try to 'combine' ODEINT and LAPACKE. Thanks for the answer.