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

Question

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.

Answer 1

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