7
$\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$ Commented 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$ Commented Jun 30, 2020 at 22:27

4 Answers 4

Reset to default
5
$\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$ Commented 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$ Commented 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$ Commented Jun 30, 2020 at 15:51
  • 1
    $\begingroup$ Thanks, feel free to hit me up if you try it and have questions. $\endgroup$ Commented Jun 30, 2020 at 16:00
3
$\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$ Commented Jun 30, 2020 at 19:19
2
$\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$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.