Skip to main content

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


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




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.


Comments

Post a Comment

Popular posts from this blog

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1....
Post a Comment
Read more