What's behind Method -> {"EquationSimplification" -> "Residual"}

Question

In order to solve a quite large system of differential equations, I have the habit to use the NDSolve command without changing any options. As I wanted more precision, I increased the number of points to integrate. There I get an error and suggestion to use Method -> {"EquationSimplification" -> "Residual"} in NDSolve... which I did and it works fine now.

Can anyone explain to me how it really works and what's behind this method ?

Answer 1

This post contains several code blocks, you can copy them easily with the help of importCode.

---

As mentioned in the comment above, the answer is hidden in this tutorial. Given the tutorial is a bit obscure, I'd like to retell the relevant part in an easier to understand way.

For illustration, let's consider the following initial value problem (IVP):

$$ 2 y''(x)=y'(x)-3 y(x)-4 $$$$ y(0)=5,\ y'(0)=7 $$

eqn = 2 y''[x] == y'[x] - 3 y[x] - 4;
ic = {y[0] == 5, y'[0] == 7};

We know this can be easily solved by NDSolve / NDSolveValue:

seq = Sequence[{eqn, ic}, y, {x, 0, 6}];
sol = NDSolveValue[seq];
ListLinePlot@sol

Nevertheless, do you know how NDSolve solves this problem? The process is quite involved of course, but the part relevant to your question is, NDSolve needs to transform the equation to some kind of standard form, which is controlled by EquationSimplification.

Currently there exist 3 possible option values for EquationSimplification: Solve, MassMatrix and Residual. By default NDSolve will try Solve first. What does Solve do? It'll transform the equation to the following form:

$$ y_1'(x)=\frac{1}{2}(-3 y_0(x)+y_1(x)-4)$$ $$ y_0'(x)=y_1(x) $$ $$y_0(0)=5,\ y_1(0)=7$$

i.e. right hand side (RHS) of the resulting system consists no derivative term, and left hand side (LHS) of every equation is only a single 1st order derivative term whose coefficient is $1$. This can be verified with the following code:

{state} = NDSolve`ProcessEquations@seq;

state[NumericalFunction][FunctionExpression]
(* Function[{x, y, NDSolve`y$1$1}, {NDSolve`y$1$1, 1/2 (-4 - 3 y + NDSolve`y$1$1)}] *)

Rule @@ (Flatten@state@# & /@ {Variables, WorkingVariables}) // Thread
(* {x -> x, y -> y, y' -> NDSolve`y$1$1, …… *)

The system is then solved with an ordinary differential equation (ODE) solver.

While when Residual is chosen (or equivalently SolveDelayed -> True is set), the equation will be simply transformed to

$$ 2 y''(x)-y'(x)+3 y(x)+4=0 $$$$ y(0)=5,\ y'(0)=7 $$

i.e. all the non-zero term in the equation is just moved to one side. This can be verified by:

{state2} = NDSolve`ProcessEquations[seq, SolveDelayed -> True];

state2[NumericalFunction][FunctionExpression]
(* Function[{x, y, NDSolve`y$2$1, y$10678, 
  NDSolve`y$2$1$10678}, {y$10678 - NDSolve`y$2$1, 
  4 + 3 y - NDSolve`y$2$1 + 2 NDSolve`y$2$1$10678}] *)

Rule @@ (Flatten@state2@# & /@ {Variables, WorkingVariables}) // Thread
(* {x -> x, y -> y, y' -> NDSolve`y$2$1, y' -> y$10678, y'' -> NDSolve`y$2$1$10678} *)

The equation is then solved with a differential algebraic equation (DAE) solver.

Apparently, the former transformation is more complicated. When the equation to be transformed is intricate, or the equation system is very large, or the equation system just can't be transform to the desired form, this transformation can be rather time consuming or even never finishes. The following is a simple way to reproduce the issue mentioned in the question with large system, we just repeat eqn for 3 10^4 times:

$HistoryLength = 0;

sys = With[{n = 3 10^4}, 
   Unevaluated@Sequence[Table[{eqn, ic} /. y -> y@i, {i, n}], y /@ Range@n, {x, 0, 6}]];
NDSolve`ProcessEquations[sys]

In v9.0.1: ndsdtc: The time constraint of 1. seconds was exceeded trying to solve for derivatives, so the system will be treated as a system of differential-algebraic equations. You can use Method->{"EquationSimplification"->"Solve"} to have the system solved as ordinary differential equations.

In v11.2: ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

As one can see, though the warning is different, NDSolve gives up transforming the equation system with Solve method in both cases. (The default TimeConstraint is 1, actually. )

So, when setting Method -> {"EquationSimplification" -> "Residual"} / SolveDelayed -> True, you're turning to a cheaper transforming process for your equations.

At this point, you may be thinking that "OK, I'll always use SolveDelayed -> True from now on", but I'm afraid it's not a good idea, because the DAE solver of NDSolve is generally weaker than the ODE solver. (Here is an example.) In certain cases, one may still need to force NDSolve to solve equation with ODE solver, which may be troublesome. (Here is an example. )

Finally, notice there exist many issues related to EquationSimplification -> Residual and I've avoided talking about most of them in order to keep this answer clean. If you want to know more about the topic, search in this site.

Answer 2

When you choose the Method "Residual", Mathematica will convert an explicit equation like

x''[t]==y'[t]+z[t]

into the form

0==y'[t]+z[t]-x''[t]

and use the solutions numerically closest to 0, since rounding errors would cause contradictions between the various explicit definitions in a set of long and numerically expensive functions when Mathematica gives the error message with the suggestion to use that Method.