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)?

  • $\begingroup$ Why not invert the matrix? $\endgroup$ – Yuriy S Oct 10 '18 at 14:57
  • 2
    $\begingroup$ 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? $\endgroup$ – Carl Christian Oct 10 '18 at 14:58
  • $\begingroup$ 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. $\endgroup$ – Chen Ziv Oct 10 '18 at 15:20
  • $\begingroup$ 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. $\endgroup$ – Carl Christian Oct 11 '18 at 9:59
  • $\begingroup$ 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. $\endgroup$ – Chris Rackauckas Oct 11 '18 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.