# 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[7] 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[1]'[t], x[2]'[t], x[3]'[t], x[4]'[t], x[5]'[t],
x[6]'[t], x[7]'[t], x[8]'[t], x[9]'[t]} ==
{-k[1] x[1][t] x[3][t] + k[2] x[5][t] + k[3] x[5][t] -
k[7] x[1][t] x[7][t] + k[8] x[8][t],
-k[4] x[2][t] x[4][t] + k[5] x[6][t] + k[6] x[6][t] -
k[9] x[2][t] x[7][t] + k[10] x[9][t],
-k[1] x[1][t] x[3][t] + k[2] x[5][t] + k[6] x[6][t],
-k[4] x[2][t] x[4][t] + k[3] x[5][t] + k[5] x[6][t],
k[1] x[1][t] x[3][t] - k[2] x[5][t] - k[3] x[5][t],
k[4] x[2][t] x[4][t] - k[5] x[6][t] - k[6] x[6][t],
-k[7] x[1][t] x[7][t] - k[9] x[2][t] x[7][t] + k[8] x[8][t] + k[10] x[9][t],
k[7] x[1][t] x[7][t] - k[8] x[8][t],
k[9] x[2][t] x[7][t] - k[10] x[9][t]};

init = {T[1], T[2], T[3], 0, 0, 0, T[4], 0, 0};

• Please include your current NDSolve[ ]code Jun 18, 2015 at 16:09
• The Subscripts, they hurt my eyes! I recommend using, for instance, x[1][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]. Jun 18, 2015 at 16:17
• These two posts may be helpful: (69756) and (84471). Jun 18, 2015 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? Jun 18, 2015 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[1]'[t], x[2]'[t], x[3]'[t], x[4]'[t], x[5]'[t],
x[6]'[t], x[7]'[t], x[8]'[t],
x[9]'[t]} == {-k[1] x[1][t] x[3][t] + k[2] x[5][t] +
k[3] x[5][t] - k[7] x[1][t] x[7][t] + k[8](*[t]*) x[8][t],
-k[4] x[2][t] x[4][t] + k[5] x[6][t] + k[6] x[6][t] -
k[9] x[2][t] x[7][t] + k[10] x[9][t],
-k[1] x[1][t] x[3][t] + k[2](*[t]*) x[5][t] +
k[6] x[6][t], -k[4] x[2][t] x[4][t] + k[3] x[5][t] + k[5] x[6][t],
k[1] x[1][t] x[3][t] - k[2] x[5][t] - k[3] x[5][t],
k[4] x[2][t] x[4][t] - k[5] x[6][t] - k[6] x[6][t],
-k[7] x[1][t] x[7][t] - k[9] x[2][t] x[7][t] + k[8] x[8][t] +
k[10] x[9][t], k[7] x[1][t] x[7][t] - k[8] x[8][t],
k[9] x[2][t] x[7][t] - k[10] x[9][t]};

init = {T[1], T[2], T[3], 0, 0, 0, T[4], 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[0]] == init,
With[{df = Through[dvars[t]]},
WhenEvent[Norm[df] < ssthreshold, x[7][t] -> x[7][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[1]["Domain"] /. sol},
PlotLegends -> Automatic}


• This seems good. Thanks. Jun 19, 2015 at 9:19
• Is it possible to get the events time ? Jun 19, 2015 at 14:51
• @LifeWorks Sure. Something like Reap@NDSolve[{..., WhenEvent[Norm[df] < ssthreshold, Sow[t]; x[7][t] -> x[7][t] + 0.1]...},...]. There are a few examples in the WhenEvent documentation; just search for Reap or Sow. Jun 19, 2015 at 21:59
• @MichaelE2 I've tried reproducing this solution and run into an error. The domain for the x[1] 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? Jun 13, 2016 at 14:28
• @tarhawk I think it's a bug.. Jun 13, 2016 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[[1]]]}] == des[[2]] /. x[n_][t] -> x[n]],
Table[x[i], {i, Length[des[[1]]]}]]


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[3] -> ((k[2] + k[3]) x[5])/(k[1] x[1]),
x[4] -> (k[3] (k[5] + k[6]) x[5])/(k[4] k[6] x[2]),
x[6] -> (k[3] x[5])/k[6], x[8] -> (k[7] x[1] x[7])/k[8],
x[9] -> (k[9] x[2] x[7])/k[10]},
{x[1] -> 0, x[4] -> 0, x[5] -> 0, x[6] -> 0, x[8] -> 0, x[9] -> (k[9] x[2] x[7])/k[10]},
{x[1] -> 0, x[2] -> 0, x[5] -> 0, x[6] -> 0, x[8] -> 0, x[9] -> 0},
{x[2] -> 0, x[3] -> 0, x[5] -> 0, x[6] -> 0, x[8] -> (k[7] x[1] x[7])/k[8],
x[9] -> 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[[2]] /. x[n_][t] -> x[n] + e dx[n], e]

(* {dx[5] k[2] + dx[5] k[3] + dx[8] k[8] - dx[3] k[1] x[1] -
dx[7] k[7] x[1] - dx[1] k[1] x[3] - dx[1] k[7] x[7],
dx[6] k[5] + dx[6] k[6] + dx[9] k[10] - dx[4] k[4] x[2] -
dx[7] k[9] x[2] - dx[2] k[4] x[4] - dx[2] k[9] x[7],
dx[5] k[2] + dx[6] k[6] - dx[3] k[1] x[1] - dx[1] k[1] x[3],
dx[5] k[3] + dx[6] k[5] - dx[4] k[4] x[2] - dx[2] k[4] x[4],
-dx[5] k[2] - dx[5] k[3] + dx[3] k[1] x[1] + dx[1] k[1] x[3],
-dx[6] k[5] - dx[6] k[6] + dx[4] k[4] x[2] + dx[2] k[4] x[4],
dx[8] k[8] + dx[9] k[10] - dx[7] k[7] x[1] - dx[7] k[9] x[2] - dx[1] k[7] x[7] -
dx[2] k[9] x[7],
-dx[8] k[8] + dx[7] k[7] x[1] + dx[1] k[7] x[7],
-dx[9] k[10] + dx[7] k[9] x[2] + dx[2] k[9] x[7]} *)


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]][[2]] // Normal) /.
Thread[Array[x, 9] -> 0]]
(* {0, 0, 0, 0, 0, -k[2] - k[3], -k[5] - k[6], -k[8], -k[10]} *)


Hence, if any of k[2] + k[3], k[5] + k[6], k[8], or k[10] 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? Jun 19, 2015 at 16:46
• @LifeWorks This is something I should have done originally. I now have. Good luck. Jun 19, 2015 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[43];
sol = NDSolve[{R i[t] + L i'[t] == V, i[0] == 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[[1]])[t], {t, 0, 10}, PlotRange -> All]