Consider the implicit ODE

$$ M(y)\dot{y} = F(t,y) $$

If $M$ is non-singular for all $y$

How to use the forward-Euler method to numerically solve for $y$ without inverting $M(y)$?

I only came out following method:

  1. Forward Euler, $y_{n+1} = y_n + h \dot{y_n}$
  2. Multiply by $M(y_n)$ to obtain $M(y_n) y_{n+1} = M(y_n) y_n + h M(y_n) \dot{y_n}$
  3. $M(y_n) y_{n+1} = M(y_n) y_n + h F(t, y_n) = \mathrm{RHS}$
  4. Solve $M(y_n) y_{n+1} = \mathrm{RHS}$ by iterative matrix solver from start point $y_n$ and get $y_{n+1}$.

But I think the procedure is not efficient since it needs to iterative solve $M(y_n) y_{n+1} = \mathrm{RHS}$.

Could anyone tell me a practical or more efficient method (need to use forward-Euler)?

  • Why not invert the matrix? – Yuriy S Oct 10 at 14:57
  • 2
    In general, you can not avoid solving a linear system with coefficent matrix $M(y_n)$. Perhaps you have a special case and some additional information which at first glance appeared irrelevant?. In particular, what can be said about the matrix $M$: size, sparsity, mathematical properties? – Carl Christian Oct 10 at 14:58
  • To Yuriy, since I want to keep the advantage of explicit method that does not need to factorize the matrix $I-\alpha \frac{\partial F}{\partial y}$ when solve ODE unlike other implicit method, so it's need to prevent to invert $M\\$. To Carl, so the only method with forward-Euler without invert $M$ is like what I proposed (i.e. need to use iterative matrix solver)? In addition, the $M$ is large sparse matrix with about $10^6$ rank. – Chen Ziv Oct 10 at 15:20
  • I imagine that it would be slightly cheaper to solve $M(y_n)z_n = F(t_n,y_n)$ for $z_n$ and then compute $y_{n+1} = y_n + h z_n$, but this is a very minor point. What matters most is how you do the solve. If a direct sparse solver fails due to excessive fill-in, you must try to reorder the system. If this is not successful, then you should combine reordering with sparse preconditioners. As the systems change only slowly, it is likely that a Krylov method like GMRES will benefit from subspace recycling. Moreover, you may be able to extrapolate a good initial guess from past values. – Carl Christian Oct 11 at 9:59
  • If this is coming from an FEM discretization, this is why people use mass lumping to change said mass matrix to a diagonal matrix. In general, there isn't a way to solve this without a linear solve. But at that point, implicit methods already need to factorize $M - \alpha J$, so using an implicit solver usually just is a better idea (if it's from a parabolic PDE discretization where stiffness will invariably occur). Otherwise, you could do an implicit method and just linear solve on $M$ each step. – Chris Rackauckas Oct 11 at 14:35

Your Answer


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.