# Solve system of ODEs for steady state symbolically

I can successfully solve a linear system of ODEs describing chemical reactions by first using DSolveValue to find the general solution and then taking the limit to infinity to find the steady state solution.

ClearAll["Global*"];
eqns = {
X'[t] == P - (k1 + k2) X[t] + A k1 X[t],
Y'[t] == k1 X[t] - A k1 X[t] - k3 Y[t],
X[0] == 0,
Y[0] == 0
}
s = Assuming[
k0 > 0 && k1 > 0 && k2 > 0 && k3 > 0 && P > 0 && A >= 0 && A < 1,
Limit[
DSolveValue[eqns,{X[t], Y[t]}, t], t -> \[Infinity]
]
]


This has the solution $$\left\{\frac{P}{-A \text{k1}+\text{k1}+\text{k2}},\frac{(A-1) \text{k1} P}{\text{k3} ((A-1) \text{k1}-\text{k2})}\right\}$$

Does mathematica provide a better/builtin way of doing this? For example by using eigenvalue method to find the steady solution directly?

• Could just set the derivatives to zero and solve for {X,Y}: In[4]:= equil = {P - (k1 + k2) xx + A k1 xx, k1 xx - A k1 xx - k3 yy}; Solve[equil == 0, {xx, yy}] Out[5]= {{xx -> -(P/(-k1 + A k1 - k2)), yy -> ((-k1 + A k1) P)/((-k1 + A k1 - k2) k3)}} Jun 30, 2020 at 21:07
• @DanielLichtblau Yes that’s very clean approach, because this defines the steady state solution I can just treat them as standard equations that are solved for zero. Jun 30, 2020 at 22:27

This will give the equilibrium, without checking convergence as t -> Infinity:

rhs = {X'[t], Y'[t]} /.
First@Solve[eqns, {X'[t], Y'[t]}, {X[0], Y[0]}];

Solve[rhs == 0, {X[t], Y[t]}]
(*
{{X[t] -> -(P/(-k1 + A k1 - k2)),
Y[t] -> ((-k1 + A k1) P)/((-k1 + A k1 - k2) k3)}}
*)


Here's way to get the convergence criteria:

Reduce[
Eigenvalues[
CoefficientArrays[rhs, {X[t], Y[t]}][[2]]
] < 0
],
{A, k1, k2, k3}
]
(*  (A | k1) \[Element] Reals && k2 > -k1 + A k1 && k3 > 0  *)

• Getting the convergence criteria is very useful! Thanks Jun 30, 2020 at 15:50
• @DanielFarrell You're welcome, and thank you for the accept. It just occurred to me I was assuming everything was real (incl. the eigenvalues). If complex eigenvalues are possible, you could take the real part of the eigenvalues: Reduce[ Flatten@{ {A, k1, k2, k3} \[Element] Reals, Thread[ 0 > Re@Eigenvalues[ CoefficientArrays[rhs, {X[t], Y[t]}][[2]] ]] }, {A, k1, k2, k3}] -- it comes out the same in this case. Jun 30, 2020 at 16:19

Although this is not an ecological model, my EcoEvo package could be useful.

(* load package *)
<< EcoEvo;

(* define model *)
SetModel[{
Aux[X] -> {Equation :> P - (k1 + k2) X[t] + A k1 X[t]},
Aux[Y] -> {Equation :> k1 X[t] - A k1 X[t] - k3 Y[t]},
Assumptions -> {k0 > 0, k1 > 0, k2 > 0, k3 > 0, 0 <= A < 1}
}]

(* find equilibria *)
eq = SolveEcoEq[]
(* {{X -> -(P/(-k1 + A k1 - k2)),
Y -> ((-k1 + A k1) P)/((-k1 + A k1 - k2) k3)}} *)

(* check stability *)
EcoStableQ[eq[[1]]]


Based on your assumptions, you will find that the equilibrium is always stable:

Simplify[EcoStableQ[eq[[1]]]]
(* True *)

• An interesting package. I’m just getting started with Mathematica so I picked the other answer. But this looks useful for when I actually understand what I’m doing! Jun 30, 2020 at 15:51
• Thanks, feel free to hit me up if you try it and have questions. Jun 30, 2020 at 16:00

Your initial approach also provides the conditions for convergence

\$Version

(* "12.1.1 for Mac OS X x86 (64-bit) (June 19, 2020)" *)

Clear["Global*"]

eqns = {X'[t] == P - (k1 + k2) X[t] + A k1 X[t],
Y'[t] == k1 X[t] - A k1 X[t] - k3 Y[t], X[0] == 0, Y[0] == 0};

sol = DSolve[eqns, {X, Y}, t][[1]];


Verifying the solution

eqns /. sol // Simplify

(* {True, True, True, True} *)


Using Limit to find the steady state also provides the conditions for convergence

ss = Limit[{X[t], Y[t]} /. sol, t -> Infinity]


cond = (And @@ (Last /@ ss)) // Simplify

(* P ∈ Reals && A k1 + k3 < k1 + k2 && k1 + k2 > A k1 && k3 > 0 *)


Simplify[cond,
k0 > 0 && k1 > 0 && k2 > 0 && k3 > 0 && P > 0 && A >= 0 && A < 1]

(* A k1 + k3 < k1 + k2 *)

• That’s true, my original code seems to take much longer. It’s very useful to see how you have manipulated the result here. Thanks Jun 30, 2020 at 19:19

Assuming that the system is stable

ClearAll["Global*"];
eqns = {X'[t] - P + (k1 + k2) X[t] + A k1 X[t], Y'[t] - k1 X[t] + A k1 X[t] + k3 Y[t]};
teqns = LaplaceTransform[eqns, t, s]
sols = Solve[teqns == 0, {LaplaceTransform[X[t], t, s], LaplaceTransform[Y[t], t, s]}][[1]]
Limit[s {LaplaceTransform[X[t], t, s], LaplaceTransform[Y[t], t, s]} /. sols, {s -> 0}]