In [2]:

    

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns


    

In [3]:

    

tmax=10.0
M=100
t=np.linspace(0,tmax,M)
t


    

    
Out[3]:





array([  0.        ,   0.1010101 ,   0.2020202 ,   0.3030303 ,
         0.4040404 ,   0.50505051,   0.60606061,   0.70707071,
         0.80808081,   0.90909091,   1.01010101,   1.11111111,
         1.21212121,   1.31313131,   1.41414141,   1.51515152,
         1.61616162,   1.71717172,   1.81818182,   1.91919192,
         2.02020202,   2.12121212,   2.22222222,   2.32323232,
         2.42424242,   2.52525253,   2.62626263,   2.72727273,
         2.82828283,   2.92929293,   3.03030303,   3.13131313,
         3.23232323,   3.33333333,   3.43434343,   3.53535354,
         3.63636364,   3.73737374,   3.83838384,   3.93939394,
         4.04040404,   4.14141414,   4.24242424,   4.34343434,
         4.44444444,   4.54545455,   4.64646465,   4.74747475,
         4.84848485,   4.94949495,   5.05050505,   5.15151515,
         5.25252525,   5.35353535,   5.45454545,   5.55555556,
         5.65656566,   5.75757576,   5.85858586,   5.95959596,
         6.06060606,   6.16161616,   6.26262626,   6.36363636,
         6.46464646,   6.56565657,   6.66666667,   6.76767677,
         6.86868687,   6.96969697,   7.07070707,   7.17171717,
         7.27272727,   7.37373737,   7.47474747,   7.57575758,
         7.67676768,   7.77777778,   7.87878788,   7.97979798,
         8.08080808,   8.18181818,   8.28282828,   8.38383838,
         8.48484848,   8.58585859,   8.68686869,   8.78787879,
         8.88888889,   8.98989899,   9.09090909,   9.19191919,
         9.29292929,   9.39393939,   9.49494949,   9.5959596 ,
         9.6969697 ,   9.7979798 ,   9.8989899 ,  10.        ])

In [5]:

    

h=t[1]-t[0]
print ("h=",h)


    

    




h= 0.10101010101

In [6]:

    

N=2
y=np.zeros((M,N))
print("N=",N)
print("M=",M)
print ("y.shape=",y.shape)


    

    




N= 2
M= 100
y.shape= (100, 2)

In [7]:

    

from scipy.integrate import odeint


    

In [8]:

    

def derivs(yvec,t,alpha,beta,delta,gamma):
    x=yvec[0]
    y=yvec[1]
    dx=alpha*x-beta*x*y
    dy=delta*x*y-gamma*y
    return np.array([dx,dy])


    

In [9]:

    

nfoxes=10
nrabbits=20
ic=np.array([nrabbits,nfoxes])
maxt=20.0
alpha=1.0
beta=0.1
delta=0.1
gamma=1.0


    

In [12]:

    

t=np.linspace(0,maxt,int(100*maxt))
soln=odeint(derivs,ic,t,args=(alpha,beta,delta,gamma),atol=1e-9,rtol=1e-8)


    

In [13]:

    

plt.plot(t, soln[:,0], label='rabbits')
plt.plot(t, soln[:,1], label='foxes')
plt.xlabel('t')
plt.ylabel('count')
plt.legend();


    

In [14]:

    

plt.plot(soln[:,0], soln[:,1])
plt.xlim(0, 25)
plt.ylim(0, 25)
plt.xlabel('rabbits')
plt.ylabel('foxes');


    

In [ ]: