# When using NDsolve, how to determine the positions of steady states?

I am currently trying to numerically solve a set of ordinary differential equations of chemical kinetics. However, I want to implement perturbations only when system reach steady states. For example increase one of the parameters or increase the concentrations of one chemical species.

For example, I want to solve the systems but want to increase, say x a small step each time when system reach steady state. Then how can I implement the trigger event (maybe with WhenEvent function)?

Resurrected code, without Subscript:

des = {x'[t], x'[t], x'[t], x'[t], x'[t],
x'[t], x'[t], x'[t], x'[t]} ==
{-k x[t] x[t] + k x[t] + k x[t] -
k x[t] x[t] + k x[t],
-k x[t] x[t] + k x[t] + k x[t] -
k x[t] x[t] + k x[t],
-k x[t] x[t] + k x[t] + k x[t],
-k x[t] x[t] + k x[t] + k x[t],
k x[t] x[t] - k x[t] - k x[t],
k x[t] x[t] - k x[t] - k x[t],
-k x[t] x[t] - k x[t] x[t] + k x[t] + k x[t],
k x[t] x[t] - k x[t],
k x[t] x[t] - k x[t]};

init = {T, T, T, 0, 0, 0, T, 0, 0};

• Please include your current NDSolve[ ]code – Dr. belisarius Jun 18 '15 at 16:09
• The Subscripts, they hurt my eyes! I recommend using, for instance, x[t] in place of Subscript[x,1][t], because Subscripts act weird in Mathematica. Also, your init: is this meant to be a set of initial conditions that you are feeding to NDSolve, because you should have == in place of =. Finally, please include values for your initial conditions Subscript[T, i]. – march Jun 18 '15 at 16:17
• These two posts may be helpful: (69756) and (84471). – march Jun 18 '15 at 16:22
• Hi, I restored your code, programmatically changing the subscripts. You can roll back if you like. (There's a "roll back" button in the edit history.) I guess it's obvious that the x's are your variables, but it's not real clear what the k's and T's are. Are they just real constants? – Michael E2 Jun 18 '15 at 22:39

[Update notice: NDSolve code below used to work without the Method option, but somewhere around V10.4, the option became necessary. Also see this.]

Based on the OP's original "set-up":

des = {x'[t], x'[t], x'[t], x'[t], x'[t],
x'[t], x'[t], x'[t],
x'[t]} == {-k x[t] x[t] + k x[t] +
k x[t] - k x[t] x[t] + k(*[t]*) x[t],
-k x[t] x[t] + k x[t] + k x[t] -
k x[t] x[t] + k x[t],
-k x[t] x[t] + k(*[t]*) x[t] +
k x[t], -k x[t] x[t] + k x[t] + k x[t],
k x[t] x[t] - k x[t] - k x[t],
k x[t] x[t] - k x[t] - k x[t],
-k x[t] x[t] - k x[t] x[t] + k x[t] +
k x[t], k x[t] x[t] - k x[t],
k x[t] x[t] - k x[t]};

init = {T, T, T, 0, 0, 0, T, 0, 0};


Use the Norm of the derivative (with WhenEvent) to determine when the system is close to a steady-state (same idea as belisarius, who answered first).

vars = Array[x, 9];

Block[{k, T, ssthreshold},
k[n_] := k[n] = (SeedRandom[n]; RandomReal[]);
T[n_] := T[n] = (SeedRandom[n + 10]; RandomReal[]);
ssthreshold = 1.*^-4;
(* Print[des]; *) (* to see the ODE *)
{sol} =
NDSolve[{des, Through[vars] == init,
With[{df = Through[dvars[t]]},
WhenEvent[Norm[df] < ssthreshold, x[t] -> x[t] + 0.1]]},
vars, {t, 0, 200}, MaxSteps -> 100000,
Method -> {"EquationSimplification" -> "Residual"}  (* needed as of V10.4 or so *)
]
];

Plot @@ {Through[vars[t]] /. sol, Flatten@{t, x["Domain"] /. sol},
PlotLegends -> Automatic} • This seems good. Thanks. – LifeWorks Jun 19 '15 at 9:19
• Is it possible to get the events time ? – LifeWorks Jun 19 '15 at 14:51
• @LifeWorks Sure. Something like Reap@NDSolve[{..., WhenEvent[Norm[df] < ssthreshold, Sow[t]; x[t] -> x[t] + 0.1]...},...]. There are a few examples in the WhenEvent documentation; just search for Reap or Sow. – Michael E2 Jun 19 '15 at 21:59
• @MichaelE2 I've tried reproducing this solution and run into an error. The domain for the x variable and all others is from {t,0.,0.} as written. Manually changing the the domain to {t, 0.,200.} throws an extrapolation error. Any chance you could take a look at the code? – tarhawk Jun 13 '16 at 14:28
• @tarhawk I think it's a bug.. – Michael E2 Jun 13 '16 at 15:36

Because the equations are nonlinear, there is no assurance that equilibrium solutions even exist for a given set of k and initial conditions. However, if they do exist, an alternative, and perhaps more informative, approach is to compute them directly:

xs = Solve[Thread[Table[0, {i, Length[des[]]}] == des[] /. x[n_][t] -> x[n]],
Table[x[i], {i, Length[des[]]}]]


Although Solve warns

Solve::svars: Equations may not give solutions for all "solve" variables. >>


playing a bit with Reduce strongly suggests that Solve gives all the equilibrium solutions here:

{* {{x -> ((k + k) x)/(k x),
x -> (k (k + k) x)/(k k x),
x -> (k x)/k, x -> (k x x)/k,
x -> (k x x)/k},
{x -> 0, x -> 0, x -> 0, x -> 0, x -> 0, x -> (k x x)/k},
{x -> 0, x -> 0, x -> 0, x -> 0, x -> 0, x -> 0},
{x -> 0, x -> 0, x -> 0, x -> 0, x -> (k x x)/k,
x -> 0}} *}


(For each of the four solutions, any x that do not appear can take any value.) With the equilibrium solutions known, one can linearize the equations des about them, Laplace transform the resulting linear equations, and Solve them to obtain the Eigensystems, completely determining the perturbed solutions in closed form.

In answer to the question posed below in a comment, the right side of des can be linearized about an equilibrium solution without difficulty. (The left side already is linear.)

desl = Coefficient[des[] /. x[n_][t] -> x[n] + e dx[n], e]

(* {dx k + dx k + dx k - dx k x -
dx k x - dx k x - dx k x,
dx k + dx k + dx k - dx k x -
dx k x - dx k x - dx k x,
dx k + dx k - dx k x - dx k x,
dx k + dx k - dx k x - dx k x,
-dx k - dx k + dx k x + dx k x,
-dx k - dx k + dx k x + dx k x,
dx k + dx k - dx k x - dx k x - dx k x -
dx k x,
-dx k + dx k x + dx k x,
-dx k + dx k x + dx k x} *)


At this point, equilibrium values of x are substituted (via Rules) into the last result and the Eigenvaluess computed. These Eigenvaluess are the growth rates of the perturbations dx. (Complex Eigenvaluess indicate oscillatory solutions.) For instance, if all equilibrium values are 0, then

Eigenvalues[(CoefficientArrays[desl, Array[dx, 9]][] // Normal) /.
(* {0, 0, 0, 0, 0, -k - k, -k - k, -k, -k} *)


Hence, if any of k + k, k + k, k, or k are negative, then perturbations grow exponentially. More generally, one would substitute one of the four equilibrium solutions xs derived earlier into the perturbed equations.

Another advantage of the procedure described here is this: In the event that an equilibrium admits exponentially growing perturbations, then attempts to solve for that equilibrium using NDSolve are likely to be unstable too.

• Thanks a lot. Do you mind expand a bit on the EigenSystem part? – LifeWorks Jun 19 '15 at 16:46
• @LifeWorks This is something I should have done originally. I now have. Good luck. – bbgodfrey Jun 19 '15 at 18:21

An example perturbing a simple RL circuit when it's about to reach the steady state.

R = 1;
L = 1;
V = 1;
SeedRandom;
sol = NDSolve[{R i[t] + L i'[t] == V, i == 0,
WhenEvent[{i'[t] ==  .01}, i[t] -> i[t] + RandomReal[{-.5, .5}]],
WhenEvent[{i'[t] == -.01}, i[t] -> i[t] + RandomReal[{-.5, .5}]]},
i, {t, 0, 10}]
Plot[(i /. sol[])[t], {t, 0, 10}, PlotRange -> All] 