How to use NDSolve to track equilibrium?

Question

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.

Answer 1

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["DifferentialEquations`InterpolatingFunctionAnatomy`"];

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.

Answer 2

You can use the EventLocator method of NDSolve.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

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!

Answer 3

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.

Answer 4

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}]