6
$\begingroup$

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ 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)}} $\endgroup$
    – Daniel Lichtblau
    Jun 30, 2020 at 21:07
  • $\begingroup$ @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. $\endgroup$
    – Daniel Farrell
    Jun 30, 2020 at 22:27

4 Answers 4

Reset to default
4
$\begingroup$

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[
 Thread[
  Eigenvalues[
   CoefficientArrays[rhs, {X[t], Y[t]}][[2]]
   ] < 0
  ],
 {A, k1, k2, k3}
 ]
(*  (A | k1) \[Element] Reals && k2 > -k1 + A k1 && k3 > 0  *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Getting the convergence criteria is very useful! Thanks $\endgroup$
    – Daniel Farrell
    Jun 30, 2020 at 15:50
  • 1
    $\begingroup$ @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. $\endgroup$
    – Michael E2
    Jun 30, 2020 at 16:19
4
$\begingroup$

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

stability conditions

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

Simplify[EcoStableQ[eq[[1]]]]
(* True *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$
    – Daniel Farrell
    Jun 30, 2020 at 15:51
  • 1
    $\begingroup$ Thanks, feel free to hit me up if you try it and have questions. $\endgroup$
    – Chris K
    Jun 30, 2020 at 16:00
2
$\begingroup$

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]

enter image description here

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

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

For your initially specified conditions, the additional requirement is

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

(* A k1 + k3 < k1 + k2 *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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 $\endgroup$
    – Daniel Farrell
    Jun 30, 2020 at 19:19
1
$\begingroup$

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}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.