Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
matheusbrf opened this issue Sep 19, 2018 · 10 comments
Open

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

matheusbrf opened this issue Sep 19, 2018 · 10 comments

Comments

@matheusbrf
Copy link

matheusbrf commented Sep 19, 2018

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
Copy link
Collaborator

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
Copy link
Author

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
Copy link
Collaborator

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
Copy link
Contributor

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

@mariomulansky
Copy link
Collaborator

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
Copy link
Owner

headmyshoulder commented Sep 23, 2018 via email

@ddemidov
Copy link
Contributor

ddemidov commented Sep 23, 2018

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
Copy link
Author

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
Copy link
Contributor

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
Copy link
Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants