Understanding NDSolve::ndmmc

When running NDSolve with Method->{"EquationSimplification" -> "MassMatrix"} I frequently encounter this error:

NDSolveValue::ndmmc: The mass matrix SparseArray[…] should be a constant numerical matrix. The option Method->{"EquationSimplification"->"Residual"} can be used to solve the equations as a system of differential-algebraic equations.

Which is quite annoying because there is zero documentation on this topic, so I have no idea how to debug the code. In fact, what is a "constant numerical matrix" is not even clear.

I've searched across Mathematica's documentation and this is the closest I've found:

https://reference.wolfram.com/language/ref/ConstantArray.html

But it is unlikely that a constant numerical matrix refers to a matrix of only 1 type of element, especially for equation solving.

I would therefore like to ask

1. What is a constant numerical matrix?
2. Why does NDSolve fail to generate constant numerical matrices under Method->{"EquationSimplification"->"MassMatrix"} under certain conditions?
3. In general, what can we do to resolve such error?

Thanks in advance to all the help offered!

Indeed, MassMatrix method isn't well documented, but this warning isn't hard to understand.

According to the tutorial Numerical Solution of Differential-Algebraic Equations, MassMatrix is a method that transforms the equation(s) to the form $$M.x'=F(x,t)$$, but what's undocumented is, the $$M$$ matrix must be a constant numerical matrix. What's constant numerical matrix? It's better to explain with a minimal example. Suppose we want to solve

$$t\ x'(t)-1=0,\ x(2)=3,\ t\in[2,4]$$

with MassMatrix method, how should we transform the ODE to the required $$M.x'=F(x,t)$$ form? Clearly the most straightforward transformation would be $$(t).(x'(t))=(1)$$ ({{t}}.{x'[t]} == {1}) i.e. $$M=(t)$$ (M = {{t}}). This is exactly what happens in NDSolve:

asso = NDSolve[{t x'[t] - 1 == 0, x[2] == 3}, x, {t, 2, 4},
Method -> EquationSimplification -> MassMatrix] // EvaluationData


NDSolve::ndmmc

asso["MessagesExpressions"][[1, 1, -1, 1]] // Normal
(* {{t}} *)


Remark

If I have to guess, NDSolve might have used something similar to

CoefficientArrays[{t x'[t] - 1 == 0}, {x'[t]}] // Last


to generate the mass matrix, but TraceInternal doesn't confirm it.

The $$M$$ involves symbolic variable $$t$$ (t) so it's not numerical. "But once the numeric ODE/DAE solver starts working, $$t$$ will be numerical at every time step!" Yeah, but it varies as time changes i.e. the matrix won't be constant. I'm not knowledge enough to explain why MassMatrix is designed not to allow such matrix, though.

"OK, then how to circumvent this? " In many cases, we can simply help the coefficient extractor a bit. For example, the example above can be transformed to

NDSolve[{x'[t] - 1/t == 0, x[2] == 3}, x, {t, 2, 4},
Method -> EquationSimplification -> MassMatrix]


"Wait, isn't this ridiculous? Why not directly set EquationSimplification -> Solve? That'll automatically transform the equation to x'[t] == 1/t!" Yeah, MassMatrix method doesn't show any advantage for this toy example. But when the system gets large (this happens quite a bit when symbolically discretizing PDE), the EquationSimplification -> Solve may be too expensive for symbolic pre-processing, while EquationSimplification -> Residual may be too weak to handle the resulting system, then MassMatrix method becomes an appreciable method. That's the reason I sometimes set this option when solving ODE systems generated by pdetoode. Typical example:

Why does NDSolve fail to solve the PDEs and spit out mconly warning?

Remark

For those who feel the toy system $$(t).(x'(t))=(1)$$ with only one equation inside uncomfortabe, here's a slightly more complicated toy example:

NDSolve[{{(t + 1) (x'[t] + y'[t]) == (t + 1) (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, y, z}, {t, 0, 1},
Method -> {"EquationSimplification" -> "MassMatrix"}]


NDSolve::ndmmc

Notice this is essentially the documented example for MassMatrix, except that I've multiplied each side of the first equation by (t + 1), yet this is enough to fool the pre-processor of MassMatrix.