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?