Reputation: 21
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
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