Ordinary Differential Equation

Last updated: July 22nd, 20202020-07-22Project preview

Before we start...

We need to install packages required to calculate. To do so, we need to use Terminal as following instruction.

It should be installed but just in case, if you attempt to execute this script at somewhere else

  1. (assuming you are in the same env as "open lab") Open Terminal from Launcher

image

  1. Click Terminal then in the prompt, type the followings

image

pip install scipy matplotlib ipywidgets ipympl

Or execute the following line:

In [1]:
!pip install scipy matplotlib
Requirement already satisfied: scipy in /usr/local/lib/python3.8/site-packages (1.5.1)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.8/site-packages (3.3.0)
Requirement already satisfied: numpy>=1.14.5 in /usr/local/lib/python3.8/site-packages (from scipy) (1.19.1)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.8/site-packages (from matplotlib) (0.10.0)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /usr/local/lib/python3.8/site-packages (from matplotlib) (2.4.7)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.8/site-packages (from matplotlib) (2.8.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.8/site-packages (from matplotlib) (1.2.0)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.8/site-packages (from matplotlib) (7.2.0)
Requirement already satisfied: six in /usr/local/lib/python3.8/site-packages (from cycler>=0.10->matplotlib) (1.14.0)
WARNING: You are using pip version 20.0.2; however, version 20.1.1 is available.
You should consider upgrading via the '/usr/local/bin/python -m pip install --upgrade pip' command.
In [2]:
## These are the python libraries you will utilize for the calculation

import matplotlib.pylab as plt 
import matplotlib.pyplot as plt; plt.rcdefaults()

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

%matplotlib inline

import numpy as np
import scipy
from scipy.optimize import curve_fit
from scipy.integrate import odeint

Ordinary Differential Equation

image

$\dfrac{d[A]}{dt}=-k_1[A]+k_2[B]$

$\dfrac{d[B]}{dt}=-(k_2+k_3)[B]+k_1[A]+k_4[C]$

$\dfrac{d[C]}{dt}=-(k_4+k_5)[C]+k_3[B]+k_6[D]$

$\dfrac{d[D]}{dt}=-k_6[D]+k_5[C]$

  1. create a function to return $\dfrac{dy}{dt}$
  2. set parameters (ex. $k_1 = 1$)
  3. set length of calculation (this is a time-dependent)
  4. execute ode solver
  5. plot them
In [3]:
def func1(y,t):
    
    k1 = 1
    k2 = 0.5
    k3 = 1.2
    k4 = 0.7
    k5 = 0.8
    k6 = 0.7
    
    dAdt = -k1*y[0] + k2*y[1]
    dBdt = -(k2+k3)*y[1] + k1*y[0] + k4*y[2]
    dCdt = -(k4+k5)*y[2] + k3*y[1] + k6*y[3]
    dDdt = -k6*y[3] + k5*y[2]
    
    dydt = [dAdt, dBdt, dCdt, dDdt]
    
    return dydt 
In [4]:
def func2(y,t,k):
    
    k1 = k[0]
    k2 = k[1]
    k3 = k[2]
    k4 = k[3]
    k5 = k[4]
    k6 = k[5]
    
    dAdt = -k1*y[0] + k2*y[1]
    dBdt = -(k2+k3)*y[1] + k1*y[0] + k4*y[2]
    dCdt = -(k4+k5)*y[2] + k3*y[1] + k6*y[3]
    dDdt = -k6*y[3] + k5*y[2]
    
    dydt = [dAdt, dBdt, dCdt, dDdt]
    
    return dydt 
In [5]:
def func3(y,t,k):
    
    A,B,C,D = y
    k1,k2,k3,k4,k5,k6 = k
    
    dAdt = -k1*A + k2*B
    dBdt = -(k2+k3)*B + k1*A + k4*C
    dCdt = -(k4+k5)*C + k3*B + k6*D
    dDdt = -k6*D + k5*C
    
    dydt = [dAdt, dBdt, dCdt, dDdt]
    
    return dydt 
In [6]:
def func4(y,t,J):
    return np.dot(J,y)
In [7]:
# Initial condition of calculation
y0 = [10,0,0,0]

# for function 1 
# the set of parameter is predetermined 

# for function 2,3
# k (reaction rate or transition rate) is one of input variables you need to specify before calculation
# function 2 and 3 handle the list of reaction rate
k = [1, 0.5, 1.2, 0.7, 0.8, 0.7]

# for function 4 
# we are using Jacobian method so we will introduce Jacobian matrix using np.array and the function will perform dot operation
k1 = k[0]
k2 = k[1]
k3 = k[2]
k4 = k[3]
k5 = k[4]
k6 = k[5]
Jacobian = np.array([[-k1, k2,       0, 0],
                     [k1, -(k2+k3), k4, 0],
                     [0, k3, -(k4+k5), k6],
                     [0, 0,       k5, -k6]])

print(Jacobian)
[[-1.   0.5  0.   0. ]
 [ 1.  -1.7  0.7  0. ]
 [ 0.   1.2 -1.5  0.7]
 [ 0.   0.   0.8 -0.7]]
In [8]:
# calculation time
time_start = 0
time_end = 10
NoOfSteps = 100
time = np.linspace(time_start,time_end,NoOfSteps)
In [9]:
y1 = odeint(func1,y0,time)
y2 = odeint(func2,y0,time,args=(k,))
y3 = odeint(func3,y0,time,args=(k,))
y4 = odeint(func4,y0,time,args=(Jacobian,))
In [10]:
print('shape of the outcome y1 =', np.shape(y1))
print('shape of the outcome y2 =', np.shape(y2))
print('shape of the outcome y3 =', np.shape(y3))
print('shape of the outcome y4 =', np.shape(y4))
shape of the outcome y1 = (100, 4)
shape of the outcome y2 = (100, 4)
shape of the outcome y3 = (100, 4)
shape of the outcome y4 = (100, 4)
In [11]:
plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y1[:,0],'r-',label='func1')
plt.plot(time,y2[:,0],'b-',label='func2')
plt.plot(time,y3[:,0],'g-',label='func3')
plt.plot(time,y4[:,0],'m-',label='func4')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[11]:
<matplotlib.legend.Legend at 0x7fe2f474e370>
In [12]:
plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y1[:,0],'r-',label='A')
plt.plot(time,y1[:,1],'b-',label='B')
plt.plot(time,y1[:,2],'g-',label='C')
plt.plot(time,y1[:,3],'m-',label='D')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[12]:
<matplotlib.legend.Legend at 0x7fe2c3e41a30>

Use of Dictionary

In [13]:
def genJacobian(k):
    k1 = k['k1']
    k2 = k['k2']
    k3 = k['k3']
    k4 = k['k4']
    k5 = k['k5']
    k6 = k['k6']
    
    Jacobian = np.array([[-k1, k2,       0, 0],
                         [k1, -(k2+k3), k4, 0],
                         [0, k3, -(k4+k5), k6],
                         [0, 0,       k5, -k6]])
    
    return Jacobian
In [14]:
kdict = {'k1':1,'k2':0.5,'k3':1.2,'k4':0.7,'k5':0.8,'k6':0.7}
J = genJacobian(kdict)

y4_1 = odeint(func4,y0,time,args=(J,))

plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y4_1[:,0],'r-',label='A')
plt.plot(time,y4_1[:,1],'b-',label='B')
plt.plot(time,y4_1[:,2],'g-',label='C')
plt.plot(time,y4_1[:,3],'m-',label='D')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[14]:
<matplotlib.legend.Legend at 0x7fe2c3ebd1c0>

Let's try something simpler

Creating your own library/package using .py file

In [15]:
import functions as fn

y1 = odeint(fn.func1,y0,time)
y2 = odeint(fn.func2,y0,time,args=(k,))
y3 = odeint(fn.func3,y0,time,args=(k,))
y4 = odeint(fn.func4,y0,time,args=(Jacobian,))
In [16]:
plt.figure(figsize=(6,4),dpi=100)
plt.tick_params(direction='in',labelsize=12)
plt.plot(time,y1[:,0],'r-',label='A')
plt.plot(time,y1[:,1],'b-',label='B')
plt.plot(time,y1[:,2],'g-',label='C')
plt.plot(time,y1[:,3],'m-',label='D')
plt.xlabel('Time [a.u.]',fontsize=12)
plt.ylabel('Conc. [a.u.]',fontsize=12)
plt.legend(loc='best',fontsize=10)
Out[16]:
<matplotlib.legend.Legend at 0x7fe2c3ea1ac0>
In [ ]:
 
Notebooks AI
Notebooks AI Profile20060