External Python Libraries for Engineering
(for internal use only)
Habiburrahman (Zulfikri), [Link]., [Link]., Ph.D.
Chemical Engineering Department,
Faculty of Engineering, Universitas Indonesia
About Me
Name : Habiburrahman
Alias : Habiburrahman Zulfikri
Birth details : Bukittinggi, November 26th 1989
Email : habib@[Link]
Education
2007–2011 Bachelor: Institut Teknologi Bandung, Indonesia
2011–2013 Master : Rijksuniversiteit Groningen, The Netherlands
Università degli Studi di Perugia, Italy
Université Paul Sabatier, Toulouse III, France
2013–2018 Doctor : Universiteit Twente, The Netherlands
Employment
2019–2021 Postdoctoral researcher at Háskóli Íslands, Iceland
2021–now Lecturer at Universitas Indonesia, Indonesia
2
Learning
Learning is the process of acquiring
new understanding, knowledge,
behaviors, skills, values, attitudes,
and preferences.
• Be brave, be confident, don't be shy
• There are no stupid questions
• Try harder, don't despair
3
Books
4
Agenda
Lecture 8: External python libraries for engineering
Lecture 9: Using Python to solve root finding problem employing the secant method
Lecture 10: Using Python for numerical integration with the Simpson's rule
Lecture 11: Using Python to solve ordinary differential equation
Lecture 12: Introduction to MATLAB/GNU Octave
Lecture 13: Application of MATLAB/GNU Octave in transport phenomena
Lecture 14: Application of MATLAB/GNU Octave in chemical reaction engineering
Note:
• The above agenda may slightly change.
• No class on May 3rd 2022.
5
Technicalities
Component of grades
Exercises = 25%
Final exam = 25%
Total = 50%
All interactions will be in MS Teams unless otherwise specified.
6
External libraries
There are many external python libraries out there, e.g., numpy,
matplotlib, scipy, sympy, pandas, seaborn, etc.
numpy: the foundation of numerical array storage in Python. The
numpy library adds powerful linear algebra data structures to Python.
It allows us to construct and manipulate vectors and tensors very
efficiently.
matplotlib: basic plotting library in python.
7
Numpy: Array and Vector Creation ([Link])
import numpy import numpy as np
myvector = [Link]([5,3,7]) myvector = [Link]([5,3,7])
myarray = [Link]([[2,3],[6.7,1.0]]) myarray = [Link]([[2,3],[6.7,1.0]])
print(myvector) print(myvector)
print(myarray) print(myarray)
print([Link]) print([Link])
print([Link]) print([Link])
• Every numpy array has a “data type” or “dtype”
• the first array, “myvector”, was constructed using
[5 3 7] a list of integers
[[2. 3. ] • The second array, “myarray”, was constructed
[6.7 1. ]]
int32
using a mixture of integers and floating-point
float64 numbers so the dtype was set as “float64”
8
Numpy: Array and Vector Creation ([Link])
• it is desirable to specify the data type when it is first created.
import numpy
myvector = [Link]([5,3,7], dtype=numpy.float64)
print(myvector)
print([Link])
[5. 3. 7.]
float64
9
Numpy: Array and Vector Creation ([Link])
[Link]() returns a new array of given shape and type, filled with zeros.
>>> [Link](5)
array([ 0., 0., 0., 0., 0.]) >>> s = (2,2)
>>> [Link](s)
>>> [Link]((2, 1)) array([[ 0., 0.],
array([[ 0.], [ 0., 0.]])
[ 0.]])
An example of allocating space for an array and then overwriting the initial values.
size = 3
myarray = [Link]([size,size])
[[1. 1. 1. ]
for i in range(size):
[1. 0.5 0.33333333]
for j in range(size):
[1. 0.33333333 0.2 ]]
myarray[i,j]=1.0/(i*j+1.0)
print(myarray) 10
Numpy: Array and Vector Creation ([Link])
size = 3
myarray = [Link]([size,size])
[[1. 1. 1. ]
for i in range(size):
[1. 0.5 0.33333333]
for j in range(size):
[1. 0.33333333 0.2 ]]
myarray[i,j]=1.0/(i*j+1.0)
print(myarray)
numpy arrays are indexed starting with zero! This is important, and forgetting how
numpy arrays are indexed leads to many troublesome bugs in the code!
11
Numpy: Array and Vector Creation ([Link])
>>> myarray = [Link](1,11)
>>> print(myarray)
[ 1 2 3 4 5 6 7 8 9 10]
• An array of length 10 that contains the integers from 1 to 10.
• The first entry in the vector (in this case “1”) is at index zero.
>>> myarray[3:6] = [Link]([300,400,500])
>>> print(myarray)
[ 1 2 3 300 400 500 7 8 9 10]
• We replace the values stored at indices 3, 4, and 5 because the slicing command – 3:6 –
does not include the last index (6) in the slice.
• Index 3 initially contains the value 4, which is replaced with 300.
12
Numpy: Array and Vector Creation (General)
>>> lin = [Link](1.0,3.0,6)
>>> print(lin)
[1. 1.4 1.8 2.2 2.6 3. ]
• The linspace() function creates a numpy vector with a starting value as the first
number and a final stopping value as the second number.
• The entries in the vector between the starting and stopping values are evenly spaced,
and the total number of entries in the vector can be specified using a third number.
>>> logger = [Link](1.0,3.0,num=5)
>>> print(logger)
[ 10. 31.6227766 100. 316.22776602 1000. ]
• We create a vector starting at 101 and ending at 103 with length 5.
• The intermediate entries in the vector are based on linearly spaced exponents,
e.g., 101.5 = 31.62.
13
Numpy: Array and Vector Creation (example)
import numpy
seqLength = 10
seq = [Link](seqLength,dtype=numpy.int32)
seq[0]=0
seq[1]=1
for i in range(2,seqLength):
seq[i]=seq[i-1]+seq[i-2]
print("Final sequence: ",seq)
Final sequence: [ 0 1 1 2 3 5 8 13 21 34]
14
Numpy: Array operations
>>> myarray = [Link](5) >>> yourarray = [Link](5)
>>> print(myarray) >>> print(yourarray)
[0 1 2 3 4] [1. 1. 1. 1. 1.]
>>> print([Link]) >>> theirarray = myarray - 3*yourarray
(5,) >>> print(theirarray)
>>> myarray = myarray*4 [-3. 1. 5. 9. 13.]
>>> print(myarray) >>> print([Link](myarray,theirarray))
[ 0 4 8 12 16] 360.0
>>> itsarray = [Link](myarray,theirarray)
>>> print(itsarray)
[[ -0. 0. 0. 0. 0.]
myarray has 5 columns an no rows. [-12. 4. 20. 36. 52.]
[-24. 8. 40. 72. 104.]
There are more possible operations [-36. 12. 60. 108. 156.]
[-48. 16. 80. 144. 208.]]
including matrix inversion, eigenvalue >>> print([Link])
calculations. (5, 5)
15
Numpy: Getting help
>>> help([Link])
logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None, axis=0)
Return numbers spaced evenly on a log scale.
In linear space, the sequence starts at ``base ** start``
(`base` to the power of `start`) and ends with ``base ** stop``
(see `endpoint` below).
.. versionchanged:: 1.16.0
Non-scalar `start` and `stop` are now supported.
Parameters
----------
start : array_like
......... (continued)
16
Numpy: Mathematical functions
>>> g = [Link](2,4,0.5)
>>> print(g)
[2. 2.5 3. 3.5]
>>> h = [Link](g)
>>> print(h)
[ 0.90929743 0.59847214 0.14112001 -0.35078323]
>>> j = [Link](g,2.5)
>>> print(j)
[ 5.65685425 9.88211769 15.58845727 22.91765149]
Comparison with the use of math library of python
import math
h = [Link](g)
print(h)
will generate an error because the sin() function in the math library is expecting a
single number and not an array of numbers.
17
Numpy: Random vectors
import numpy as np [6 5 5 6 8]
[[0.6188863 0.67311716]
a = [Link](0,10, size = 5) [0.07445555 0.70210533]]
print(a) [[6.18886305 6.7311716 ]
b = [Link](size=(2,2)) [0.74455553 7.02105329]]
print(b) [2 3 0 1 4]
print(b*10.0)
c=[Link](5)
[Link](c)
[3 1 8 5 4]
print(c)
[[0.97847571 0.97470508]
[0.72859013 0.30731679]]
Note: [[9.78475706 9.74705081]
[7.2859013 3.0731679 ]]
• Random integers can and do repeat [0 4 3 1 2]
• In [Link](), a new array is not
returned, but, instead, the original array passed into
the function is forever shuffled.
18
Numpy: Sorting and Searching
>>> import numpy as np
>>> a = [Link](4)
>>> print(a)
[0.48896493 0.80166636 0.20071366 0.0535259 ]
>>> b = [Link](a)
>>> print(b)
[0.0535259 0.20071366 0.48896493 0.80166636]
>>> print([Link](a))
0.8016663616788768
>>> print([Link](a))
1 Why does this return 1?
>>> print([Link](a))
0.0535259
>>> print([Link](a))
(array([0, 1, 2, 3], dtype=int64),)
19
Numpy: Polynomials
Polynomials can be represented using the [Link] package.
Polynomials must be written with the zero-order term first and then progressing
sequentially to higher-order terms.
>>> import [Link] as np
>>> import [Link] as plt
>>> f = [Link]([-3., -2., 1.])
>>> print([Link]())
[-1. 3.]
>>> (x,y) = [Link](8,domain=[-2,5])
>>> [Link](x,y)
>>> [Link]('[Link]')
20
Numpy: Loading and saving arrays
Numpy includes extensive support for writing vectors and arrays into files and then
loading those files at another time.
>>> import numpy as np
>>> x = [Link](5)
>>> print(x) The content of [Link]
[0 1 2 3 4]
>>> [Link]('binFile',x) # .npy extension auto added 0.000000000000000000e+00
>>> y=[Link]('[Link]') # binary file load 1.000000000000000000e+00
>>> print(2.0*y) 2.000000000000000000e+00
[0. 2. 4. 6. 8.] 3.000000000000000000e+00
>>> [Link]('[Link]',x) # extension required 4.000000000000000000e+00
>>> z = [Link]('[Link]')
>>> print()
[Link]() function saves the array to a binary file.
[Link]() function is used to save the original array to an ASCII text file.
21
Matplotlib library
it is recommended that the user import [Link]
>>> import numpy >>> import numpy
>>> import [Link] as plt >>> import [Link] as plt
>>> x = [Link](0,10, num=100) >>> x = [Link](0,10, num=100)
>>> y = [Link](x) >>> y = [Link](x)
>>> [Link](x,y) >>> [Link](x,y,'bo’)
>>> [Link]() >>> [Link]()
22
Matplotlib library: multiple plot
import [Link] as plt
import numpy as np
x = [Link](0, 2, 20)
[Link](x, x, label='linear')
[Link](x, x**2, '.’, label='quadratic')
[Link](x, x**3, '--', label='cubic')
[Link]('x-axis label')
[Link]('y-axis label')
[Link]("Polynomials")
[Link]()
[Link]('[Link]',dpi=150)
23
Matplotlib library: contour plot
import pylab
import numpy
def f(x,y):
return (1-x/2+x**2+y**3)*[Link](-x**2-y**2)
n = 256
x = [Link](-2,4,n)
y = [Link](-2,4,n)
X,Y = [Link](x,y)
C = [Link](X, Y, f(X,Y), 8)
[Link](C,inline=1)
[Link](C,orientation='vertical')
[Link]('[Link]',dpi=150)
The numpy meshgrid function extends the one-
dimensional vectors over a two-dimensional array.
24
Application: Gillespie algorithm
The Gillespie algorithm is based on counting the
exact number of molecules in the system for a
preset number of reactions.
The algorithm is based on the generation of
random numbers (i.e., analogous to rolling a dice)
for two calculations:
(1) using the reaction rate and a random number,
the algorithm determines whether or not a
reaction occurred for a random molecular collision
(2) using a random number to discretely
approximate the time until the next collision.
25
Application: Gillespie algorithm
import numpy
import [Link] as plt
k1 = 0.1 # forward rate constant
k2 = 0.01 # reverse rate constant
maxReact = 1000 # maximum number of reaction to simulate
numMol = [Link]((2,maxReact),dtype = [Link])
timePt = [Link](maxReact, dtype = [Link])
numMol[0,0] = 175 # initial number of A's
numMol[1,0] = 25 # initial number of B's
timePt[0] = 0.0 # initial time
rands = [Link](2,maxReact)
26
Application: Gillespie algorithm
for i in range(maxReact-1):
proB = k1*numMol[0,i] # probability of forming B
proA = k2*numMol[1,i] # probability of forming A
dt = -[Link](rands[0,i])/(proB+proA) # time till next reaction
timePt[i+1] = timePt[i] + dt
if rands[1,i] < (proB/(proA+proB)): # check to see if we form B
numMol[0,i+1] = numMol[0,i] - 1.0
numMol[1,i+1] = numMol[1,i] + 1.0
else: # else we form A
numMol[0,i+1] = numMol[0,i] + 1.0
numMol[1,i+1] = numMol[1,i] - 1.0
[Link](timePt,numMol[0,:], label="A")
[Link](timePt,numMol[1,:], label="B")
[Link]('time')
[Link]('number of molecules')
[Link]()
[Link]('[Link]',dpi=150)
27
Assignment
1. Create smooth curves of A and B obtained with the deterministic method and
superimposed them with the Gillespie’s curves.
2. Create the profile of the number of molecules of A, I and P for the following
consecutive reaction! Assume the following: ka = 0.1, kb = 0.01, maximum
number of reaction is 1000, and there are initially 200 A molecules.
Superimpose also the smooth curves of each species resulted from the
application of the deterministic method
3. Repeat assignment 2 for using ka = 0.01 and kb = 0.2.
28
THANK YOU
29