In [1]:

    

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
%matplotlib inline

def dX_dt(t, X):
    R = X[0]
    F = X[1]
    alpha, beta, gamma, delta = parameters
    dRdt = alpha * R - beta * F * R
    dFdt = gamma * R * F - delta * F
    return np.array([dRdt,dFdt])

# initial condition
X0 = np.array([20.,10.])
parameters = (.2, .01, .001, .1)
extra = 2
t = np.linspace(0,300,1000)
dX_dt(0., X0)


    

    
Out[1]:





array([ 2. , -0.8])

In [2]:

    

r = integrate.ode(dX_dt).set_integrator('zvode', method='bdf')
r.set_initial_value(X0, 0.)
dt = 1
x = np.empty((300, 2))
while r.successful() and r.t < 300:
    x[int(r.t) - 1, :] =r.integrate(r.t+dt)

plt.plot(x)


    

    




/home/lorenzo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: ComplexWarning: Casting complex values to real discards the imaginary part
  






    
Out[2]:





[<matplotlib.lines.Line2D at 0x7f8acd184588>,
 <matplotlib.lines.Line2D at 0x7f8acd184748>]






    

In [3]:

    

X0 = np.vstack((X0, np.eye(X0.shape[0]))).flatten()
X0


    

    
Out[3]:





array([ 20.,  10.,   1.,   0.,   0.,   1.])

In [43]:

    

def dfx(t, X):
    R = X[0]
    F = X[1]
    alpha, beta, gamma, delta = parameters
    return np.dot(np.array([alpha - beta * F, - beta * R, gamma * F, gamma * R -  delta]).reshape((2, 2)), X[2:].reshape((2, 2))).flatten()

def df(t, X):
    var = X[:2]
    return np.hstack((dX_dt(t, var), dfx(t, X)))
dfx(0, X0)


    

    
Out[43]:





array([ 0.1 , -0.2 ,  0.01, -0.08])

In [44]:

    

df(0., X0)


    

    
Out[44]:





array([ 2.  , -0.8 ,  0.1 , -0.2 ,  0.01, -0.08])

In [45]:

    

r = integrate.ode(df).set_integrator('zvode', method='bdf')
r.set_initial_value(X0, 0.)
dt = 1
x = np.zeros((300, 6))
i=1
x[0, :] = X0
print(r.integrate(r.t+dt))
while r.successful() and r.t < 300:
    x[i, :] = r.integrate(r.t+dt)
    i+=1


    

    




[  2.21889697e+01+0.j   9.24097480e+00+0.j   1.10835942e+00+0.j
  -2.13247893e-01+0.j   9.72848818e-03+0.j   9.23133472e-01+0.j]






    




/home/lorenzo/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:9: ComplexWarning: Casting complex values to real discards the imaginary part
  if __name__ == '__main__':

In [46]:

    

x.shape


    

    
Out[46]:





(300, 6)

In [52]:

    

plt.plot(x[:, 2:])


    

    
Out[52]:





[<matplotlib.lines.Line2D at 0x7f8accccc748>,
 <matplotlib.lines.Line2D at 0x7f8accfa7eb8>,
 <matplotlib.lines.Line2D at 0x7f8accfa7588>,
 <matplotlib.lines.Line2D at 0x7f8accfa7208>]






    

In [50]:

    

x[-1, 2:]


    

    
Out[50]:





array([ 0.72820157,  1.69434038, -2.4751122 , -4.0916437 ])

In [ ]: