Skip to main content
2 of 4
added another idea
VoB
  • 560
  • 4
  • 14

Mass Matrix and how to handle it (ODEs) - References

I'm interested in a good reference/paper about how to handle numerically a mass-matrix system as

\begin{align} \mathbf{M}(t,y)\dot{y} =F(y,t) \end{align}

I know that such a problem can be solved by using classical ODEsolvers in MatLab such as ode15s, and others. Of course, due to copiright issues, it's not possible to read the source code and, IMHO, there's a lack of reference.

I've found something interesting in Hairer & Wanner's books (both volumes), but honestly there's not so much stuff. In such problems, the non-singularity of $M(t,y)$ plays a fundamental role, and I'm quite interested in some techniques.


Suppose to have a not-time-dependent and non-singular mass matrix $M(t,y)=M(y)$ and a system \begin{align} M(y)\dot{y} =F(y,t) \end{align}

and say I want to compute the numerical solution with Backward Euler (or any implicit method, in order to avoid numerical instabilities)

So, one can write $\dot{y} = M^{-1}(y) f(y,t):= \tilde{F}(y,t)$ and formally we have

\begin{align} y_{n+1}=y_n + \Delta t \tilde{F}(t_{n+1},y_{n+1})=y_n+ M(y_{n+1})^{-1} \Delta t F(t,y_{n+1}) \end{align}

and hence

\begin{align} M(y_{n+1}) [y_{n+1}-y_n]=\Delta t F(t_{n+1},y_{n+1}) \end{align}

Another interesting problem is how to solve with Newton's nethod this non-linear system of equations. As far I know, a simplified Newton's method is used, but I can't found any reference about it.

I'd be very grateful if someone has references/hints or something else !

EDIT (after Bill's comment) [REMARK: I'd like not to copy explicitely MatLab's approach] For the simplfied Newton's method, a possible approach could be the following.

Starting from the functional we have to set to zer written above, I'd like to differentiate w.r.t. $y_{n+1}$, namely

\begin{align} J_n=\frac{\partial (M(y_{n+1}) \cdot y_{n+1})}{\partial y_{n+1})}- \Delta t \frac{\partial F(t_{n+1},y_{n+1})}{\partial y_{n+1}} \end{align}

and consider the constant matrix $M(y_{n+1})$ at each step, so the jacobian would become

\begin{align} J_n=M(y_{n+1})-\Delta t \frac{\partial F(t_{n+1},y_{n+1})}{\partial y_{n+1}} \end{align}

Could it be a good starting point?

VoB
  • 560
  • 4
  • 14