# How to use NDSolve to track equilibrium?

I am looking for an extension of NDSolve where integration runs until certain variables are settled at an equilibrium. Now I have a working solution in my sleeves but I would rather not show it yet, as I want answers to be unbiased and original, since I'm not convinced that my solution is optimal at all.

Consider the following putative example (note that EquilibriumNDSolve is an undefined hypothetical function):

eqns = {
Derivative[1][a][t] == -a[t] - 0.2 a[t]^2 + 2.1 b[t],
Derivative[1][b][t] == a[t] + 0.1 a[t]^2 - 1.1 b[t],
a[0] == 0.5,
b[0] == 0.5
};

steps = {};

sol = EquilibriumNDSolve[eqns, {a, b}, {t, 0, 1000}, a + b,
EquilibriumThreshold -> 10^-5,
EquilibriumStepMonitor :> AppendTo[steps, t]
];


This reads as: "numerically solve eqns for a and b variables while t goes from 0, until (a + b) is settled at an equilibrium. If no equilibrium is reached until t = 1000", terminate. The solution should be a set of {y.i -> InterpolatingFunction[...]} (just like in case of NDSolve), and the result should be something like this:

Plot[Evaluate[{a[t], b[t]} /. sol], {t, 0, Last[steps]},
PlotStyle -> AbsoluteThickness[2], ImageSize -> 400,
GridLines -> {steps, {}}, GridLinesStyle -> {Dashed, GrayLevel[.7]}]


Vertical gridlines indicate positions where an equilibrium test was performed. The syntax and usage should vaguely work like this:

EquilibriumNDSolve[eqns, {y1, y2, ...}, {t, t0, tmax}, z}] numerically solves the differential equation system eqns for all variables $y_i$, with independent variable t running from $t_0$ to a maximum of $t_{max}$. Iteration stops when variable z settles at an equilibrium or no equilibrium was found in the given range of t. Variable z can be a single variable, a list of variables, or any sensible combination of them, e.g. $a+b$.

EquilibriumThreshold should define a threshold value above which no equilibrium is registered. EquilibriumStepMonitor should work like StepMonitor: each time an equilibrium-test is performed, the rhs of EquilibriumStepMonitor :> func is evaluated as well.

-
 I'm curious to see what solution you have up your sleeve :) – JxB Jan 23 '12 at 6:13 @JxB: I had the very same solution everyone posted below, so no unexpected solutions have pop up. Mine has a fine tuned sampling method and some necessary additions to be able to package the function. If someone is interested, I can post the function in its entirety. – István Zachar Jan 27 '12 at 14:34

This approach finds equilibrium by checking that all derivatives up to the order of the differential equation are below a threshold. Following the template (defined below) suggested by the OP, here is an example for a damped harmonic oscillator:

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];

eqns1 = {a''[t] == Pi^2/2500 - (Pi^2*a[t])/2500 - 0.02*a'[t],
a[0] == 0., a'[0] == 0};

steps1 = {};
sol1 = equilibriumNDSolve[eqns1, {a}, {t, 0, 1000}, a,
equilibriumThreshold -> 1*^-5,
equilibriumStepMonitor :> AppendTo[steps1, t]];

end1 = InterpolatingFunctionDomain[a /. sol1[[1]]][[1, 2]];
Plot[a[t] /. sol1, {t, 0, end1}, PlotRange -> {0, 2},
Prolog -> {Thin, Dashed, Line[{{#, 0}, {#, 2}}] & /@ steps1[[1 ;; -1 ;; 2]]},
PlotStyle -> Thick, AxesLabel -> {t, a[t]}]


stops early (t=656.59), giving

The dashed vertical lines show the step monitor times. That is a 2nd order differential equation. The OP's example is a first order differential equation:

eqns2 = {Derivative[1][a][t] == -a[t] - 0.2 a[t]^2 + 2.1 b[t], Derivative[1][b][t] == a[t] + 0.1 a[t]^2 - 1.1 b[t], a[0] == 0.5, b[0] == 0.5};

steps2 = {};
sol2 = equilibriumNDSolve[eqns2, {a, b}, {t, 0, 1000}, a + b,
equilibriumThreshold -> 1*^-3,
equilibriumStepMonitor :> AppendTo[steps2, t]];
end2 = InterpolatingFunctionDomain[a /. sol2[[1]]][[1, 2]];
Plot[Evaluate[{a[t], b[t]} /. sol2], {t, 0, end2}, PlotRange -> Automatic,
Prolog -> {Thin, Dashed, Line[{{#, 0}, {#, 1100}}] & /@ steps2},
PlotStyle -> Thick, AxesLabel -> {t, Row[{a[t], ", ", b[t]}]}]


uses a threshold of 1*^-3, and stops at t=465.234:

Here is the definition of equilibriumNDSolve[]:

Clear[equilibriumNDSolve];
Options[equilibriumNDSolve] = {equilibriumThreshold :> 1*^-5, equilibriumStepMonitor -> None};
equilibriumNDSolve[eqns_, vars_, {t_, start_, finish_}, equilibriumexpr_, opts : OptionsPattern[]] :=
Module[{threshold, order},
threshold = OptionValue[equilibriumThreshold];
order = Max[Cases[eqns, Derivative[n_][_][_] :> n, Infinity]];
NDSolve[eqns, vars, {t, start, finish}, Method -> {"EventLocator",
"Event" -> And @@ ((Distribute@Abs[Through[
Distribute[Derivative[#][equilibriumexpr]][t], Plus]] <threshold) & /@    Range[order])},
StepMonitor :> OptionValue[equilibriumStepMonitor]]]


The key part is the "EventLocator" method of NDSolve, as pointed out by Sjoerd and Szabolcs.

The function expects the stopping criterion (equilibriumexpr) to involve at most the addition of the NDSolve variables (more complicated expressions do not work as-is). The transformation of equilibriumexpr into an Event is not clean (i.e., not easy to follow), and may not be robust, but it works for the two cases above.

-
 Nice solution checking the derivatives and dealing with DE of any order. – István Zachar Jan 27 '12 at 14:37 (+1), but one should point out that higher-order equations like the one in this example can be reduced to a system of first-order equations, so that the other methods based on systems of first-order equation also apply to this example and have more generality. – Jens May 8 '12 at 3:11 @Jens feel free to edit my answer to address your comment. – JxB May 8 '12 at 3:41

My variant of Szabolcs code. It doesn't need an extra package:

sol = First[
NDSolve[eqns, {a, b}, {t, 0, 1000}, Method -> {"EventLocator",
"Event" -> Abs[a'[t]] +Abs[b'[t]] < 10^-5,
"EventAction" :> Throw[end = t, "StopIntegration"]}]];

Plot[Evaluate[{a[t], b[t]} /. sol], {t, 0, end}]


As you can see it makes use of the "EventAction" option which stops the integration and also returns the end point in the end variable.

-
 Mine doesn't need an extra package either as you can extract the domain apparently by passing "Domain" to the interpolating function. But I can't find where this functionality is documented, so I didn't feel so secure about it. – Szabolcs Jan 22 '12 at 15:58

You can use the EventLocator method of NDSolve.

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];

eqns = {Derivative[1][a][t] == -a[t] - 0.2 a[t]^2 + 2.1 b[t],
Derivative[1][b][t] == a[t] + 0.1 a[t]^2 - 1.1 b[t],
a[0] == 0.5, b[0] == 0.5};

sol = First@
NDSolve[eqns, {a, b}, {t, 0, 1000},
Method -> {"EventLocator",
"Event" -> Abs[a'[t]] + Abs[b'[t]] > 0.1}]

end = InterpolatingFunctionDomain[a /. sol][[1, 2]]

Plot[{a[t], b[t]} /. sol // Evaluate, {t, 0, end}]


I used the condition Abs[a'[t]] + Abs[b'[t]] <= 0.1 to stop the integration.

I am not really familiar with EventLocator, so I had to check the documentation, and I might have missed something. There might be cleaner ways to do this.

At least in version 8, you can also get the interpolating function's domain using a["Domain"] /. sol, but I don't know where this is documented. I'd appreciate a pointer!

-
Bummer. I was working along the same lines. BTW Shouldn't the Abs[a'[t]] + Abs[b'[t]] > 0.1 in your code be <0.1? – Sjoerd C. de Vries Jan 22 '12 at 15:47
@Sjoerd Interesting point about the inequality. First I assumed the event must be a numerical quantity, and that the integration would stop when it crosses zero. So I actually used Boole[Abs[a'[t]] + Abs[b'[t]] > 0.1]. The I removed the Boole and noted that it still works, so I posted this version. The interesting thing is that it works no matter if we use >, < or ==. I think it just detects a change in truth value in this case. – Szabolcs Jan 22 '12 at 16:06
The documentation says "Event functions are expected to be real- or Boolean-valued, so if there is a change, there must be an event in the step interval" – Sjoerd C. de Vries Jan 22 '12 at 16:20

You can explicitly find the stable values, then integrate the ODE system until they are nearly attained.

sys = {Derivative[1][a][t] == -a[t] - 0.2 a[t]^2 + 2.1 b[t],
Derivative[1][b][t] == a[t] + 0.1 a[t]^2 - 1.1 b[t]};
inits = {a[0] == 0.5, b[0] == 0.5};
vars = {a[t], b[t]};

stopvals =
vars /. NSolve[Join[sys, Thread[D[vars, t] == 0]],
Join[vars, D[vars, t]]]

Out[279]= {{100., 1000.}, {0., 0.}}

thresh = .01;

endtime =
Reap[soln =
First@NDSolve[Join[sys, inits], vars, {t, 0, 10^4},
Method -> {"EventLocator",
"Event" ->
Abs[a[t] - stopvals[[1, 1]]] +
Abs[b[t] - stopvals[[1, 2]]] <= thresh,
"EventAction" :> Sow[t]}]][[2, 1, 1]]

Out[301]= 525.745

Plot[Evaluate[vars /. soln], {t, 0, endtime}]


-
 How come I did not think of this? Though there is the open question: how much extra time is required for running NSolve? My calculations show that the NSolve and NDSolve lines use the same magnitude of time, meaning that this code practically doubles the time of finding equilibrium (so it might be slower for more complicated situations). – István Zachar Jun 21 '12 at 23:03 @Istvan Zachar I would expect that NSolve won't be a bottleneck unless the equations are nonlinear in a[t] et al. – Daniel Lichtblau Jun 21 '12 at 23:33