3
$\begingroup$

I am trying to to find steady states to a system and for that I tried the following code. This is a system of 5 ODEs (b,be,ei,eu,n are the variables).

Assuming[{x1 > 0, x2 > 0, x3 > 0, x4 > 0, x5 > 0, x6 > 0, x7 > 0, 
  x8 > 0, x9 > 0, x10 > 0, x11 > 0, y > 0}, 
 Simplify[Solve[
   b[t]*(1 - (b[t] + be[t])) - x1*b[t]*eu[t] - 
      x2*b[t]*n[t]/(n[t] + b[t]) - x3*b[t] + 
      ei[t]*be[t]*(1 - (b[t] + be[t])) == 
     0 && (1 - ei[t])*be[t]*(1 - (b[t] + be[t])) + x1*b[t]*eu[t] - 
      x4*be[t]*n[t]/(n[t] + be[t]) == 0 && 
    x5*b[t]*eu[t] - x6*ei[t] + 
      x7*(1 - ei[t])*be[t]*(1 - (b[t] + be[t])) == 0 && 
    x8 - x5*b[t]*eu[t] - x8*eu[t] - 
      x7*(1 - ei[t])*be[t]*(1 - (b[t] + be[t])) == 
     0 && (x9 - x10*n[t])*(b[t] + be[t]) - x11*n[t] == 0 && 
    0 <= b[t] <= 1 && 0 <= be[t] <= 1 && 0 <= ei[t] <= 1 && 
    0 <= eu[t] <= 1 && 0 <= n[t] <= y, {b[t], be[t], ei[t], eu[t], 
    n[t]}, Reals]]]

But, Mathematica would run for a very long time without giving a solution (tried the above using Reduce and NSolve as well).

So, I thought of solving the equations by hand.
By setting (x9 - x10*n[t])*(b[t] + be[t]) - x11*n[t] == 0 I can find n and then the plan is to substitute this to (b + be)*(1 - (b + be)) - n*((x2*b/(n + b)) + (x4*be/(n + be))) - x3*b == 0 and solve this for b, so I would get a solution for b in terms of be.

For this I tried out,

n = x9*(b + be)/(x11 + x10*(b + be));   
Assuming[{x9 > 0, x11 > 0,x10 > 0, x3 > 0, x2 > 0, x4 > 0}, 
     Simplify[Solve[(b + be)*(1 - (b + be)) - 
          n*((x2*b/(n + b)) + (x4*be/(n + be))) - x3*b == 0 && 
        0 <= b <= 1 && 0 <= be <= 1, {b}, Reals]]]

and

n = x9*(b + be)/(x11 + x10*(b + be));
subdomain = 0 <= b <= 1 && 0 <= be <= 1;
NSolve[subdomain && (b + be)*(1 - (b + be)) - 
    n*((x2*b/(n + b)) + (x4*be/(n + be))) - x3*b == 0, {b}, Reals]

But in both ways, Mathematica would run for a very long time and still wouldn't produce a result.

Can someone please tell me how to find an analytical solution to this system.


The solution outputs are based on the answer by @bbgodfrey

How can I read this solution obtained from Reduce.
So, is one solution is (b == 0 && ei == (be (-1 + b + be) x7)/(-x6 + be (-1 + b + be) x7) && eu == (be (-1 + b + ei - b ei) x7 + be^2 (x7 - ei x7) + x8)/( b x5 + x8)) with n == ((b + be) x9)/(b x10 + be x10 + x11).
Then the other solution is

(((be == 0 && b > 0 && 
      ei x6 (b x5 + x8) == b x5 x8) || (0 < be < 
       1 && ((ei == (((-1 + be) be x7 + 
              b (-x5 + be x7)) x8)/(-(x6 - (-1 + be) be x7) x8 + 
            b (-x5 x6 + be x7 x8)) && (x5 == (be (-1 + b + be) x7)/b ||
             x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || 
            x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8)) && (b + be >= 
             1 || b > 0) && (b + be < 1 || 
            x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || b + be > 1 ||
             x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8))) || (x5 == (
           be (-1 + b + be) x7)/b && 
          x6 == (be (-1 + b + be) x7 x8)/(
           b x5 + x8) && ((b > 0 && b + be < 1) || 
            b + be > 1)))) || (be == 1 && 
      b > 0 && (ei == (b (x5 - x7) x8)/(b x5 x6 + x6 x8 - b x7 x8) || 
        x6 (b x5 + x8) == b x7 x8) && (x5 == x7 || 
        x6 (b x5 + x8) != b x7 x8))) && 
  eu == (be (-1 + b + be) x7 + ei (x6 - be (-1 + b + be) x7))/(b x5)) with same `n`?

I cannot understand how to read this part of the solution

((be == 0 && b > 0 && 
    ei x6 (b x5 + x8) == b x5 x8) || (0 < be < 
     1 && ((ei == (((-1 + be) be x7 + 
            b (-x5 + be x7)) x8)/(-(x6 - (-1 + be) be x7) x8 + 
          b (-x5 x6 + be x7 x8)) && (x5 == (be (-1 + b + be) x7)/b || 
          x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || 
          x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8)) && (b + be >= 1 ||
           b > 0) && (b + be < 1 || 
          x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || b + be > 1 || 
          x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8))) || (x5 == (
         be (-1 + b + be) x7)/b && 
        x6 == (be (-1 + b + be) x7 x8)/(
         b x5 + x8) && ((b > 0 && b + be < 1) || 
          b + be > 1)))) || (be == 1 && 
    b > 0 && (ei == (b (x5 - x7) x8)/(b x5 x6 + x6 x8 - b x7 x8) || 
      x6 (b x5 + x8) == b x7 x8) && (x5 == x7 || 
      x6 (b x5 + x8) != b x7 x8)))

In this part does it mean that,

eu == (be (-1 + b + be) x7 + ei (x6 - be (-1 + b + be) x7))/(b x5))) &&
  n == ((b + be) x9)/(b x10 + be x10 + x11)

and the other possible solutions can be
1) (be == 0 && b > 0 && ei x6 (b x5 + x8) == b x5 x8) or
2)

(be == 1 && 
  b > 0 && (ei == (b (x5 - x7) x8)/(b x5 x6 + x6 x8 - b x7 x8) || 
    x6 (b x5 + x8) == b x7 x8) && (x5 == x7 || 
    x6 (b x5 + x8) != b x7 x8)) or  

3)

(0 < be < 
   1 && ((ei == (((-1 + be) be x7 + 
          b (-x5 + be x7)) x8)/(-(x6 - (-1 + be) be x7) x8 + 
        b (-x5 x6 + be x7 x8)) && (x5 == (be (-1 + b + be) x7)/b || 
        x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || 
        x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8)) && (b + be >= 1 || 
        b > 0) && (b + be < 1 || 
        x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || b + be > 1 || 
        x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8))) || (x5 == (
       be (-1 + b + be) x7)/b && 
      x6 == (be (-1 + b + be) x7 x8)/(
       b x5 + x8) && ((b > 0 && b + be < 1) || b + be > 1))))  

Inside this 3rd solution there are several other possibilities that can happen,is it? But, how do I read off those solutions. Specially I don't understand the part of

(x5 == (be (-1 + b + be) x7)/b || 
  x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || 
  x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8))  
$\endgroup$
  • $\begingroup$ Reduce almost certain can obtain a solution. However, because this is a fifteenth order system, the computation will be very slow, and the results is unlikely to be simple. $\endgroup$ – bbgodfrey Nov 13 '17 at 19:14
4
$\begingroup$

As I noted in an earlier comment, Reduce almost certainly can obtain a solution. However, because this is a fourteenth order system, the computation will be very slow, and the results is unlikely to be simple. I have been running the general problem for about seven hours so far.

Here is a partial solution. The five equations with [t] removed (unnecessary) and denominators cleared are,

eq = 
  b^3 + (-1 + be) be ei n + b^2 (-1 + be + be ei + n + eu x1 + x3) + 
    b (be^2 ei + be (ei (-1 + n) + n) + n (-1 + eu x1 + x2 + x3)) == 0 &&
  be^3 (-1 + ei) + be^2 (-1 + ei) (-1 + b + n) + b be eu x1 + b eu n x1 + 
    be n (1 + b (-1 + ei) - ei - x4) == 0 && 
  b eu x5 + b be (-1 + ei) x7 + (-1 + be) be ei x7 == ei x6 + (-1 + be) be x7 && 
  be (-1 + ei) x7 + be^2 (x7 - ei x7) + x8 == b eu x5 + b be (-1 + ei) x7 + eu x8 &&
  (b + be) (-n x10 + x9) == n x11

Now, solve for the three values, {ei, eu, n},

seieun = Simplify[Reduce[eq[[3 ;; 5]], {ei, eu, n}, Reals], 
    {x1 > 0, x2 > 0, x3 > 0, x4 > 0, x5 > 0, x6 > 0, 
    x7 > 0, x8 > 0, x9 > 0, x10 > 0, x11 > 0, y > 0, 0 <= b <= 1, 0 <= be <= 1}]

which yields after several minutes,

(* ((b == 0 && ei == (be (-1 + b + be) x7)/(-x6 + be (-1 + b + be) x7) &&
   eu == (be (-1 + b + ei - b ei) x7 + be^2 (x7 - ei x7) + x8)/(b x5 + x8)) || 
  (((be == 0 && b > 0 && ei x6 (b x5 + x8) == b x5 x8) || (0 < be < 1 && 
  ((ei == (((-1 + be) be x7 + b (-x5 + be x7)) x8)/(-x6 x8 + (-1 + be) be x7 x8 + 
  b (-x5 x6 + be x7 x8)) && (x5 == (be (-1 + b + be) x7)/b || 
 x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8)) &&
  (b + be >= 1 || b > 0) && (b + be < 1 || x6 < (be (-1 + b + be) x7 x8)/(b x5 + x8) || 
  b + be > 1 || x6 > (be (-1 + b + be) x7 x8)/(b x5 + x8))) || 
  (x5 == (be (-1 + b + be) x7)/b && x6 == (be (-1 + b + be) x7 x8)/(b x5 + x8) &&
  ((b > 0 && b + be < 1) || b + be > 1)))) || (be == 1 && b > 0 && 
  (ei == (b (x5 - x7) x8)/(b x5 x6 + x6 x8 - b x7 x8) || x6 (b x5 + x8) == b x7 x8) &&
  (x5 == x7 || x6 (b x5 + x8) != b x7 x8))) && 
  eu == (be (-1 + b + be) x7 + ei (x6 - be (-1 + b + be) x7))/(b x5))) && 
  n == ((b + be) x9)/(b x10 + be x10 + x11) *)

Examining this maze of options indicates that

n == ((b + be) x9)/(b x10 + be x10 + x11)

in all cases, as expected. There are two solutions for eu:

eu == (be (-1 + b + ei - b ei) x7 + be^2 (x7 - ei x7) + x8)/(b x5 + x8)

with corresponding ei == (be (-1 + b + be) x7)/(-x6 + be (-1 + b + be) x7) and b == 0 Vanishing b appears to over-constrain the remaining two equations, although one cannot be certain without substituting into them the expressions for {ei, eu, n} and attempting to Reduce the resulting equations. The second solution is

eu == (be (-1 + b + be) x7 + ei (x6 - be (-1 + b + be) x7))/(b x5)

with several solutions for ei and multiple constraints on the constants. To proceed further, it is necessary to substitute each of these alternative solutions into the first two equations and attempt to solve them.

Addendum - Quick Solution

After thirty hours, the brute-force computation mentioned at the beginning of this answer still is running. However, because the OP is interested only in solutions for arbitrary values of the constants (see comments below), assigning random values to the twelve constants and solving the resulting equations produces the correct answer in seconds. For instance,

Thread[Rule[{x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, y}, RandomInteger[100, 12]/10]]
Simplify[Reduce[eq /. %, {b, be, ei, eu, n}, Reals], {0 <= b <= 1, 0 <= be <= 1}]
(* {x1 -> 49/5, x2 -> 5/2, x3 -> 27/5, x4 -> 73/10, x5 -> 6, x6 -> 8, 
    x7 -> 29/5, x8 -> 83/10, x9 -> 7/10, x10 -> 14/5, x11 -> 12/5, y -> 33/10} *)
(* b == 0 && be == 0 && ei == 0 && eu == 1 && n == 0 *)

Repeating this computation several times with different random numbers produces the same results, as expected. To emphasize, this result is not merely a correct answer, it is the only generally correct answer.

For completeness, this same result can be obtained from the first two lines of the block of solutions named seieun above:

nrule = seieun // Last // ToRules
eq1 = Map[Numerator, Together[eq[[1 ;; 2]] /. %], {2}] // Simplify;
(* {n -> ((b + be) x9)/(b x10 + be x10 + x11)} *)

erule = seieun // First // First // ToRules
Map[Numerator, Together[eq1 //. %], {2}] // Simplify;
Simplify[Reduce[%, be, Reals], {x1 > 0, x2 > 0, x3 > 0, x4 > 0, x5 > 0, x6 > 0, x7 > 0, 
    x8 > 0, x9 > 0, x10 > 0, x11 > 0, y > 0, 0 <= b <= 1, 0 <= be <= 1}] // ToRules
(* {b -> 0, ei -> (be (-1 + b + be) x7)/(-x6 + be (-1 + b + be) x7), 
    eu -> (be (-1 + b + ei - b ei) x7 + be^2 (x7 - ei x7) + x8)/(b x5 + x8)} *)
(* be -> 0 *)

Flatten@{nrule, erule[[2 ;; 3]]} //. {b -> 0, be -> 0}
(* {n -> 0, ei -> 0, eu -> 1} *)

as desired. Finally, to confirm after all these computations that this solution is correct,

(List @@ eq) //. {n -> 0, ei -> 0, eu -> 1, b -> 0, be -> 0}
(* {True, True, True, True, True} *)
$\endgroup$
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Kuba Nov 16 '17 at 22:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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