Why Is Sympy Nsolve Giving Me Different Values Than What Is Graphically Shown?

I am trying to find the values of parameters that solve a system of differential equations (find when the rates are zero, and the values no longer change). I iterated the system wi

Solution 1:

Nonlinear systems are notorious for their sensitivity to initial guesses when using a numeric solver. In this case you have done nothing wrong. If you use initial guesses closer to what you see on the screen you might get a solution where P does not go to zero:

guess N,R,H,P= (1.5, 0.7, 0.6, 0.1)
N = 1.50905265512538
R = 0.749587544253425
H = 0.571428571428571
P = 0.129589598408824

guess N,R,H,P= (1.51, 0.7, 0.6, 0.1)
N = 1.66676383218043
R = 0.571428571428571
H = 0.597439764807826
P = -1.21722045507780e-18

In cases like these it is nice to be able to make a systematic tweaking of the input variables. For that, something like tweak shown below can be used:

def tweak(v, p=.1, d=None):
    """Given a vector of initial guesses, return
    vectors with values that are modified by a factor of +/-p
    along with the unmodified values until all possibilities
    are given. d controls number of decimals shown in returned
    values; all are shown when d isNone (default).

    EXAMPLES
    ========

    >>> for i intweak((1,2),.1):
    ...     print(i)
    [1.1, 2.2]
    [1.1, 2]
    [1.1, 1.8]
    [1, 2.2]
    [1, 2]
    [1, 1.8]
    [0.9, 2.2]
    [0.9, 2]
    [0.9, 1.8]
    """"
    from sympy.utilities.iterables import cartes
    if d is None:
        d = 15
    c = []
    for i in v:
        c.append([i*(1+p),i,i*(1-p)])
    for i incartes(*c):
        yield [round(i,d) for i in i]

eqs = [
    I - E*N + (r[0]*(m[0]*R)**q[0])/(s[0]**q[0]+(m[0]*R)**q[0]) - a[0]*N*R/(b[0] + N),
    a[0]*N*R/(b[0] + N) - m[0]*R - l[0]*R - a[1]*H*R/(b[1] + R),
    a[1]*H*R/(b[1] + R) - m[1]*H - a[1]*H*P/(b[2] + H),
    a[1]*H*P/(b[2] + H) - m[2]*P]

saw = set()
vv = 1.5, 0.7, 0.6, 0.1for v0 intweak(vv, d=3):
    equilibrium_values = nsolve(eqs, (N, R, H, P), v0)
    ok = round(equilibrium_values[-1], 3)
    if ok notin saw:
        saw.add(ok)
        print('guess N,R,H,P=',v0,ok)

which gives

guess N,R,H,P= [1.65, 0.77, 0.66, 0.11] 0.130
guess N,R,H,P= [1.65, 0.77, 0.6, 0.1] 0.0