I want to solve the following system of nonlinear equations. Is it possible to put condition that all variables are greater or equal zero and all parameters are positive? The variables are (x1,x2,x3,x4,y1,y2 ) and others are just parameters.
And is Maple better to solve this system than sympy?
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
x1, x2, x3, x4, y1, y2 = sp.symbols('x1, x2, x3, x4, y1, y2')
N, c1, c2, c3, c4 = sp.symbols('N, c1, c2, c3, c4')
r1, r2, r3, r4 = sp.symbols('r1, r2, r3, r4')
f11, f21, f31, f41 = sp.symbols('f11, f21, f31, f41')
f12, f22, f32, f42 = sp.symbols('f12, f22, f32, f42')
eta11, eta12, eta13, eta14 = sp.symbols('eta11, eta12, eta13, eta14')
eta21, eta22, eta23, eta24 = sp.symbols('eta21, eta22, eta23, eta24')
eta31, eta32, eta33, eta34 = sp.symbols('eta31, eta32, eta33, eta34')
eta41, eta42, eta43, eta44 = sp.symbols('eta41, eta42, eta43, eta44')
epsilon1, epsilon2, K11, K22 = sp.symbols('epsilon1, epsilon2, K11, K22')
omega1, omega2, gamma12, g12 = sp.symbols('omega1, omega2, gamma12, g12')
beta11, beta21, beta31, beta41 = sp.symbols('beta11, beta21, beta31, beta41')
beta12, beta22, beta32, beta42 = sp.symbols('beta12, beta22, beta32, beta42')
F2 = x1 * (r1 * (1 - (eta11 * x1 + eta12 * x2 + eta13 * x3 + eta14 * x4) / N) - \
f11 * y1 - f12 * y2)
F3 = x2 * (r2 * (1 - (eta21 * x1 + eta22 * x2 + eta23 * x3 + eta24 * x4) / N) - \
f21 * y1 - f22 * y2)
F4 = x3 * (r3 * (1 - (eta31 * x1 + eta32 * x2 + eta33 * x3 + eta34 * x4) / N) - \
f31 * y1 - f32 * y2)
F5 = x4 * (r4 * (1 - (eta41 * x1 + eta42 * x2 + eta43 * x3 + eta44 * x4) / N) - \
f41 * y1 - f42 * y2)
F6 = y1 * (-epsilon1 * (1 + (y1 + omega2 * y2) / K22) - g12 * y2 + beta11 * f11 * x1 + \
beta21 * f21 * x2 + beta31 * f31 * x3 + beta41 * f41 * x4)
F7 = y2 * (-epsilon2 * (1 + (omega1 * y1 + y2) / K11) +gamma12 * g12 * y1 + \
beta12 * f12 * x1 + beta22 * f22 * x2 + beta32 * f32 * x3 + beta42 * f42 * x4)
equ = (F2, F3, F4, F5, F6, F7)
sol = nonlinsolve(equ, x1, x2, x3, x4, y1, y2)
print(sol)
This is a polynomial system and we can put into standard form with
In [2]: equ = [eq.as_numer_denom()[0].expand() for eq in equ]
In [3]: for eq in equ: pprint(eq)
2
-N⋅f₁₁⋅x₁⋅y₁ - N⋅f₁₂⋅x₁⋅y₂ + N⋅r₁⋅x₁ - η₁₁⋅r₁⋅x₁ - η₁₂⋅r₁⋅x₁⋅x₂ - η₁₃⋅r₁⋅x₁⋅x₃ - η₁₄⋅r₁⋅x₁⋅x₄
2
-N⋅f₂₁⋅x₂⋅y₁ - N⋅f₂₂⋅x₂⋅y₂ + N⋅r₂⋅x₂ - η₂₁⋅r₂⋅x₁⋅x₂ - η₂₂⋅r₂⋅x₂ - η₂₃⋅r₂⋅x₂⋅x₃ - η₂₄⋅r₂⋅x₂⋅x₄
2
-N⋅f₃₁⋅x₃⋅y₁ - N⋅f₃₂⋅x₃⋅y₂ + N⋅r₃⋅x₃ - η₃₁⋅r₃⋅x₁⋅x₃ - η₃₂⋅r₃⋅x₂⋅x₃ - η₃₃⋅r₃⋅x₃ - η₃₄⋅r₃⋅x₃⋅x₄
2
-N⋅f₄₁⋅x₄⋅y₁ - N⋅f₄₂⋅x₄⋅y₂ + N⋅r₄⋅x₄ - η₄₁⋅r₄⋅x₁⋅x₄ - η₄₂⋅r₄⋅x₂⋅x₄ - η₄₃⋅r₄⋅x₃⋅x₄ - η₄₄⋅r₄⋅x₄
2
K₂₂⋅β₁₁⋅f₁₁⋅x₁⋅y₁ + K₂₂⋅β₂₁⋅f₂₁⋅x₂⋅y₁ + K₂₂⋅β₃₁⋅f₃₁⋅x₃⋅y₁ + K₂₂⋅β₄₁⋅f₄₁⋅x₄⋅y₁ - K₂₂⋅ε₁⋅y₁ - K₂₂⋅g₁₂⋅y₁⋅y₂ - ε₁⋅ω₂⋅y₁⋅y₂ - ε₁⋅y₁
2
K₁₁⋅β₁₂⋅f₁₂⋅x₁⋅y₂ + K₁₁⋅β₂₂⋅f₂₂⋅x₂⋅y₂ + K₁₁⋅β₃₂⋅f₃₂⋅x₃⋅y₂ + K₁₁⋅β₄₂⋅f₄₂⋅x₄⋅y₂ - K₁₁⋅ε₂⋅y₂ - K₁₁⋅g₁₂⋅γ₁₂⋅y₁⋅y₂ - ε₂⋅ω₁⋅y₁⋅y₂ - ε₂⋅y₂
SymPy will try to solve this using the Groebner basis but it takes a long time to compute that:
In [4]: groebner(equ, [x1,x2,x3,x4,y1,y2]) # Not sure how long this takes
I expect though that even if it does complete the result will not admit analytic solutions because solving will likely lead to polynomials of greater than order 4.
If you replace all of the parameters with concrete rational numbers then it is possible that a solution could be found but otherwise in terms of arbitrary symbols (r3 etc) I don't expect that a closed form solution will exist - if that's true then it doesn't matter whether you use Maple or SymPy or anything else.
EDIT: I now see what your system is. Each equation is of the form x1 * (a*x1 + b*x2 + ...) so it's a linear equation multiplied by one of the unknowns. That means that there are two possibilities: x1 = 0 or the linear equation is satisfied. So one solution is that x1 = x2 = ... = 0 and then there is another where none are zero. With 6 unknowns there are 64 possible solutions except some might not satisfy the non-negativity assumptions. You can find them all with
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
x1, x2, x3, x4, y1, y2 = sp.symbols('x1, x2, x3, x4, y1, y2', nonnegative=True)
N, c1, c2, c3, c4 = sp.symbols('N, c1, c2, c3, c4', positive=True)
r1, r2, r3, r4 = sp.symbols('r1, r2, r3, r4', positive=True)
f11, f21, f31, f41 = sp.symbols('f11, f21, f31, f41', positive=True)
f12, f22, f32, f42 = sp.symbols('f12, f22, f32, f42', positive=True)
eta11, eta12, eta13, eta14 = sp.symbols('eta11, eta12, eta13, eta14', positive=True)
eta21, eta22, eta23, eta24 = sp.symbols('eta21, eta22, eta23, eta24', positive=True)
eta31, eta32, eta33, eta34 = sp.symbols('eta31, eta32, eta33, eta34', positive=True)
eta41, eta42, eta43, eta44 = sp.symbols('eta41, eta42, eta43, eta44', positive=True)
epsilon1, epsilon2, K11, K22 = sp.symbols('epsilon1, epsilon2, K11, K22', positive=True)
omega1, omega2, gamma12, g12 = sp.symbols('omega1, omega2, gamma12, g12', positive=True)
beta11, beta21, beta31, beta41 = sp.symbols('beta11, beta21, beta31, beta41', positive=True)
beta12, beta22, beta32, beta42 = sp.symbols('beta12, beta22, beta32, beta42', positive=True)
F2 = (r1 * (1 - (eta11 * x1 + eta12 * x2 + eta13 * x3 + eta14 * x4) / N) - \
f11 * y1 - f12 * y2)
F3 = (r2 * (1 - (eta21 * x1 + eta22 * x2 + eta23 * x3 + eta24 * x4) / N) - \
f21 * y1 - f22 * y2)
F4 = (r3 * (1 - (eta31 * x1 + eta32 * x2 + eta33 * x3 + eta34 * x4) / N) - \
f31 * y1 - f32 * y2)
F5 = (r4 * (1 - (eta41 * x1 + eta42 * x2 + eta43 * x3 + eta44 * x4) / N) - \
f41 * y1 - f42 * y2)
F6 = (-epsilon1 * (1 + (y1 + omega2 * y2) / K22) - g12 * y2 + beta11 * f11 * x1 + \
beta21 * f21 * x2 + beta31 * f31 * x3 + beta41 * f41 * x4)
F7 = (-epsilon2 * (1 + (omega1 * y1 + y2) / K11) - gamma12 * g12 * y1 + \
beta12 * f12 * x1 + beta22 * f22 * x2 + beta32 * f32 * x3 + beta42 * f42 * x4)
equ = ((x1, F2), (x2, F3), (x3, F4), (x4, F5), (y1, F6), (y2, F7))
from itertools import product
for eqs in product(*equ):
sol = solve(eqs, [x1, x2, x3, x4, y1, y2])
pprint(sol)
That gives:
$ python t.py
{x₁: 0, x₂: 0, x₃: 0, x₄: 0, y₁: 0, y₂: 0}
[]
[]
⎧ ε₂⋅(K₁₁⋅(K₂₂⋅g₁₂ + ε₁⋅ω₂) - K₂₂⋅ε₁) ε₁⋅(-K₁₁⋅ε₂ + K₂₂⋅(K₁₁⋅g₁₂⋅γ₁₂ + ε₂⋅ω₁
⎨x₁: 0, x₂: 0, x₃: 0, x₄: 0, y₁: ───────────────────────────────────────────────, y₂: ──────────────────────────────────────────
⎩ ε₁⋅ε₂ - (K₂₂⋅g₁₂ + ε₁⋅ω₂)⋅(K₁₁⋅g₁₂⋅γ₁₂ + ε₂⋅ω₁) ε₁⋅ε₂ - (K₂₂⋅g₁₂ + ε₁⋅ω₂)⋅(K₁₁⋅g₁₂⋅γ₁₂ + ε
)) ⎫
─────⎬
₂⋅ω₁)⎭
⎧ N ⎫
⎨x₁: 0, x₂: 0, x₃: 0, x₄: ───, y₁: 0, y₂: 0⎬
⎩ η₄₄ ⎭
⎧ N⋅ε₂⋅(K₁₁⋅f₄₂ + r₄) K₁₁⋅r₄⋅(N⋅β₄₂⋅f₄₂ - ε₂⋅η₄₄)⎫
⎪x₁: 0, x₂: 0, x₃: 0, x₄: ──────────────────────────, y₁: 0, y₂: ───────────────────────────⎪
⎨ 2 2 ⎬
⎪ K₁₁⋅N⋅β₄₂⋅f₄₂ + ε₂⋅η₄₄⋅r₄ K₁₁⋅N⋅β₄₂⋅f₄₂ + ε₂⋅η₄₄⋅r₄⎪
⎩ ⎭
⎧ N⋅ε₁⋅(K₂₂⋅f₄₁ + r₄) K₂₂⋅r₄⋅(N⋅β₄₁⋅f₄₁ - ε₁⋅η₄₄) ⎫
⎪x₁: 0, x₂: 0, x₃: 0, x₄: ──────────────────────────, y₁: ───────────────────────────, y₂: 0⎪
⎨ 2 2 ⎬
⎪ K₂₂⋅N⋅β₄₁⋅f₄₁ + ε₁⋅η₄₄⋅r₄ K₂₂⋅N⋅β₄₁⋅f₄₁ + ε₁⋅η₄₄⋅r₄ ⎪
⎩ ⎭
... (continues)
The empty solutions [] correspond to the cases that are known not to satisfy the non-negativity requirement.