illylilly
illylilly

Reputation: 21

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.

Upvotes: 1

Views: 524

Answers (1)

illylilly
illylilly

Reputation: 21

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.

Upvotes: 1

Related Questions