Solving ODEs
We already know how to do basic operations such as variable operations, loops, conditions, numpy operations, plotting.
Let's take that further and learn to solve Ordinary differential equations (ODEs).
If you don't know how to solve these problems it's ok, start small and make your way up.
As the name implies, scipy is designed to bring science and python together.
We can import the whole library
import scipy
Or just the part we need
from scipy.integrate import odeint
A part that is of particular concern to use is the integration of ordinary differential equations. Something we will be doing quite a bit of in this and the upcoming tutorials.
Take for example $\dfrac{dX}{dt} = -k \cdot X(t)$
Without scipy, we would need to implement code to handle the finite difference ourselves.
def f(X, t):
dXdt = - k*X
return dXdt
Here dXdt is being treated as a new variable, and we will need to set the value of $k$.
Let's say $t$ represents time, we need an array of different instances of time from start to end. Our old friend linspace comes in handy
time_start = 0
time_finish = 1000
N_points = 50
times = np.linspace(time_start,time_finish,N_points)
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')
from scipy.integrate import odeint
def f(X, t):
dXdt = - k*X
return dXdt
k = -0.01
time_start = 0
time_finish = 1000
N_points = 50
times = np.linspace(time_start,time_finish,N_points)
initial_dependent_var = 3. #Some initial condition
calculated_dependent_var = odeint(f, initial_dependent_var, times)
plt.plot(times, calculated_dependent_var, 'ro-')
print(calculated_dependent_var[49])
plt.xlabel('Independent Variable')
plt.ylabel('Dependent variable')
plt.show()
This method of solving ODEs automatically is what we call the integrated
method, next we will learn how to solve ODEs manually with the Euler
or analytical
method.
We will be needing numpy here now more than ever:
import numpy as np
We will be taking the same example as before $\dfrac{dX}{dt} = -k \cdot X(t)$
With the same initial conditions
k = -0.01
time_start = 0
time_finish = 1000
N_points = 50
And the same linspace
times = np.linspace(time_start,time_finish,N_points)
Remember how earlier we said "without scipy, we would need to implement code to handle the finite difference ourselves". We'll let's handle the finite difference.
def f(X, t):
dXdt = - k*X
return dXdt
Here dXdt is being treated as a new variable, and we will need to set the value of $k$.
Let's say $t$ represents time, we need an array of different instances of time from start to end. Our old friend linspace comes in handy
time_start = 0
time_finish = 1000
N_points = 50
times = np.linspace(time_start,time_finish,N_points)
%reset
%matplotlib inline
import matplotlib.pyplot as plt # plotting modules
plt.style.use('ggplot') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
k = -0.01
time_start = 0
time_finish = 1000
N_points = 50
initial_dependent_var = 3. #Some initial condition
times = np.linspace(time_start,time_finish,N_points)
calculated_dependent_var = initial_dependent_var*np.exp(-k*times)
plt.plot(times, calculated_dependent_var, 'ro-')
print(calculated_dependent_var[49])
plt.xlabel('Independent Variable')
plt.ylabel('Dependent variable')
plt.show()
Let's say we have a nonlinear equation
$x^2 - log_{10}(x) = 1$
The equation can be easily solved using the method of "thorough observation"
x = 1 is the root
How do we find it numerically?
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
def func(x):
return x**2 - np.log10(x) - 1
fsolve(func, 0.5)
And that is it! Nothing more! Imagine how useful that would be for the Thermo 346 course =)
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
def func(x):
return x**2 - np.log10(x) - 1
answer = fsolve(func, 0.5)
print(answer[0])
With basic Python and ODEs, we now have the tools we need to begin to solve systems of ODEs.
Again, if you don't know how to solve these problems it's ok, start small and make your way up.
There are many ways of solving systems of equations on paper, and additionally there are many ways of solving them in python. Here we will go over two very similar hypothetical cases.
Let's say we have two ODEs with $x1$ and $x2$
$\frac{dx1}{dt} = f(x1,t) = -k \cdot x1$
and
$\frac{dx2}{dt} = g(x2,t) = \frac{-q}{k} \cdot x2$
We can treat this as a matrix with $X$ being $\begin{pmatrix} x1 \\ x2 \end{pmatrix}$
We also need our initial conditions vector $X0 = \begin{pmatrix} 10 \\ 15 \end{pmatrix}$
We first need to define our function
def func(X, t, k, q):
dXdt = np.zeros(2)
dXdt[0] = - k*X[0]
dXdt[1] = - q/k*X[1]
return dxdt
We had to initialize dXdt before using it, here the numpy.zeros
function was very helpful
And we need to define $k$ and $q$
k = 1
q = 2
We define our time intervals just as before:
time_start = 0.
time_finish = 10000.
N_points = 100
time_array = np.linspace(0, 10., N_points)
And just like that we are done:
C_num = odeint(func,x0, time_array, args=(k,q))
$k$ and $q$ were put in as arguments, but there are ways of implementing the code to avoid using the "args" parameter
$x1$ is found in the first column, and $x2$ in the second
x1 = C_num[:,0]
x2 = C_num[:,0]
import matplotlib.pyplot as plt # plotting modules
%matplotlib inline
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
#step 1 definition of the function:
def func(X, t, k, q):
dXdt = np.zeros(2)
dXdt[0] = - k*X[0]
dXdt[1] = - (q/k)*X[1]
return dXdt
k = 1
q = 2
X0 = [10.,15.]
time_start = 0.
time_finish = 10000.
N_points = 100
time_array = np.linspace(0, 10., N_points)
C_num = odeint(func, X0, time_array, args=(k,q))
plt.plot(time_array, C_num[:,0], 'r-')
plt.plot(time_array, C_num[:,0], 'go', alpha=0.8, fillstyle='none') # alpha - 80% transparency
# fillstyle - none - means no filling of the circles
plt.plot(time_array, C_num[:,1], 'b-')
plt.plot(time_array, C_num[:,1], 'k*')
plt.xlabel('t, s')
plt.ylabel('$C_X(t)$, mole/L')
plt.show()
Let's say we have two ODEs with Concentration, $C_A$, and volume, $V$
$\frac{dC_A}{dt} = f(C_A,t) = -k \cdot C_A$
and
$\frac{dV}{dt} = g(V,t) = \frac{-q}{k} \cdot V$
we can treat this as a matrix $\begin{pmatrix} C_A \\ V \end{pmatrix}$
We also need our initial conditions vector $\begin{pmatrix} C_{A0} \\ V_0 \end{pmatrix} = \begin{pmatrix} 10 \\ 15 \end{pmatrix}$
We first need to define our function
def func(X, t, k, q):
C_A = X[0]
V = X[1]
dC_Adt = -k*C_A
dVdt = -q/k*V
return [dC_Adt,dVdt]
and then we put in our values as a vector into the solver:
values = odeint(func,[C_A0, V0], time_array, args=(k,q))
And seperating them is fairly trivial
C_A = values[:,0]
V = values[:,1]
%reset
import matplotlib.pyplot as plt # plotting modules
%matplotlib inline
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
#step 1 definition of the function:
def func(X, t, k, q):
C_A = X[0]
V = X[1]
dC_Adt = -k*C_A
dVdt = -q/k*V
return [dC_Adt,dVdt]
k = 1.
q = 2.
[C_A0, V0] = [10.,15.]
#V0 = 15.
time_start = 0.
time_finish = 10000.
N_points = 100
time_array = np.linspace(0, 10., N_points)
values = odeint(func, [C_A0, V0], time_array, args=(k,q))
C_A = values[:,0]
V = values[:,1]
plt.plot(time_array, C_A, 'r-')
plt.plot(time_array, C_A, 'go', alpha=0.8, fillstyle='none') # alpha - 80% transparency
# fillstyle - none - means no filling of the circles
plt.xlabel('t, s')
plt.ylabel('$C_A(t)$, mole/L')
plt.show()
plt.plot(time_array, V, 'b-')
plt.plot(time_array, V, 'k*')
plt.xlabel('t, s')
plt.ylabel('$V(t)$, L')
plt.show()
We can look at our last two values with the following:
print(C_A[-1])
print(V[-1])
Great so now we should be ready to do some problems.
However, these problems at this stage are a bit difficult and you may find it easier to go over part 2 before coming to them.
Originally Problem #6
Imagine we live in the world of where a part of human population is infected and could take over the world.
S - Living people are Susceptible (S) victims to infection, and they can become a zombie
Z - Number of zombies, they could come from either infected people or 'resurrected' dead people
R - rate by which people die
with the following notations:
S: the number of susceptible victims
Z: the number of zombies
R: the number of people "killed"
P: the population birth rate
d: the chance of a natural death
B: the chance the "zombie disease" is transmitted (an alive person becomes a zombie)
G: the chance a dead person is resurrected into a zombie
A: the chance a zombie is totally destroyed
There are different reactions happening at the same time:
Rate of accumulation of the living people (Susceptible victims)
birth rate -> S (something is born -> S zeroth order reaction)
S -> Infected -> Z ( S + Z -> Z second order reaction with a constant B => $-r_S = B \cdot S \cdot Z$ )
S -> Dead ( S -> Dead, first with the rate $-r_S = d \cdot S$)
So the corresponding design equation will look like:
$\dfrac{\delta S}{dt} = P - B \cdot S \cdot Z - d \cdot S$
Rate of accumulation of zombies
S -> Infected -> Z ( S + Z -> Z second order reaction with a constant B => $-r_S = B \cdot S \cdot Z$ )
Dead person -> Z (resurrection of a dead person into a zombie, first order reaction with a constant R)
Z -> 0 (Killing of zombies by living people second order reaction S + Z -> 0 $-r_Z = A \cdot S \cdot Z$
$\dfrac{\delta Z}{dt} = B \cdot S \cdot Z + G \cdot R - A \cdot S \cdot Z$
Rate of dying
S -> Dead
S + Z -> 0
Dead -> Resurrected
$\dfrac{\delta R}{dt} = d \cdot S + A \cdot S \cdot Z - G \cdot R$
Initial parameters
P = 0 # birth rate
d = 0.0001 # natural death percent (per day)
B = 0.0095 # transmission percent (per day)
G = 0.0001 # resurect percent (per day)
A = 0.0001 # destroy percent (per day)
#zombie apolcalypse
import matplotlib.pyplot as plt # plotting modules
%matplotlib inline
#plt.style.use('presentation')
plt.style.use('ggplot')
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
P = 0. # birth rate
d = 0.0001 # natural death percent (per day)
B = 0.0095 # transmission percent (per day)
# # B = 0.001 # transmission percent (per day)
G = 0.0001 # resurect percent (per day)
# # G = 0.1 # resurect percent (per day)
A = 0.001 # destroy percent (per day)
# solve the system dy/dt = f(y, t)
def f(y, t):
Si = y[0]
Zi = y[1]
Ri = y[2]
# the model equations (see Munz et al. 2009)
f0 = P - B*Si*Zi - d*Si
f1 = B*Si*Zi + G*Ri - A*Si*Zi
f2 = d*Si + A*Si*Zi - G*Ri
return [f0, f1, f2]
# initial conditions
S0 = 500. # initial population
Z0 = 0 # initial zombie population
R0 = 0 # initial death population
y0 = [S0, Z0, R0] # initial condition vector
time_start = 0
time_finish = 10.
N_points = 1000
t = np.linspace(time_start, time_finish, N_points) # time grid
# solve the DEs
soln = odeint(f, y0, t)
S = soln[:, 0]
Z = soln[:, 1]
R = soln[:, 2]
# plot results
plt.figure()
plt.plot(t, S,'g--', label='Living')
plt.plot(t, Z,'r--', label='Zombies')
plt.plot(t, R,'k-', label='Dead')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
# plt.title('Zombie Apocalypse - No Init. Dead Pop.; No New Births.')
plt.legend(loc=0)
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import odeint
plt.style.use('ggplot')
from ipywidgets import interact, interactive, fixed, interact_manual
P = 0 # 10
# d = 0.0001
B = 0.0095
G = 0.001
A = 0.001
def russians_are_not_scared_of_zombies(d):
def f(C_list, time_array):
Si = C_list[0]
Zi = C_list[1]
Ri = C_list[2]
dsdt = P - B*Si*Zi - d*Si
dzdt = B*Si*Zi + G*Ri - A*Si*Zi
drdt = d*Si + A*Si*Zi - G*Ri
return [dsdt, dzdt, drdt]
S0 = 500
Z0 = 0
R0 = 0
C_list_0 = [S0, Z0, R0]
time_start = 0
time_finish = 10.
N_points = 100
time_array = np.linspace(time_start, time_finish, N_points)
C_all_num = odeint(f, C_list_0, time_array)
C_S_num = C_all_num[:,0]
C_Z_num = C_all_num[:,1]
C_R_num = C_all_num[:,2]
plt.plot(time_array, C_S_num, 'r--',label='living')
plt.plot(time_array, C_Z_num, 'b--',label='zombie')
plt.plot(time_array, C_R_num, 'g--',label='dead')
plt.xlabel('days after the outbreak')
plt.ylabel('population')
plt.show()
mymax = 0.1
mymin = 0.00000
interact(russians_are_not_scared_of_zombies, d=(mymin, mymax, (mymax-mymin)/10))
Originally problem #7
Calculate the functions $T(t), C_{A}(t)$ in a mixer:
Stirred Tank reactor:
$T_{in}, C_{A0} \to T_f, C_A$ mixer $q = \dot{V} = 100 m^3/hr$ $V = 100 m^3$
energy balance:
$V \dfrac{dT}{dt} = q (T_{f} - T(t))$
mass balance:
$V \dfrac{dC_A}{dt} = q (C_{f} - C_{A}(t))$
You need to give the program parameters $C_{A0}=0, C_{Af} = 1$ $T_0 = 350$, $T_f = 300$
you can call the python function:
def mixer(X_list, t, Tf, Cf):
...
return [dTdt, dCa/dt]
The video solution could be found here: