I have a system of 5 non linear ODEs and I am trying to find the steady states of the system by solving the following. The variables are a,b,c,d,e.
eq=a^3 + (-1 + b) b c e + a^2 (-1 + b + b c + e + d x1 + x3) +
a (b^2 c + b (c (-1 + e) + e) + e (-1 + d x1 + x2 + x3)) == 0
b^3 (-1 + c) + b^2 (-1 + c) (-1 + a + e) + a b d x1 + a d e x1 +
b e (1 + a (-1 + c) - c - x4) == 0
a d x5 + a b (-1 + c) x7 + (-1 + b) b c x7 == c x6 + (-1 + b) b x7
b (-1 + c) x7 + b^2 (x7 - c x7) + x8 == a d x5 + a b (-1 + c) x7 + d x8
(a + b) (-e x10 + x9) == e x11
with conditions on parameters and variables as,
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 <= a <= 1, 0<= b <= 1, 0 <= c <= 1, 0<= d <= 1, 0 <= e <= y
As it takes a long time to solve the full system analytically, I am solving the last three equations in eq first and then substituting these solutions to the first 2 equations to get the final solution the entire system.
This is the Matlab code that I used.
syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 y a b c d e;
eq1=a^3 + (-1 + b)*b*c*e + a^2 *(-1 + b + b*c + e + d*x1 + x3)+a *(b^2 *c + b *(c *(-1 + e) + e)...
+ e *(-1 + d *x1 + x2 + x3));
eq2=b^3*(-1 + c) + b^2 *(-1 + c)* (-1 + a + e) + a* b* d* x1 + a *d* e *x1 + b*e*(1 + a *(-1 + c) ...
- c - x4);
eq3=a*d *x5 + a* b *(-1 + c)* x7 + (-1 + b)* b* c *x7 - c* x6 - (-1 + b)* b* x7;
eq4= b* (-1 + c)* x7 + b^2 *(x7 - c *x7) + x8 -a *d* x5 + a *b* (-1 + c)* x7 + d* x8;
eq5=(a + b)* (-e *x10 + x9) -e *x11;
assume(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<=a &a<=1 & 0<=e & e<=1 & 0<=b & b<=1 & 0<=c & c<=1 & 0<=d & d<=1)
sol = solve([eq3,eq4,eq5],[c,d,e],'ReturnConditions', true);
cval=sol.c;
dval=sol.d;
eval=sol.e;
bqe=subs(eq1,[c,d,e],[cval,dval,eval]);
beqe=subs(eq2,[c,d,e],[cval,dval,eval]);
sol2= solve([bqe,beqe],[a,b],'ReturnConditions', true);
aval=sol2.a
bval=sol2.b
But what it outputs for c,d and e are
cval =
-(- 2*x5*x7*a^2*b + x7*x8*a*b + x5*x8*a + x7*x8*b^2 - x7*x8*b)/(2*x5*x7*a^2*b - x7*x8*a*b - x5*x6*a - x7*x8*b^2 + x7*x8*b + x6*x8)
dval =
(b^2*x7*x8 - b^2*x6*x7 - x6*x8 + b*x6*x7 - b*x7*x8 + a*b*x6*x7 + a*b*x7*x8)/(2*x5*x7*a^2*b - x7*x8*a*b - x5*x6*a - x7*x8*b^2 + x7*x8*b + x6*x8)
eval =
(x9*(a + b))/(x11 + a*x10 + b*x10)
and for a and b its just sym variables as z and z1.
The Mathematica output equivalent to c,d and e is
((a == 0 && c == (b (-1 + a + b) x7)/(-x6 + b (-1 + a + b) x7) &&
d == (b (-1 + a + c - a c) x7 + b^2 (x7 - c x7) + x8)/(
a x5 + x8)) || (((b == 0 && a > 0 &&
c x6 (a x5 + x8) == a x5 x8) || (0 < b <
1 && ((c == (((-1 + b) b x7 +
a (-x5 + b x7)) x8)/(-(x6 - (-1 + b) b x7) x8 +
a (-x5 x6 + b x7 x8)) && (x5 == (b (-1 + a + b) x7)/a ||
x6 < (b (-1 + a + b) x7 x8)/(a x5 + x8) ||
x6 > (b (-1 + a + b) x7 x8)/(a x5 + x8)) && (a + b >=
1 || a > 0) && (a + b < 1 ||
x6 < (b (-1 + a + b) x7 x8)/(a x5 + x8) || a + b > 1 ||
x6 > (b (-1 + a + b) x7 x8)/(a x5 + x8))) || (x5 == (
b (-1 + a + b) x7)/a &&
x6 == (b (-1 + a + b) x7 x8)/(
a x5 + x8) && ((a > 0 && a + b < 1) ||
a + b > 1)))) || (b == 1 &&
a > 0 && (c == (a (x5 - x7) x8)/(a x5 x6 + x6 x8 - a x7 x8) ||
x6 (a x5 + x8) == a x7 x8) && (x5 == x7 ||
x6 (a x5 + x8) != a x7 x8))) &&
d == (b (-1 + a + b) x7 + c (x6 - b (-1 + a + b) x7))/(a x5))) &&
e == ((a + b) x9)/(a x10 + b x10 + x11)
1) Why do Matlab and Mathematica solutions differ for c,d,e. Without imposing conditions on parameters (as I am looking for general solutions) Mathematica gives two solutions for c,d,e while only one set of solutions are given by Matlab. Also, how can I simplify the Matlab output (similar to Mathematica)
2) In regards to Matlab outputs z and z1I read in here that "If the expression within the first parameter to RootOf() has z to a power 5 or higher, you will not be able to find analytic solutions (not unless some other part of the expression can be induced to cancel out other parts so that an analytical can be found.) In such a case, you will need to do numeric solving, or (if the situation warrants) convert the RootOf() into an appropriate call to roots()."
After substituting the values to of c,d,e to the first two equations in eq it turns out to be a polynomial of 5th order. So, is it possible to find an analytical solution to this system (in Matlab or Mathematica) ? What is meant by
convert the RootOf() into an appropriate call to roots()