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

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more