What's the solver behind Method -> {“EquationSimplification” -> “MassMatrix”}

Question

Related: What's behind Method -> {“EquationSimplification” -> “Residual”}

By reading the document carefully, I've learnt that there are 3 possible options in "EquationSimplification" for DAE systems: Solve, Residual and MassMatrix. Solve and Residual have been explained in the above answer detailedly. As far as I understand, Residual converts the DAE system into the residual form $f(t, x, \dot{x}) = 0$ and solves the system with IDA from Sundials.

On the other hand, MassMatrix converts the system into the mass matrix form $M.\dot{x} = f(t, x)$. But I cannot find in the documentation which solver will be used (maybe I missed something?). Can anyone explain what's the solver behind it?

I tried some examples with NDSolve`ProcessEquations with "EquationSimplification" -> "MassMatrix" and "TimeIntegration" -> "StateSpace". The NumericalFunction behind them are different. My guess is that it does not use the method StateSpace. It probably uses some ODE solvers, since other softwares like Matlab does use ODE solvers to solve DAE problems in mass matrix form?

Answer 1

I tried some examples with NDSolve`ProcessEquations with "EquationSimplification" -> "MassMatrix" and "TimeIntegration" -> "StateSpace". The NumericalFunction behind them are different.

I'm not sure how you checked their equivalence. (Notice Experimental`NumericalFunction involves temporal variables so you cannot check the equivalence with ===. ) But Method -> {EquationSimplification -> MassMatrix, TimeIntegration -> StateSpace} seems to be the answer. This can be checked with

test[method_] := 
 NDSolveValue[{{x'[t] + y'[t] == Sin[t] - z'[t], 
                y'[t] + z'[t] == Sin[t] + x[t], 
                x'[t] + z'[t] == y[t] + Cos[t]}, 
               {x[0] == 1, y[0] == 1, z[0] == 1}}, 
              x, {t, 0, 1}, 
              Method -> {EquationSimplification -> MassMatrix, 
                         TimeIntegration -> method}]

test[Automatic] === test[StateSpace]
(* True *)

Further check shows the following methods can also be combined with MassMatrix method, but none of them reproduces the default output:

methodlist = {{"DoubleStep", Method -> Automatic}, 
              {"Extrapolation", Method -> Automatic}, 
              {"FixedStep", Method -> Automatic}, 
             "LinearlyImplicitEuler", 
             "LinearlyImplicitMidpoint", 
             "LinearlyImplicitModifiedMidpoint", 
             "ImplicitRungeKutta"};

test /@ methodlist // Map@SameAs@test@Automatic
(* {False, False, False, False, False, False, False} *)