odeint-v2
odeint-v2 copied to clipboard
headmyshoulder
Solve ODEs in the form [M] {x}' = {F}
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.
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.
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.
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.
One could also solve for y in My = F(x) using an iterative solver. Depending on M this could be more memory efficient.
One could also solve for
yinMy = F(x)using an iterative solver. Depending onMthis 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.
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 [email protected] 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.
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.
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
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.
I'm currently using LAPACKE for LU decomposition. I'll try to 'combine' ODEINT and LAPACKE. Thanks for the answer.