1

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.

illylilly
  • 21
  • 4
  • you could try testing if `(1-F)` is zero and replace it with something like [machine epsilon](https://stackoverflow.com/a/19141711/7540911) to have it as small as possible. I'm not a mathematician, so you'll have to figure out if that makes sense – Nullman May 15 '22 at 08:08
  • @Nullman Hi, thanks so much for your reply! I do have to admit, though, that I am not quite sure how to test whether or not `(1-F)` is zero... I have only ever worked with arrays before, this is my first time working with variables that are just kind of... there? If you could possibly give me some pointers in that direction, I'd be ever so grateful :) – illylilly May 15 '22 at 08:15
  • I have to admit, I understand nothing of the matter, but to check whether 1-F is zero, you can just put a ```print(1-F)``` in your fNullcline function just before the calculation. (I hope I got the question right. Sorry if not) – ductTapeIsMagic May 15 '22 at 08:48
  • @ductTapeIsMagic Hi! Thank so much for your comment. I actually tried your suggestion, got back some curious value, and then had an epiphany! I will post about this solution as an answer. (I think that's how it's done here? But feel free to correct me!) Thanks again, printing (1-F) may have opened up my third eye to a possible solution :) – illylilly May 15 '22 at 09:24

1 Answers1

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.

illylilly
  • 21
  • 4
  • If you append these errors, we might be able to help you also with them :). If you decide to do that, I would suggest you put them also in your original question, so it is easy for others to keep track of your progress. – ductTapeIsMagic May 15 '22 at 09:36
  • `set` is a keyword in python. I would HIGHLY recommend not overwriting keywords – Nullman May 15 '22 at 09:43
  • @ductTapeIsMagic Hi again, and thanks for the input - I've added the current progress in the original question, so if you and/or anyone else wants to help - which, again, thanks so much for that! - they can track the problem better! I hope my phrasing so far is understandable. – illylilly May 15 '22 at 10:20
  • 1
    @Nullman Hi again and also thank you so much for pointing that out - I had no idea! I've instead named the array `array` so no keywords are overridden! – illylilly May 15 '22 at 10:21