# Solve a nonlinear system of equations with condition symbolically

Consider the following system of nonlinear non-dimensional ODE : \begin{align*} \frac{d\bar{x}(\tau)}{d\tau}=\quad&\bar{x}\left[1-\bar{x}-\frac{\bar {y}_{1}}{\bar{m}_{1}\bar {y}_{1}+\bar{x}}-\frac{\bar{y} _{2}}{\bar{m}_{2}y_{2}+\bar{x}} \right] +\bar{s}_{1}(\tau),\quad \bar{x}_{1}(0)>0\\ \frac{d\bar{y}_{1}(\tau)}{d\tau}=\quad& \bar{y}_{1} \Bigg[-\bar{\epsilon}_{1}\left(1+\bar{\omega}_{11}\bar {y}_{1}+\bar{\omega}_{12}\bar{y}_{2}\right)-\frac{\bar{b}_{12}}{\bar{p}_{12}\bar{y}_{2}+\bar{y}_{1}}\bar{y}_{2}+\frac{\bar{\beta}_{1}}{\bar{m}_{1}\bar {y}_{1}+\bar{x}}\bar{x}\Bigg],\quad \bar{y}_{1}(0)>0\\ \frac{d\bar{y}_{2}(\tau)}{d\tau}=\quad& \bar{y}_{2} \Bigg[-\bar{\epsilon}_{2}\left(1+\bar{\omega}_{21}\bar {y}_{1}+\bar{\omega}_{22}\bar{y}_{2}\right)+\frac{\bar{\gamma}_{12}}{\bar{p}_{12}\bar{y}_{2}+\bar{y}_{1}}\bar{y}_{1}+\frac{\bar{\beta}_{2}}{\bar{m}_{2}\bar {y}_{2}+\bar{x}}\bar{x}\Bigg],\quad \bar{y}_{2}(0)>0 \end{align*}

I tried to solve the steady sate system by matlab to find equilibrium points with s=0 as a special case:

the Matlab code:

syms x y1 y2 real
assume(x >= 0 & y1 >= 0 & y2 >= 0)
syms N R c m1 m2 s  positive;
syms epsilon1 epsilon2 beta1 beta2   positive ;
syms omega11 omega12 omega21 omega22  positive ;
syms g12 b12 p12 gamma12  positive;

F2=x*(1-x-y1/(m1*y1+x)-y2/(m2*y2+x))==0;
F6=y1*(-epsilon1*(1+omega11*y1+omega12*y2)-(b12*y2)/(p12*y2+y1)+(beta1*x)/(m1*y1+x))==0;
F7=y2*(-epsilon2*(1+omega21*y1+omega22*y2)+(gamma12*y1)/(p12*y2+y1)+(beta2*x)/(m2*y2+x))==0;
S = solve(eqns,[x,y1,y2], 'returnconditions', true);
diary('SolutionOfEquation_Subsystem_NonDim_RDFR_Condition.txt')
S.x
S.y1
S.y2

And got the follwing solution:

ans =
0
1
z
-(beta1 - epsilon1 - beta1*m1 + (epsilon1*omega11*(2*epsilon1^(3/2)*omega11 + beta1*epsilon1^(1/2)*m1^2 + beta1*m1*(4*epsilon1*omega11 - 4*beta1*omega11 + epsilon1*m1^2 + epsilon1*omega11^2 + 4*beta1*m1*omega11 - 2*epsilon1*m1*omega11)^(1/2) - 2*beta1*epsilon1^(1/2)*omega11 + beta1*epsilon1^(1/2)*m1*omega11))/(2*(epsilon1^(3/2)*omega11^2 + beta1*epsilon1^(1/2)*m1^2*omega11)))/(beta1*m1)
-(beta1 - epsilon1 - beta1*m1 + (epsilon1*omega11*(2*epsilon1^(3/2)*omega11 + beta1*epsilon1^(1/2)*m1^2 - beta1*m1*(4*epsilon1*omega11 - 4*beta1*omega11 + epsilon1*m1^2 + epsilon1*omega11^2 + 4*beta1*m1*omega11 - 2*epsilon1*m1*omega11)^(1/2) - 2*beta1*epsilon1^(1/2)*omega11 + beta1*epsilon1^(1/2)*m1*omega11))/(2*(epsilon1^(3/2)*omega11^2 + beta1*epsilon1^(1/2)*m1^2*omega11)))/(beta1*m1)
-(beta2 - epsilon2 - beta2*m2 + (epsilon2*omega22*(2*epsilon2^(3/2)*omega22 + beta2*epsilon2^(1/2)*m2^2 + beta2*m2*(4*epsilon2*omega22 - 4*beta2*omega22 + epsilon2*m2^2 + epsilon2*omega22^2 + 4*beta2*m2*omega22 - 2*epsilon2*m2*omega22)^(1/2) - 2*beta2*epsilon2^(1/2)*omega22 + beta2*epsilon2^(1/2)*m2*omega22))/(2*(epsilon2^(3/2)*omega22^2 + beta2*epsilon2^(1/2)*m2^2*omega22)))/(beta2*m2)
-(beta2 - epsilon2 - beta2*m2 + (epsilon2*omega22*(2*epsilon2^(3/2)*omega22 + beta2*epsilon2^(1/2)*m2^2 - beta2*m2*(4*epsilon2*omega22 - 4*beta2*omega22 + epsilon2*m2^2 + epsilon2*omega22^2 + 4*beta2*m2*omega22 - 2*epsilon2*m2*omega22)^(1/2) - 2*beta2*epsilon2^(1/2)*omega22 + beta2*epsilon2^(1/2)*m2*omega22))/(2*(epsilon2^(3/2)*omega22^2 + beta2*epsilon2^(1/2)*m2^2*omega22)))/(beta2*m2)

ans =
0
0
z1
-(2*epsilon1^(3/2)*omega11 + beta1*epsilon1^(1/2)*m1^2 + beta1*m1*(4*epsilon1*omega11 - 4*beta1*omega11 + epsilon1*m1^2 + epsilon1*omega11^2 + 4*beta1*m1*omega11 - 2*epsilon1*m1*omega11)^(1/2) - 2*beta1*epsilon1^(1/2)*omega11 + beta1*epsilon1^(1/2)*m1*omega11)/(2*(epsilon1^(3/2)*omega11^2 + beta1*epsilon1^(1/2)*m1^2*omega11))
-(2*epsilon1^(3/2)*omega11 + beta1*epsilon1^(1/2)*m1^2 - beta1*m1*(4*epsilon1*omega11 - 4*beta1*omega11 + epsilon1*m1^2 + epsilon1*omega11^2 + 4*beta1*m1*omega11 - 2*epsilon1*m1*omega11)^(1/2) - 2*beta1*epsilon1^(1/2)*omega11 + beta1*epsilon1^(1/2)*m1*omega11)/(2*(epsilon1^(3/2)*omega11^2 + beta1*epsilon1^(1/2)*m1^2*omega11))
0
0

ans =
0
0
z2
0
0
-(2*epsilon2^(3/2)*omega22 + beta2*epsilon2^(1/2)*m2^2 + beta2*m2*(4*epsilon2*omega22 - 4*beta2*omega22 + epsilon2*m2^2 + epsilon2*omega22^2 + 4*beta2*m2*omega22 - 2*epsilon2*m2*omega22)^(1/2) - 2*beta2*epsilon2^(1/2)*omega22 + beta2*epsilon2^(1/2)*m2*omega22)/(2*(epsilon2^(3/2)*omega22^2 + beta2*epsilon2^(1/2)*m2^2*omega22))
-(2*epsilon2^(3/2)*omega22 + beta2*epsilon2^(1/2)*m2^2 - beta2*m2*(4*epsilon2*omega22 - 4*beta2*omega22 + epsilon2*m2^2 + epsilon2*omega22^2 + 4*beta2*m2*omega22 - 2*epsilon2*m2*omega22)^(1/2) - 2*beta2*epsilon2^(1/2)*omega22 + beta2*epsilon2^(1/2)*m2*omega22)/(2*(epsilon2^(3/2)*omega22^2 + beta2*epsilon2^(1/2)*m2^2*omega22))

I checked if the solution is correct via sympy simplify(F.subs()) but found out that just the first and second one satisfies the equations and not the others. I tried again to solve the system via sympy nonlinsove or transferring the system to polynomial system via as_numer_denom() then solve via solve_poly_system or groebner()but no output after several hours.

Any suggestions?

Sympy code:

python
from sympy.interactive import printing
printing.init_printing(use_latex='mathjax')
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
x,y1,y2 = sp.symbols('x,y1,y2',real=True,nonnegative=True)
N,R,c,m1,m2,s = sp.symbols('N,R,c,m1,m2,s',real=True,positive=True)
epsilon1,epsilon2,beta1,beta2= sp.symbols('epsilon1,epsilon2,beta1,beta2',real=True,positive=True)
omega11,omega12,omega21,omega22 = sp.symbols('omega11,omega12,omega21,omega22',real=True,positive=True)
g12,b12,p12,gamma12=sp.symbols('g12,b12,p12,gamma12',real=True,positive=True)

F2=x*(1-x-(y1)/(m1*y1+x)-(y2)/(m2*y2+x))
F6=y1*(-epsilon1*(1+omega11*y1+omega12*y2)-(b12*y2)/(p12*y2+y1)+(beta1*x)/(m1*y1+x))
F7=y2*(-epsilon2*(1+omega21*y1+omega22*y2)+(gamma12*y1)/(p12*y2+y1)+(beta2*x)/(m2*y2+x))

equ=[F2,F6,F7]
sol = nonlinsolve(equ, [x, y1, y2])

other way(converting the system into multivariate polynomial system):

F_2=-m1*m2*x**2*y1*y2 + m1*m2*x*y1*y2 - m1*x**3*y1 + m1*x**2*y1 - m1*x*y1*y2 - m2*x**3*y2 + m2*x**2*y2 - m2*x*y1*y2 - x**4 + x**3 - x**2*y1 - x**2*y2
F_6=-b12*m1*y1**2*y2 - b12*x*y1*y2 + beta1*p12*x*y1*y2 + beta1*x*y1**2 - epsilon1*m1*omega11*p12*y1**3*y2 - epsilon1*m1*omega11*y1**4 - epsilon1*m1*omega12*p12*y1**2*y2**2 - epsilon1*m1*omega12*y1**3*y2 - epsilon1*m1*p12*y1**2*y2 - epsilon1*m1*y1**3 - epsilon1*omega11*p12*x*y1**2*y2 - epsilon1*omega11*x*y1**3 - epsilon1*omega12*p12*x*y1*y2**2 - epsilon1*omega12*x*y1**2*y2 - epsilon1*p12*x*y1*y2 - epsilon1*x*y1**2

F_7=beta2*p12*x*y2**2 + beta2*x*y1*y2 - epsilon2*m2*omega21*p12*y1*y2**3 - epsilon2*m2*omega21*y1**2*y2**2 - epsilon2*m2*omega22*p12*y2**4 - epsilon2*m2*omega22*y1*y2**3 - epsilon2*m2*p12*y2**3 - epsilon2*m2*y1*y2**2 - epsilon2*omega21*p12*x*y1*y2**2 - epsilon2*omega21*x*y1**2*y2 - epsilon2*omega22*p12*x*y2**3 - epsilon2*omega22*x*y1*y2**2 - epsilon2*p12*x*y2**2 - epsilon2*x*y1*y2 + gamma12*m2*y1*y2**2 + gamma12*x*y1*y2
equs=[F_2,F_6,F_7]
sol=groebner(equs, [x,y1,y2])
print(sol)
sol1=solve_poly_system(equs,x,y1,y2)
print(sol1)
• Your system of equations seems to be an ODE system, is that right? – nicoguaro Jan 4 at 18:43
• @ nicoguaro yes it is – F.O Jan 4 at 20:30