# 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 method this non-linear system of equations. As far as I know, a simplified Newton's method is used, but I can't find 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?

EDIT$$^2$$ [02/03/19']

I tried to use the previous approach to solve simple autonomous and non-singular problems like

\begin{align} \begin{bmatrix} y(1) & 0 \\ 0 & y(2) \end{bmatrix} \dot{y} = [y(1),\sin(y(2))] \end{align}

with $$y(1)(0)=y(2)(0)=1$$ and it worked in the right way.

Also with other non-linear and more complicated problems, the numerical solution was ok.

I also tried to generalize this approach to the trapezoidal rule: startinf from $$M(y) \dot{y}=f(y)$$, in the hypotesis that $$M$$ is non-singular, I have \begin{align} \dot{y}= M^{-1}(y)f(y)=\tilde{f}(y) \end{align}

So, the standard trapezoidal rule leads to : \begin{align} y_{n+1}=y_{n} +\frac{\Delta t}{2}(M^{-1}(y_{n+1})f(y_{n+1})+M^{-1}(y_n)f(y_n)) \end{align} and multiplying first for one and then for the other inverse, one gets

\begin{align} M(y_n)M(y_{n+1})[y_{n+1}-y_n]=\frac{\Delta t}{2} (M(y_n)f(y_{n+1})+M(y_{n+1})f(y_n)) \end{align}

The jacobian has been computed with the same logic as before, "freezing" the $$M(y_{n+1})$$ term.

A rapid numerical experiment with the same problem above shows the correct order of convergence.

Any other hints, suggestions, comments, are really appreciated.

• On the contrary, the matlab ode solvers including ode15s are well documented (IMHO). I suggest starting with this reference: mathworks.com/help/pdf_doc/otherdocs/ode_suite.pdf – Bill Greene Feb 28 '19 at 16:47
• @BillGreene you're right! Honestly, I read it some time ago, and probably too much quickly! At page 8 there's an approach which is quite familiar with mine – VoB Feb 28 '19 at 16:53
• Edited my post, adding a possible way to handle Newton's method – VoB Feb 28 '19 at 17:23
• A quick question, don’t you need commutativity between $M(y_n)$ and $M(y_{n+1}$ in your trapezoidal rule derivation (last equation in Edit 2)? Is there a reason to assume that such commutativity holds? – VorKir Mar 2 '19 at 20:11
• With the test problem I used above, the order was right with both the expressions, but your question makes sense: there's nothing that ensures commutativity for that product – VoB Mar 3 '19 at 9:52

Ignoring Newton's method here is the wrong approach! The fact that you're using Newton's method is what makes this cheap to add, and is what makes singular mass matrices possible. Essentially look at Implicit Euler. You'll notice that the only matrix you need to use in your quasi-Newton steps is

$$W = (I - \gamma J)$$

and if you have a mass matrix, this becomes

$$W = (M - \gamma J)$$

This is true for any BDF/SDIRK/Rosenbrock method, since you can simplify any of their rootfinding equations to

$$M u_{n+1} - \gamma f(u_{n+1},p,t) - C = 0$$

where $$C$$ is all of the stuff that doesn't depend on $$u_{n+1}$$. So with that small change to your Newton method, you now converge to the solution with the mass matrix for any of these (semi-)implicit schemes. Newton's method just uses this for a linear solve, so you only ever have to LU it, never inverse. If $$M$$ is singular, you can still have that $$W$$ isn't, and also just use something robust like QR-factorization, and then this will still converge. A singular mass matrix is another representation of a DAE (use the terms without a derivative to define a constraint relation).

(Note that for higher order methods, you may have to handle FSAL separately with an $$M$$ linear solve.)

• The equations you wrote are OK, (I assume $\gamma$ is the step-size $\Delta t$). Looking at mine, I see that the third your equation (the one I have to set to zero) is the same I have written in the first edit, i.e. : $M(y_{n+1}) [y_{n+1}-y_n]=\Delta t F(t_{n+1},y_{n+1})$. So actually the matrix I need in my quasi-Newton's method is $M(y_{n})- \Delta t J$, where $J=\frac{\partial f(t_{n+1},y_{n+1})}{\partial y_{n+1}}$. This seems to be, basically, the same thing you wrote. As you said, I could LU it at each time step in order to save computational cost. Am I wright? – VoB Mar 4 '19 at 8:47
• Well, actual ODE solvers are much more complicated than that. $\gamma$ is proportional to $\Delta t$ and is dependent on the diagonal of the tableau (which happens to be 1 for Implicit Euler), and you don't necessarily have to recompute $W$ every step, but yes the general idea is there that there's no reason to invert $M$ and then $W$ since it can be done simultaneously, and doing it simultaneously as numerically sound and handles the singular $M$ case. – Chris Rackauckas Mar 4 '19 at 10:48