# Codes for practice No. 10

###### tags: `mathbio`

## Code for plotting static phase portrait

%matplotlib inline

A = 1.0; B = 3.0; k1 = 2.0; km1 = 2.0; k2 = 1.8; k3 = 3.0

P = lambda x,y: k1*A - (km1 + k2*B)*x + k3*x**2*y
Q = lambda x,y: k2*B*x - k3*x**2*y

fun = lambda p,t: [ P(p[0],p[1]), Q(p[0],p[1]) ]

plt.figure(figsize=(10,10))

x0_ = 1.0
tt = np.linspace(0,100, 1000)
for y0_ in [0,1,2,3,4]:
    xy0 = [x0_, y0_]
    zz = odeint(fun, xy0, tt)
    plt.plot(zz[:,0], zz[:,1])
    plt.plot([xy0[0]], [xy0[1]], 'o', mec='k', ms=4)

pt = zz[-1,:]
plt.plot([pt[0]], [pt[1]], 'o', mec='k', ms=4)

# visualize the vector field
xx = np.linspace(0,3,20)
yy = np.linspace(0,3,20)
Xm,Ym = np.meshgrid(xx,yy)
plt.quiver(Xm,Ym, P(Xm,Ym), Q(Xm,Ym), scale_units='xy', angles='xy')
    
plt.legend()
plt.xlim(0,3)
plt.ylim(0,3)
plt.xlabel('X')
plt.ylabel('Y')

## Code for animating the bifuraction

%matplotlib notebook
from matplotlib.animation import FuncAnimation

A = 1.0; B = 3.0; k1 = 2.0; km1 = 2.0; k3 = 3.0
# k2 = 1.8

P = lambda x,y,k2: k1*A - (km1 + k2*B)*x + k3*x**2*y
Q = lambda x,y,k2: k2*B*x - k3*x**2*y

fun = lambda p,t,k2: [ P(*p, k2), Q(*p, k2) ]

fig, ax = plt.subplots()

line, = ax.plot([], '-')


ax.set_xlim(0,3)
ax.set_ylim(0,3)
ax.set_xlabel('X')
ax.set_ylabel('Y')

def animate(fn):
    tt = np.linspace(0,100, 1000)
    zz = odeint(fun, [1.0, 4.0], tt, args=(.8+1./100*fn,))
    line.set_data(zz[:,0], zz[:,1])
    return line

    
anim = FuncAnimation(fig, animate, frames=100, interval=50)

plt.show()

## Jacobian code

Jxx = lambda x,y: -(km1+k2*B)+2*k3*x*y
Jxy = lambda x,y: k3*x**2
Jyx = lambda x,y: k2*B - 2*k3*x*y
Jyy = lambda x,y: -k3*x**2

JJ = lambda x,y: np.array([
    [Jxx(x,y), Jxy(x,y)],
    [Jyx(x,y), Jyy(x,y)]
])

curJ = JJ(pt[0], pt[1])

display(curJ)

eivals, eivecs = np.linalg.eig(curJ)
display(eivals)

# Imports for 3D

%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
ax = fig.add_subplot(111, projection='3d')