Notebook source code: examples/stability/stability.ipynb
Run the notebook yourself on binder Binder badge

Stability analysis of a bacterial growth model through computer algebra

A. Yabo, M. Safey El Din, J.-B. Caillau and J.-L. Gouzé

This notebook is the log of the run of the original Maple worksheet. The solution is based on the RAGlib Maple package developed by the second author.

Thumbnail

Auxiliary functions

> SmallMidPoint:=proc(a,b)
> local d, n, delta, r, i, it;
> if a =b then return a; fi;
> if b<a then error "bad interval for constructing rationals"; fi;
> it:=2:
> while true do
> r:=convert((a+b)/2,confrac,'m', it);
> for i from 1 to nops(m) do
>    if m[i]> a and m[i]<b then
>       return m[i];
>    fi;
> od:
> it:=2*it:
> od:
> end;
SmallMidPoint := proc(a, b)
local d, n, delta, r, i, it;
    if a = b then return a end if;
    if b < a then error "bad interval for constructing rationals" end if;
    it := 2;
    do
        r := convert(1/2*a + 1/2*b, confrac, 'm', it);
        for i to nops(m) do
            if a < m[i] and m[i] < b then return m[i] end if
        end do;
        it := 2*it
    end do
end proc


> ConstructSmallRational:=proc(a,b)
> local d, n, delta;
> if a =b then return a; fi;
> if b<a then error "bad interval for construction rationals"; fi;
> return SmallMidPoint(a, b);
> end;
ConstructSmallRational := proc(a, b)
local d, n, delta;
    if a = b then return a end if;
    if b < a then error "bad interval for construction rationals" end if;
    return SmallMidPoint(a, b)
end proc

> find_rationals:=proc(listrat)
> local rationals,expo,apps,boo,den,newapps,i, newrats:
>  if nops(listrat)=0 then return [0] fi;
>  newrats:=[];
>  for i from 2 to nops(listrat) do
>     newrats:= [op(newrats),ConstructSmallRational(listrat[i-1], listrat[i])];
>  od;
>  newrats:=[floor(listrat[1]-1),op(newrats), ceil(listrat[nops(listrat)]+1)];
> return newrats;
> end:

Original problem. Equilibrium conditions (i.e. dsdt = dVdt = dpdt = 0)

> dsdt = V*k2*s*(u - 1)/(K2 + s) - d*(s - 1);
                              V k2 s (u - 1)
                       dsdt = -------------- - d (s - 1)
                                  K2 + s

> dVdt = -d + p*u/(K + p);
                                            p u
                               dVdt = -d + -----
                                           K + p

> dpdt = -d*(p + 1) + k1*p*(u - 1)/(K1 + p) - k2*s*(u - 1)/(K2 + s);
                                    k1 p (u - 1)   k2 s (u - 1)
                dpdt = -d (p + 1) + ------------ - ------------
                                       K1 + p         K2 + s

Positivity of the variables / parameters (i.e. ineq1 > 0, ineq2 > 0)

> inequalities1 = [k1, K1, u, 1-u, k2, K2, d, K];
                inequalities1 = [k1, K1, u, 1 - u, k2, K2, d, K]

> inequalities2 = [s,V,p];
                           inequalities2 = [s, V, p]

Positivity of the (strict) negation of Routh-Hurwitz condition (i.e. pol > 0)

The condition to be checked is

\[\phi_1 - {\left(2 \, D + \phi_2 + \phi_3 + \phi_4\right)} {\left(D^{2} + D \phi_2 + D \phi_3 + D \phi_4 + \phi_2 \phi_4 + \phi_3 \phi_4 \right)} \gt 0 \tag{POL}\]

where

\[\begin{split}\phi_1(\theta) = \dfrac{\mathcal{V}^* k_{2}^{2} s^* {\left(u^* - 1\right)}^{2} u^* {\left(\frac{p^*}{K + p^*} - 1\right)} {\left(\frac{s^*}{K_{2} + s^*} - 1\right)}}{{\left(K + p^*\right)} {\left(K_{2} + s^*\right)}^{2} }\,,\\ \phi_2(\theta) = -\dfrac{{\left(p^* + 1\right)} u^* {\left(\frac{p^*}{K + p^*} - 1\right)}}{K + p^*}\,,\\ \phi_3(\theta) = \dfrac{k_{1} {\left(u^* - 1\right)} {\left(\frac{p^*}{K_{1} + p^*} - 1\right)}}{K_{1} + p^*}\,,\\ \phi_4(\theta) = \dfrac{\mathcal{V}^* k_{2} {\left(u^* - 1\right)} {\left(\frac{s^*}{K_{2} + s^*} - 1\right)}}{{\left(K_{2} + s^*\right)} } \cdot\end{split}\]

It is the strengthening of the negation of the Routh-Hurwitz condition of order \(3\). The negation of this condition actually is a large inequality, while the method checks emptyness / non-emptyness of each connected component of the open set defined by (POL).

> pol := V*k2*s*(k2*(u - 1)/(K2 + s) - k2*s*(u - 1)/(K2 + s)^2)*(u - 1)*(u/(K + p) - p*u/(K + p)^2)/(K2 + s) + (V*k2*(u - 1)/(K2 + s) - V*k2*s*(u - 1)/(K2 + s)^2 - d + k1*(u - 1)/(K1 + p) - k1*p*(u - 1)/(K1 + p)^2 - (p + 1)*u/(K + p) - p*u/(K + p) + (p + 1)*p*u/(K + p)^2)*(V*k2*(u - 1)/(K2 + s) - V*k2*s*(u - 1)/(K2 + s)^2 - d)*(k1*(u - 1)/(K1 + p) - k1*p*(u - 1)/(K1 + p)^2 - (p + 1)*u/(K + p) - p*u/(K + p) + (p + 1)*p*u/(K + p)^2);
              /k2 (u - 1)   k2 s (u - 1)\         /  u       p u   \
       V k2 s |---------- - ------------| (u - 1) |----- - --------|
              |  K2 + s              2  |         |K + p          2|
              \              (K2 + s)   /         \        (K + p) /   /
pol := ------------------------------------------------------------- + |
                                  K2 + s                               |
                                                                       \

    V k2 (u - 1)   V k2 s (u - 1)       k1 (u - 1)   k1 p (u - 1)   (p + 1) u
    ------------ - -------------- - d + ---------- - ------------ - ---------
       K2 + s                2            K1 + p              2       K + p
                     (K2 + s)                         (K1 + p)

        p u    (p + 1) p u\ /V k2 (u - 1)   V k2 s (u - 1)    \
     - ----- + -----------| |------------ - -------------- - d|
       K + p           2  | |   K2 + s                2       |
                (K + p)   / \                 (K2 + s)        /

    /k1 (u - 1)   k1 p (u - 1)   (p + 1) u    p u    (p + 1) p u\
    |---------- - ------------ - --------- - ----- + -----------|
    |  K1 + p              2       K + p     K + p           2  |
    \              (K1 + p)                           (K + p)   /

The latter means that, if pol>0, the equilibrium is unstable (we want to proove this never occurs, i.e. the problem has no solution).

Re-parametrized problem

We replaced K1 = K1p - p and K2 = K2p - s and we replaced the equilibrium conditions in every vector. Positivity of the variables / parameters (i.e. newineq1 > 0, newineq2 > 0).

> newineq1 := [-K2p*d*p - k2*s*u - K2p*d + k2*s, K1p - p, u, -u + 1, k2, K2p - s, d, -d + u];
newineq1 := [

    -K2p d p - k2 s u - K2p d + k2 s, K1p - p, u, 1 - u, k2, K2p - s, d, -d + u

    ]

> newineq2 := [s, -s + 1, p]:

Positivity of the (numerator of the) inverse Routh-Hurwitz condition (i.e. numpol > 0)

The latter because we already saw the denominator is always positive.

> numpol := -K1p^2*K2p^2*d^4*p^2*s^3 + 2*K1p^2*K2p^2*d^3*p^2*s^3*u + 2*K1p*K2p^2*d^3*p^3*s^3*u + K1p^2*K2p*d^3*p^2*s^4*u - K1p^2*K2p^2*d^2*k2*p*s^3*u^2 - K1p^2*K2p^2*d^2*p^2*s^3*u^2 - 2*K1p*K2p^2*d^2*p^3*s^3*u^2 - K2p^2*d^2*p^4*s^3*u^2 - K1p^2*K2p*d^2*k2*p*s^4*u^2 - K1p^2*K2p*d^2*p^2*s^4*u^2 + 2*K1p*K2p*d^2*k2*p^2*s^4*u^2 - K1p*K2p*d^2*p^3*s^4*u^2 + K1p^2*K2p^2*d*k2*p*s^3*u^3 + K1p^2*K2p*d*k2*p*s^4*u^3 - 2*K2p*d*k2*p^3*s^4*u^3 + K1p^2*d*k2*p*s^5*u^3 - K1p*d*k2*p^2*s^5*u^3 - K1p^2*k2^2*s^5*u^4 + 2*K1p*k2^2*p*s^5*u^4 - k2^2*p^2*s^5*u^4 - K1p^2*K2p^3*d^4*p^2*s + K1p^2*K2p^2*d^4*p^2*s^2 - 2*K1p^2*K2p^2*d^4*p*s^3 + 2*K1p^2*K2p^3*d^3*p^2*s*u + 2*K1p*K2p^3*d^3*p^3*s*u - 2*K1p*K2p^2*d^3*p^3*s^2*u + 2*K1p^2*K2p^2*d^3*p*s^3*u + K1p^2*K2p^2*d^2*k2*p*s^3*u - 2*K1p^2*K2p*d^3*p^2*s^3*u + 4*K1p*K2p^2*d^3*p^2*s^3*u + K1p^2*K2p*d^3*p*s^4*u + K1p^2*K2p*d^2*k2*p*s^4*u - 2*K1p*K2p*d^2*k2*p^2*s^4*u - K1p^2*K2p^3*d^2*p^2*s*u^2 - 2*K1p*K2p^3*d^2*p^3*s*u^2 - K2p^3*d^2*p^4*s*u^2 - K1p^2*K2p^2*d^2*k2*p*s^2*u^2 - K1p^2*K2p^2*d^2*p^2*s^2*u^2 + 2*K1p*K2p^2*d^2*k2*p^2*s^2*u^2 + K2p^2*d^2*p^4*s^2*u^2 - K1p^2*K2p^2*d*k2*p*s^3*u^2 + K1p^2*K2p*d^2*k2*p*s^3*u^2 + 2*K1p^2*K2p*d^2*p^2*s^3*u^2 - 2*K1p*K2p^2*d^2*p^2*s^3*u^2 - 2*K1p*K2p*d^2*k2*p^2*s^3*u^2 + 2*K1p*K2p*d^2*p^3*s^3*u^2 - 2*K2p^2*d^2*p^3*s^3*u^2 - 2*K1p^2*K2p*d^2*k2*s^4*u^2 - K1p^2*K2p*d*k2*p*s^4*u^2 + 2*K1p*K2p*d^2*k2*p*s^4*u^2 - K1p*K2p*d^2*p^2*s^4*u^2 + 2*K2p*d*k2*p^3*s^4*u^2 - K1p^2*d*k2*p*s^5*u^2 + K1p*d*k2*p^2*s^5*u^2 + K1p^2*K2p^2*d*k2*p*s^2*u^3 - 2*K2p^2*d*k2*p^3*s^2*u^3 + K1p^2*K2p*d*k2*p*s^3*u^3 - 2*K1p*K2p*d*k2*p^2*s^3*u^3 + 2*K2p*d*k2*p^3*s^3*u^3 - 2*K1p^2*d*k2*p*s^4*u^3 + 2*K1p*K2p*d*k2*p*s^4*u^3 + 2*K1p*d*k2*p^2*s^4*u^3 - 2*K2p*d*k2*p^2*s^4*u^3 + 2*K1p^2*k2^2*s^5*u^3 - 4*K1p*k2^2*p*s^5*u^3 + 2*k2^2*p^2*s^5*u^3 - K1p^2*K2p*k2^2*s^3*u^4 + 2*K1p*K2p*k2^2*p*s^3*u^4 - K2p*k2^2*p^2*s^3*u^4 + K1p^2*k2^2*s^4*u^4 - 2*K1p*k2^2*p*s^4*u^4 + k2^2*p^2*s^4*u^4 - 2*K1p^2*K2p^3*d^4*p*s + 2*K1p^2*K2p^2*d^4*p*s^2 - K1p^2*K2p^2*d^4*s^3 + K1p^2*K2p^3*d^3*p^2*u + 2*K1p^2*K2p^3*d^3*p*s*u - 2*K1p^2*K2p^2*d^3*p^2*s*u + 4*K1p*K2p^3*d^3*p^2*s*u + K1p^2*K2p^2*d^2*k2*p*s^2*u + K1p^2*K2p*d^3*p^2*s^2*u - 4*K1p*K2p^2*d^3*p^2*s^2*u - 2*K1p*K2p^2*d^2*k2*p^2*s^2*u - 2*K1p^2*K2p*d^3*p*s^3*u + 2*K1p*K2p^2*d^3*p*s^3*u - K1p^2*K2p*d^2*k2*p*s^3*u + 2*K1p*K2p*d^2*k2*p^2*s^3*u + 2*K1p^2*K2p*d^2*k2*s^4*u - 2*K1p*K2p*d^2*k2*p*s^4*u - K1p^2*K2p^3*d^2*p^2*u^2 - K1p*K2p^3*d^2*p^3*u^2 + 2*K1p^2*K2p^2*d^2*p^2*s*u^2 - 2*K1p*K2p^3*d^2*p^2*s*u^2 + 2*K1p*K2p^2*d^2*p^3*s*u^2 - 2*K2p^3*d^2*p^3*s*u^2 - 2*K1p^2*K2p^2*d^2*k2*s^2*u^2 - K1p^2*K2p^2*d*k2*p*s^2*u^2 + 2*K1p*K2p^2*d^2*k2*p*s^2*u^2 - K1p^2*K2p*d^2*p^2*s^2*u^2 - K1p*K2p*d^2*p^3*s^2*u^2 + 2*K2p^2*d^2*p^3*s^2*u^2 + 2*K2p^2*d*k2*p^3*s^2*u^2 + 2*K1p^2*K2p*d^2*k2*s^3*u^2 - K1p^2*K2p*d*k2*p*s^3*u^2 - 2*K1p*K2p*d^2*k2*p*s^3*u^2 + 2*K1p*K2p*d^2*p^2*s^3*u^2 - K2p^2*d^2*p^2*s^3*u^2 + 2*K1p*K2p*d*k2*p^2*s^3*u^2 - 2*K2p*d*k2*p^3*s^3*u^2 + 2*K1p^2*d*k2*p*s^4*u^2 - 2*K1p*K2p*d*k2*p*s^4*u^2 - 2*K1p*d*k2*p^2*s^4*u^2 + 2*K2p*d*k2*p^2*s^4*u^2 - K1p^2*k2^2*s^5*u^2 + 2*K1p*k2^2*p*s^5*u^2 - k2^2*p^2*s^5*u^2 + K1p^2*K2p^2*d*k2*p*s*u^3 - K1p*K2p^2*d*k2*p^2*s*u^3 - 2*K1p^2*K2p*d*k2*p*s^2*u^3 + 2*K1p*K2p^2*d*k2*p*s^2*u^3 + 2*K1p*K2p*d*k2*p^2*s^2*u^3 - 2*K2p^2*d*k2*p^2*s^2*u^3 + 2*K1p^2*K2p*k2^2*s^3*u^3 + K1p^2*d*k2*p*s^3*u^3 - 2*K1p*K2p*d*k2*p*s^3*u^3 - 4*K1p*K2p*k2^2*p*s^3*u^3 - K1p*d*k2*p^2*s^3*u^3 + 2*K2p*d*k2*p^2*s^3*u^3 + 2*K2p*k2^2*p^2*s^3*u^3 - 2*K1p^2*k2^2*s^4*u^3 + 4*K1p*k2^2*p*s^4*u^3 - 2*k2^2*p^2*s^4*u^3 - K1p^2*K2p^3*d^4*s + K1p^2*K2p^2*d^4*s^2 + K1p^2*K2p^3*d^3*p*u - 2*K1p^2*K2p^2*d^3*p*s*u + 2*K1p*K2p^3*d^3*p*s*u + 2*K1p^2*K2p^2*d^2*k2*s^2*u + K1p^2*K2p*d^3*p*s^2*u - 2*K1p*K2p^2*d^3*p*s^2*u - 2*K1p*K2p^2*d^2*k2*p*s^2*u - 2*K1p^2*K2p*d^2*k2*s^3*u + 2*K1p*K2p*d^2*k2*p*s^3*u - K1p*K2p^3*d^2*p^2*u^2 - K1p^2*K2p^2*d*k2*p*s*u^2 + 2*K1p*K2p^2*d^2*p^2*s*u^2 - K2p^3*d^2*p^2*s*u^2 + K1p*K2p^2*d*k2*p^2*s*u^2 + 2*K1p^2*K2p*d*k2*p*s^2*u^2 - 2*K1p*K2p^2*d*k2*p*s^2*u^2 - K1p*K2p*d^2*p^2*s^2*u^2 + K2p^2*d^2*p^2*s^2*u^2 - 2*K1p*K2p*d*k2*p^2*s^2*u^2 + 2*K2p^2*d*k2*p^2*s^2*u^2 - K1p^2*K2p*k2^2*s^3*u^2 - K1p^2*d*k2*p*s^3*u^2 + 2*K1p*K2p*d*k2*p*s^3*u^2 + 2*K1p*K2p*k2^2*p*s^3*u^2 + K1p*d*k2*p^2*s^3*u^2 - 2*K2p*d*k2*p^2*s^3*u^2 - K2p*k2^2*p^2*s^3*u^2 + K1p^2*k2^2*s^4*u^2 - 2*K1p*k2^2*p*s^4*u^2 + k2^2*p^2*s^4*u^2;
              2    2  4  2  3        2    2  3  2  3
numpol := -K1p  K2p  d  p  s  + 2 K1p  K2p  d  p  s  u

          2    2  2       3  2      2    2  2  2  3  2      2    2         3  3
     - K1p  K2p  d  k2 p s  u  - K1p  K2p  d  p  s  u  + K1p  K2p  d k2 p s  u

          2      3  2  4        2      2       4  2      2      2  2  4  2
     + K1p  K2p d  p  s  u - K1p  K2p d  k2 p s  u  - K1p  K2p d  p  s  u

          2             4  3      2         5  3      2   2  5  4
     + K1p  K2p d k2 p s  u  + K1p  d k2 p s  u  - K1p  k2  s  u

                2  3  3  3              2  2  3  3  2
     + 2 K1p K2p  d  p  s  u - 2 K1p K2p  d  p  s  u

                  2     2  4  2            2  3  4  2             2  5  3
     + 2 K1p K2p d  k2 p  s  u  - K1p K2p d  p  s  u  - K1p d k2 p  s  u

               2    5  4      2  2  4  3  2               3  4  3
     + 2 K1p k2  p s  u  - K2p  d  p  s  u  - 2 K2p d k2 p  s  u

         2  2  5  4      2    3  4  2          2    3  3  2
     - k2  p  s  u  - K1p  K2p  d  p  s + 2 K1p  K2p  d  p  s u

          2    3  2  2    2      2    2  4  2  2        2    2  4    3
     - K1p  K2p  d  p  s u  + K1p  K2p  d  p  s  - 2 K1p  K2p  d  p s

            2    2  3    3        2    2  2       3        2    2  2       2  2
     + 2 K1p  K2p  d  p s  u + K1p  K2p  d  k2 p s  u - K1p  K2p  d  k2 p s  u

          2    2  2  2  2  2      2    2         3  2      2    2         2  3
     - K1p  K2p  d  p  s  u  - K1p  K2p  d k2 p s  u  + K1p  K2p  d k2 p s  u

            2      3  2  3        2      3    4        2      2       4
     - 2 K1p  K2p d  p  s  u + K1p  K2p d  p s  u + K1p  K2p d  k2 p s  u

          2      2       3  2        2      2     4  2        2      2  2  3  2
     + K1p  K2p d  k2 p s  u  - 2 K1p  K2p d  k2 s  u  + 2 K1p  K2p d  p  s  u

          2             4  2      2             3  3      2       2  3  4
     - K1p  K2p d k2 p s  u  + K1p  K2p d k2 p s  u  - K1p  K2p k2  s  u

          2         5  2        2         4  3        2   2  5  3
     - K1p  d k2 p s  u  - 2 K1p  d k2 p s  u  + 2 K1p  k2  s  u

          2   2  4  4            3  3  3                3  2  3    2
     + K1p  k2  s  u  + 2 K1p K2p  d  p  s u - 2 K1p K2p  d  p  s u

                2  3  3  2              2  3  2  3
     - 2 K1p K2p  d  p  s  u + 4 K1p K2p  d  p  s  u

                2  2     2  2  2            2  2  2  3  2
     + 2 K1p K2p  d  k2 p  s  u  - 2 K1p K2p  d  p  s  u

                  2     2  4                2     2  3  2
     - 2 K1p K2p d  k2 p  s  u - 2 K1p K2p d  k2 p  s  u

                  2       4  2              2  3  3  2            2  2  4  2
     + 2 K1p K2p d  k2 p s  u  + 2 K1p K2p d  p  s  u  - K1p K2p d  p  s  u

                       2  3  3                     4  3               2    3  4
     - 2 K1p K2p d k2 p  s  u  + 2 K1p K2p d k2 p s  u  + 2 K1p K2p k2  p s  u

                 2  5  2               2  4  3           2    5  3
     + K1p d k2 p  s  u  + 2 K1p d k2 p  s  u  - 4 K1p k2  p s  u

               2    4  4      3  2  4    2      2  2  4  2  2
     - 2 K1p k2  p s  u  - K2p  d  p  s u  + K2p  d  p  s  u

            2  2  3  3  2        2       3  2  3               3  4  2
     - 2 K2p  d  p  s  u  - 2 K2p  d k2 p  s  u  + 2 K2p d k2 p  s  u

                   3  3  3               2  4  3         2  2  3  4
     + 2 K2p d k2 p  s  u  - 2 K2p d k2 p  s  u  - K2p k2  p  s  u

           2  2  5  3     2  2  4  4        2    3  4          2    3  3  2
     + 2 k2  p  s  u  + k2  p  s  u  - 2 K1p  K2p  d  p s + K1p  K2p  d  p  u

            2    3  3            2    3  2  2  2        2    2  4    2
     + 2 K1p  K2p  d  p s u - K1p  K2p  d  p  u  + 2 K1p  K2p  d  p s

          2    2  4  3        2    2  3  2          2    2  2       2
     - K1p  K2p  d  s  - 2 K1p  K2p  d  p  s u + K1p  K2p  d  k2 p s  u

            2    2  2     2  2        2    2  2  2    2
     - 2 K1p  K2p  d  k2 s  u  + 2 K1p  K2p  d  p  s u

          2    2         2  2      2    2           3      2      3  2  2
     - K1p  K2p  d k2 p s  u  + K1p  K2p  d k2 p s u  + K1p  K2p d  p  s  u

            2      3    3        2      2       3          2      2     4
     - 2 K1p  K2p d  p s  u - K1p  K2p d  k2 p s  u + 2 K1p  K2p d  k2 s  u

            2      2     3  2      2      2  2  2  2      2             3  2
     + 2 K1p  K2p d  k2 s  u  - K1p  K2p d  p  s  u  - K1p  K2p d k2 p s  u

            2             2  3        2       2  3  3        2         4  2
     - 2 K1p  K2p d k2 p s  u  + 2 K1p  K2p k2  s  u  + 2 K1p  d k2 p s  u

          2         3  3      2   2  5  2        2   2  4  3
     + K1p  d k2 p s  u  - K1p  k2  s  u  - 2 K1p  k2  s  u

                3  3  2              3  2  3  2            3  2  2    2
     + 4 K1p K2p  d  p  s u - K1p K2p  d  p  u  - 2 K1p K2p  d  p  s u

                2  3  2  2              2  3    3              2  2     2  2
     - 4 K1p K2p  d  p  s  u + 2 K1p K2p  d  p s  u - 2 K1p K2p  d  k2 p  s  u

                2  2       2  2            2  2  3    2          2       2    3
     + 2 K1p K2p  d  k2 p s  u  + 2 K1p K2p  d  p  s u  - K1p K2p  d k2 p  s u

                2         2  3              2     2  3
     + 2 K1p K2p  d k2 p s  u  + 2 K1p K2p d  k2 p  s  u

                  2       4                2       3  2            2  3  2  2
     - 2 K1p K2p d  k2 p s  u - 2 K1p K2p d  k2 p s  u  - K1p K2p d  p  s  u

                  2  2  3  2                   2  3  2
     + 2 K1p K2p d  p  s  u  + 2 K1p K2p d k2 p  s  u

                       2  2  3                     4  2
     + 2 K1p K2p d k2 p  s  u  - 2 K1p K2p d k2 p s  u

                         3  3               2    3  3               2  4  2
     - 2 K1p K2p d k2 p s  u  - 4 K1p K2p k2  p s  u  - 2 K1p d k2 p  s  u

                 2  3  3           2    5  2           2    4  3
     - K1p d k2 p  s  u  + 2 K1p k2  p s  u  + 4 K1p k2  p s  u

            3  2  3    2        2  2  3  2  2      2  2  2  3  2
     - 2 K2p  d  p  s u  + 2 K2p  d  p  s  u  - K2p  d  p  s  u

            2       3  2  2        2       2  2  3               3  3  2
     + 2 K2p  d k2 p  s  u  - 2 K2p  d k2 p  s  u  - 2 K2p d k2 p  s  u

                   2  4  2               2  3  3           2  2  3  3
     + 2 K2p d k2 p  s  u  + 2 K2p d k2 p  s  u  + 2 K2p k2  p  s  u

         2  2  5  2       2  2  4  3      2    3  4        2    3  3
     - k2  p  s  u  - 2 k2  p  s  u  - K1p  K2p  d  s + K1p  K2p  d  p u

          2    2  4  2        2    2  3              2    2  2     2
     + K1p  K2p  d  s  - 2 K1p  K2p  d  p s u + 2 K1p  K2p  d  k2 s  u

          2    2           2      2      3    2          2      2     3
     - K1p  K2p  d k2 p s u  + K1p  K2p d  p s  u - 2 K1p  K2p d  k2 s  u

            2             2  2      2       2  3  2      2         3  2
     + 2 K1p  K2p d k2 p s  u  - K1p  K2p k2  s  u  - K1p  d k2 p s  u

          2   2  4  2            3  3                3  2  2  2
     + K1p  k2  s  u  + 2 K1p K2p  d  p s u - K1p K2p  d  p  u

                2  3    2              2  2       2              2  2  2    2
     - 2 K1p K2p  d  p s  u - 2 K1p K2p  d  k2 p s  u + 2 K1p K2p  d  p  s u

              2       2    2            2         2  2              2       3
     + K1p K2p  d k2 p  s u  - 2 K1p K2p  d k2 p s  u  + 2 K1p K2p d  k2 p s  u

                2  2  2  2                   2  2  2                     3  2
     - K1p K2p d  p  s  u  - 2 K1p K2p d k2 p  s  u  + 2 K1p K2p d k2 p s  u

                   2    3  2             2  3  2           2    4  2
     + 2 K1p K2p k2  p s  u  + K1p d k2 p  s  u  - 2 K1p k2  p s  u

          3  2  2    2      2  2  2  2  2        2       2  2  2
     - K2p  d  p  s u  + K2p  d  p  s  u  + 2 K2p  d k2 p  s  u

                   2  3  2         2  2  3  2     2  2  4  2
     - 2 K2p d k2 p  s  u  - K2p k2  p  s  u  + k2  p  s  u


> ineqs := [op(newineq1), op(newineq2), numpol]:

> infolevel[Groebner:-Basis]:=0:

We expect the semi-algebraic set defined by map(_p->_p > 0, ineqs) to be empty This computations looks hard. sols:=RAG[HasRealSolutions](map(_p->_p > 0, ineqs)): numpol has degree 2 in K1p ; then one eliminates K1p like in a weak CAD. We are lucky: these polynomials can be factored.

> res:=factor(resultant(numpol, newineq1[2], K1p))/(K2p*d*p^2);
             3  2  3          2  2  3                 3  2            2  3  2
res := -K2p d  p  s  + 4 K2p d  p  s  u - K2p d k2 p s  u  - 4 K2p d p  s  u

                 3  3    2  2  4             4  2        2  4  2         4  3
     + K2p k2 p s  u  + d  p  s  u + d k2 p s  u  - 2 d p  s  u  - k2 p s  u

          2  3  2          2  2  2            2    2    2        3  2  2
     - K2p  d  p  s + 4 K2p  d  p  s u - 4 K2p  d p  s u  + K2p d  p  s

              3    3          2  2  2            2    3                 3
     - 2 K2p d  p s  - 2 K2p d  p  s  u + 6 K2p d  p s  u + K2p d k2 p s  u

                   2  2              3  2             3  2             2  3
     + K2p d k2 p s  u  - 4 K2p d p s  u  - K2p k2 p s  u  - K2p k2 p s  u

          2  2  3      2    4             4             3  2        2  3  2
     - 2 d  p  s  u + d  p s  u - d k2 p s  u - d k2 p s  u  + 4 d p  s  u

            4  2         4  2         3  3        2  3          2  2  2
     - d p s  u  + k2 p s  u  + k2 p s  u  - 2 K2p  d  p s + K2p  d  p  u

            2  2              2    2  2        2        2          3    2
     + 6 K2p  d  p s u - 2 K2p  d p  u  - 4 K2p  d p s u  + 2 K2p d  p s

            3  3          2  2              2    2            2  3
     - K2p d  s  - 2 K2p d  p  s u - 4 K2p d  p s  u + 2 K2p d  s  u

                   2              2    2              2  2          3  2
     - K2p d k2 p s  u + 4 K2p d p  s u  + 2 K2p d p s  u  - K2p d s  u

                 2  2    2  2  2        2    3             3          2  2  2
     + K2p k2 p s  u  + d  p  s  u - 2 d  p s  u + d k2 p s  u - 2 d p  s  u

              3  2         3  2      2  3        2  2            2  2
     + 2 d p s  u  - k2 p s  u  - K2p  d  s + K2p  d  p u + 2 K2p  d  s u

          2      2      2      2        3  2          2                2  2
     - K2p  d p u  - K2p  d s u  + K2p d  s  - 2 K2p d  p s u - 2 K2p d  s  u

                    2          2  2    2    2          2  2
     + 2 K2p d p s u  + K2p d s  u  + d  p s  u - d p s  u

> disc:=factor(discrim(numpol, K1p))/((s^2+K2p-s)*(K2p*d*p+k2*s*u+K2p*d-k2*s)^2*d*p^3*u^3):

One can get rid of factors which we know are sign invariant because of the other inequalities.

> newineqs:=map(_p->if not(member(K1p, indets(_p))) then _p fi, [op(newineq1), op(newineq2)]):

res, disc and polynomials in newineqs have degree 1 in k2. Then, we also eliminate it. Those polynomials must be positive (we omit k2>0)

> newineqs2:=[u, 1 - u, K2p - s, d, -d + u, s, -s + 1, p];
           newineqs2 := [u, 1 - u, K2p - s, d, -d + u, s, -s + 1, p]

It remains to deal with newineqs[1] which has degree 1 in k2. Its leading coeff dominant is -s*(u-1) and is sign invariant. Its constant coeff has degree 0 and is -K2p*d*(p+1). p+1 is >0 (because p>0), s also and K2p also, hence no constraint coming from newineqs[1].

It remains to treat res and disc. We start with res. leading coeff is -p*s^2*u*(u-1)*(s-1)*(d-u)*(K2p-s). All its factors are sign invariant modulo the above inequalities. Constant term is more involved and brings new constraints.

> qres:=map(_p->_p[1], factors(factor(coeff(res, k2, 0))/(d*(s^2+K2p-s)))[2]);
qres := [d p - 2 p u + d - u,

                                 2
    K2p d p s - 2 K2p p s u - p s  u + K2p d s - K2p p u - K2p s u + p s u]

One can check that d*p-2*p*u+d-u is negative modulo newineqs2.

> qres:={op(qres)} minus {d*p-2*p*u+d-u};
qres :=

                                  2
    {K2p d p s - 2 K2p p s u - p s  u + K2p d s - K2p p u - K2p s u + p s u}

> qres:=[op(qres)];
qres :=

                                  2
    [K2p d p s - 2 K2p p s u - p s  u + K2p d s - K2p p u - K2p s u + p s u]

We study now disc as above. Leading coeff -4*K2p*s^3*(u-1)*(s-1)*(d-u)*(K2p-s) is sign invariant modulo newineqs2. Same for constant term which is d*p*u*(s^2+K2p-s)^3.

Finally it remains to study the resultant de res and disc (eliminating k2). factor(resultant(res, disc, k2)); is -d*s^2*(u-1)*(s-1)*(d-u)*(s^2+K2p-s)*(K2p-s)*(2*K2p*d*p*s-4*K2p*p*s*u-p*s^2*u+2*K2p*d*s-K2p*p*u-2*K2p*s*u+p*s*u)^2 Sign invariant modulo newineqs2. Then it remains to solve:

> sols:=RAG[PointsPerComponents]([op(map(_p->_p>0, newineqs2)), op(map(_p->_p<>0, qres))]):
RAGlib release 3.54 (2021)
Mohab Safey El Din (Mohab.Safey@lip6.fr)
Sorbonne University, Paris, France
<-Path set to /home/safey/Recherche/Software/raglib/fgb/linux/current/FGblib/../libfgbunid.so
    FGb/Maple interface package Version 1.70
    JC Faugere (Jean-Charles.Faugere@inria.fr)
    Type ?FGb for documentation
remove directly 0 pols (elim=0)
[...]
<-remove directly 0 pols (elim=0)

>> lprint(sols);
[[K2p = 1, d = 149521988969/274877906944, p = 126991278303/274877906944, s =
170417339309/274877906944, u = 149934875065/274877906944]]

## One finds for instance
> sols:=[[K2p = 1, d = 149521988969/274877906944, p = 126991278303/274877906944, s = 170417339309/274877906944, u = 149934875065/274877906944]];
                       149521988969      126991278303      170417339309
sols := [[K2p = 1, d = ------------, p = ------------, s = ------------,
                       274877906944      274877906944      274877906944

        149934875065
    u = ------------]]
        274877906944


> newsols:=[];
                                 newsols := []

> newsols2:=[];
                                 newsols2 := []

> for i from 1 to nops(sols) do
> pol:=primpart(mul(_p, _p in subs(sols[i], [op(newineqs), k2, res, disc])));
> rroots:=map(rhs, RootFinding[Isolate](pol, output=midpoint));
> rats:=find_rationals(rroots);
> newsols:=[op(newsols), op(map(r->[k2=r, op(sols[i])], rats))];
>   for j from 1 to nops(newsols) do
>     pol:=primpart(mul(_p, _p in subs(newsols[j], ineqs)));
>     rroots:=map(rhs, RootFinding[Isolate](pol, output=midpoint));
>     rats:=find_rationals(rroots);
>     newsols2:=[op(newsols2), op(map(r->[K1p=r, op(newsols[j])], rats))];
>   od:
> od;
                                                                2
pol := (21292459058018746831611 k2 - 60088279883482951540343) k2  (285213291\
    006574533704412395051472525699287293777309797839895367981988329552238557\

    66399031967744000 k2 - 1105414254127672263994281646220111978336705860396\

    33747299112285425789279213032630168085055465869614227) (5613539793088465\
    868905040971147827178373438318436231856832532228576624635928760908590522\

    5708994560 k2 + 36566302589113111167332848769911798446046750710304325338\

    862355195441751332124457614231119867694065417)

rroots :=

     -6008055972219972459393     213227699562065608954835  2130713862655685
    [-----------------------, 0, ------------------------, ----------------]
       9223372036854775808       75557863725914323419136     549755813888

                      rats := [-653, -326, 1, 1939, 3877]

                                     149521988969      126991278303
newsols := [[k2 = -653, K2p = 1, d = ------------, p = ------------,
                                     274877906944      274877906944

        170417339309      149934875065                            149521988969
    s = ------------, u = ------------], [k2 = -326, K2p = 1, d = ------------,
        274877906944      274877906944                            274877906944

        126991278303      170417339309      149934875065
    p = ------------, s = ------------, u = ------------], [k2 = 1, K2p = 1,
        274877906944      274877906944      274877906944

        149521988969      126991278303      170417339309      149934875065
    d = ------------, p = ------------, s = ------------, u = ------------], [
        274877906944      274877906944      274877906944      274877906944

                            149521988969      126991278303      170417339309
    k2 = 1939, K2p = 1, d = ------------, p = ------------, s = ------------,
                            274877906944      274877906944      274877906944

        149934875065                            149521988969      126991278303
    u = ------------], [k2 = 3877, K2p = 1, d = ------------, p = ------------,
        274877906944                            274877906944      274877906944

        170417339309      149934875065
    s = ------------, u = ------------]]
        274877906944      274877906944

> finalsols:=[];
                                finalsols := []

> for i from  1 to nops(newsols2) do
> if not(member(-1, map(sign, subs(newsols2[i], ineqs)))) then
>   finalsols:=[op(finalsols), newsols2[i]];
> fi;
> od;
> quit
memory used=2246.9MB, alloc=145.3MB, time=55.74

Conclusion

finalsols is empty: the semi-algebraic set is empty, and stability is checked.