If you look at your system you will see that the 4th and 5th equations have two solutions since solving the 4th for A and substituting into the 5th gives an expression that factors as I*f(S) -- giving, for the value of A, I = 0 and S such that f(S) = 0. Judicious selection of which equation(s) to solve next and taking time to lump constants together so you don't bog down the solver gives both solutions in about 10 seconds with relatively small operation counts (relative to the results of nonlinsolve above -- 10 and 5192 operations). The process gives the same solutions for the representative values above:
def condense(eq, *x, reps=None):
"""collapse additive/multiplicative constants into single
variables, returning condensed expression and replacement
values.
Note: use of the replacement dictionary may require topological sorting
if values depend on the keys.
"""
from sympy.core.traversal import bottom_up
from sympy.simplify.radsimp import collect
from sympy.utilities.iterables import numbered_symbols
if reps is None:
reps = {}
else:
reps = {v:k for k,v in reps.items()}
con = numbered_symbols('c')
free = eq.free_symbols
def c():
while True:
rv = next(con)
if rv not in free:
return rv
def do(e):
if not e.args:
return e
e = e.func(*[do(i) for i in e.args])
isAdd=e.is_Add
if not (isAdd or e.is_Mul):
return e
if isAdd:
ee = collect(e, x, exact=None)
if ee != e:
e = do(ee)
co, id = e.as_coeff_Add() if isAdd else e.as_coeff_Mul()
i, d = id.as_independent(*x, as_Add=isAdd)
if not i.args:
return e
return e.func(co, reps.get(i, reps.setdefault(i, c())), d)
rv = do(bottom_up(eq, do))
return rv, {v: k for k, v in reps.items()}
def repsort(*replace):
"""Return sorted replacement tuples `(o, n)` such that `(o_i, n_i)`
will appear before `(o_j, n_j)` if `o_j` appears in `n_i`. An error
will be raised if `o_j` appears in `n_i` and `o_i` appears in `n_k`
if `k >= i`.
Examples
========
>>> from sympy.abc import x, y, z, a
>>> repsort((x, y + 1), (z, x + 2))
[(z, x + 2), (x, y + 1)]
>>> repsort((x, y + 1), (z, x**2))
[(z, x**2), (x, y + 1)]
>>> repsort(*Tuple((x,y+z),(y,a),(z,1/y)))
[(x, y + z), (z, 1/y), (y, a)]
Any two of the following 3 tuples will not raise an error,
but together they contain a cycle that raises an error:
>>> repsort((x, y), (y, z), (z, x))
Traceback (most recent call last):
...
raise ValueError("cycle detected")
"""
from itertools import permutations
from sympy import default_sort_key, topological_sort
free = {i for i,_ in replace}
defs, replace = sift(replace,
lambda x: x[1].is_number or not x[1].has_free(*free),
binary=True)
edges = [(i, j) for i, j in permutations(replace, 2) if
i[1].has(j[0]) and (not j[0].is_Symbol or j[0] in i[1].free_symbols)]
rv = topological_sort([replace, edges], default_sort_key)
rv.extend(ordered(defs))
return rv
def dupdate(d, s):
"""update values in d with values from s and return the combined dictionaries"""
rv = {k: v.xreplace(s) for k,v in d.items()}
rv.update(s)
return rv
# Variables of the system
syms=S, V, E, A, I, R = symbols('S, V, E, A, I, R')
# Parameters of the system
const = var('a:j k')
system = [
-A*S*c/a - I*S*c/a + R + S*(-h - j) + a*h,
A*(V*c*d/a - V*c/a) + I*(V*c*d/a - V*c/a) + S*j - V*h,
A*(-S*c*k/a + S*c/a - V*c*d/a + V*c/a) - E*h +
I*(-S*c*k/a + S*c/a - V*c*d/a + V*c/a),
A*(S*b*c*k/a - e - f - h) + I*S*b*c*k/a,
A*(-S*b*c*k/a + S*c*k/a + f) + I*(-S*b*c*k/a + S*c*k/a - g - h),
A*e + I*g + R*(-h - 1)
]
import sympy as sym
# Variables of the system
syms = S, V, E, A, I, R = sym.symbols('S, V, E, A, I, R')
# Parameters of the system
N = sym.Symbol("N", positive = True)
mi = sym.Symbol("mi", positive = True)
v = sym.Symbol("v", positive = True)
epsilon = sym.Symbol("epsilon", positive = True)
alpha = sym.Symbol("alpha", positive = True)
gamma_as = sym.Symbol("gamma_as", positive = True)
gamma_s = sym.Symbol("gamma_s", positive = True)
gamma_a = sym.Symbol("gamma_a", positive = True)
lamb = sym.Symbol("lamb", positive = True)
tau = sym.Symbol("tau", positive = True)
beta = sym.Symbol("beta", positive = True)
# Declaration of the system equations
system = [
mi*N - v*S + R - (beta*(A+I)/N)*S - mi*S,
v*S - (1-epsilon)*(beta*(A+I)/N)*V - mi*V,
(beta*(A+I)/N)*S + (1-epsilon)*(beta*(A+I)/N)*V -
sym.exp(-mi*tau)*(beta*(A+I)/N)*S - mi*E,
alpha*sym.exp(-mi*tau)*(beta*(A+I)/N)*S - (gamma_as + gamma_a + mi)*A,
(1-alpha)*sym.exp(-mi*tau)*(beta*(A+I)/N)*S + gamma_as*A - (gamma_s + mi)*I,
gamma_a*A + gamma_s*I - (1+mi)*R]
system, srep = condense(Tuple(*system), *syms)
asol = solve(system[3], A, dict=True)[0]
aeq=Tuple(*[i.xreplace(asol) for i in system])
si = solve(aeq[4], *syms, dict=True)
sol1 = dupdate(asol, si[0])
sol1 = dupdate(sol1, solve(Tuple(*system).xreplace(sol1),syms,dict=1)[0]); sol1
aeqs4 = Tuple(*[i.xreplace(si[1]) for i in aeq])
ceq, crep = condense(Tuple(*aeqs4),*syms,reps=srep)
ir = solve([ceq[0], ceq[-1]], I, R, dict=1)[0]
ve = solve([i.simplify() for i in Tuple(*ceq).xreplace(ir)], syms, dict=True)[0] # if we don't simplify to make first element 0 we get no solution -- bug?
sol2 = dupdate(asol, si[1])
sol2 = dupdate(sol2, ir)
sol2 = dupdate(sol2, ve)
crep = repsort(*crep.items())
sol1 = Dict({k:v.subs(crep) for k,v in sol1.items()}) # 10 ops
sol2 = Dict({k:v.subs(crep) for k,v in sol2.items()}) # 5192 ops
Test for specific values (as above):
>>> rep = {alpha: 1, beta: 2, epsilon: 3, gamma_as: 4, gamma_s: 5,
... gamma_a: 6, exp(-mi*tau): 7, N: 8, v: 9, mi: 10}
...
>>> sol1.xreplace(rep)
{A: 0, E: 0, I: 0, R: 0, S: 80/19, V: 72/19}
>>> sol2.xreplace(rep)
{A: -960/23, E: 89280/851, I: -256/23,
R: -640/23, S: 1200/133, V: -124200/4921}
Of course, it took time to find this path to the solution. But if the solver could make better selections of what to solve (rather than trying to get the Groebner basis of the whole system) the time for obtaining a solution from SymPy could be greatly reduced.
