11
$\begingroup$

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 ?

$\endgroup$
  • 8
    $\begingroup$ I wish it was more visible in the documentation, but I strongly recommend going through this tutorial. In particular this section. $\endgroup$ – chuy Jan 24 '17 at 17:52
  • 2
    $\begingroup$ @chuy It's the first tutorial in the doc-center search for EquationSimplification (second hit, after a "featured example"). It would be nice if it were easy to find everything, but everything is quite a lot of stuff, even just everything related to NDSolve. $\endgroup$ – Michael E2 Jan 24 '17 at 18:06
  • 1
    $\begingroup$ It's just my experience that the vast majority of users never know this tutorial exists. Even ones who use the documentation regularly. $\endgroup$ – chuy Jan 24 '17 at 18:11
  • $\begingroup$ @chuy Thanks for the link ! $\endgroup$ – Valacar Jan 25 '17 at 7:56
15
$\begingroup$

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

Mathematica graphics

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.

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ What do you mean by "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"? If you mean, the warning suggesting to use Method -> {"EquationSimplification" -> "Residual"} appears because NDSolve can't solve the system accurately with the default method, I'm afraid you're wrong. This warning is more related to pre-processing. $\endgroup$ – xzczd Oct 21 '17 at 3:48
  • 1
    $\begingroup$ @xzczd More specifically, it solves the second equation using the IDA DAE-solver, using Newton's method (I assume) to approximate the solution to the equation. I interpreted the "since..." to mean that the 2nd equation is solved approximately. I didn't think that Simon meant that the original system cannot be solved, but I can see now that yours is a possible interpretation. $\endgroup$ – Michael E2 Oct 21 '17 at 12:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.