'Python Error: "divide by zero encountered in log" in function, how to adjust? (Beginner here)

For my current assignment, I am to establish the stability of intersection/equilibrium points between two nullclines, which I have defined as follows:

def fNullcline(F):
    P = (1/k)*((1/beta)*np.log(F/(1-F))-c*F+v)
    return P

def pNullcline(P):
    F = (1/delta)*(pD-alpha*P+(r*P**2)/(m**2+P**2))
    return F

I also have a method "stability" that applies the Hurwitz criteria on the underlying system's Jacobian:

def dPdt(P,F):
    return pD-delta*F-alpha*P+(r*P**2)/(m**2+P**2)

def dFdt(P,F):
    return s*(1/(1+sym.exp(-beta*(-v+c*F+k*P)))-F)

def stability(P,F):
    
    x = sym.Symbol('x')
    
    ax = sym.diff(dPdt(x, F),x) 
    ddx = sym.lambdify(x, ax)
    a = ddx(P)
    
    # shortening the code here: the same happens for b, c, d
    
    matrix = [[a, b],[c,d]]
    
    eigenvalues, eigenvectors = np.linalg.eig(matrix) 
    e1 = eigenvalues[0]
    e2 = eigenvalues[1]
    
    if(e1 >= 0 or e2 >= 0):
        return 0  
    else:
        return 1

The solution I was looking for was later provided. Basically, values became too small! So this code was added to make sure no too small values are being used for checking the stability:

set={0}

for j in range(1,210):
    for i in range(1,410):
        x=i*0.005
        y=j*0.005
        x,y=fsolve(System,[x,y])
        nexist=1
        for i in set:
            if(abs(y-i))<0.00001:
                nexist=0
        if(nexist):
            set.add(y)
    set.discard(0)

I'm still pretty new to coding so the function in and on itself is still a bit of a mystery to me, but it eventually helped in making the little program run smoothly :) I would again like to express gratitude for all the help I have received on this question. Below, there are still some helpful comments, which is why I will leave this question up in case anyone might run into this problem in the future, and can find a solution thanks to this thread.



Solution 1:[1]

After a bit of back and forth, I came to realise that to avoid the log to use unwanted values, I can instead define set as an array:

set = np.arange(0, 2, 0.001)

I get a list of values within this array as output, complete with their according stabilities. This is not a perfect solution as I still get runtime errors (in fact, I now get... three error messages), but I got what I wanted out of it, so I'm counting that as a win?

Edit: I am further elaborating on this in the original post to improve the documentation, however, I would like to point out again here that this solution does not seem to be working, after all. I was too hasty! I apologise for the confusion. It's a very rocky road for me. The correct solution has since been provided, and is documented in the original question.

Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source
Solution 1