Notebook source code:
examples/stability/stability.ipynb
Run the notebook yourself on binder
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.
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.