# How to solve an implicit ODE with forward Euler?

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 '18 at 14:57
• 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 '18 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 '18 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 '18 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 '18 at 14:35