SciPy#
SciPy is the most used Scientific library in Python. It is well documented and easy to use, and it is fast.
It is useful for different operations, for example:
Integration
Optimization
Special functions
Interpolation
Linear Algebra
Signal Processing
Statistics
Spatial data structures and algorithms.
See all the possibilities here.
In this lesson, we will look at some examples where you use Python to solve different mathematical operations, such as integrations, equations and linear algebra.
We will mostly use the SciPy library, so please take some time to get familiar with it.
Statistics#
The scipy.stats module contains a large number of probability distributions, frequency statistics, correlation functions and statistical tests, and more.
We will not go through examples of this module at this time, but it’s a good idea to get familiar with the list of functions available.
Constants#
Instead of looking for a constant and manually assigning it to a value, you can use SciPy to find most of them!
The module you will use is scipy.constants, check it out here.
# uncomment the following line
# to print the list of available constants
import scipy.constants
scipy.constants.find()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 3
1 # uncomment the following line
2 # to print the list of available constants
----> 3 import scipy.constants
4 scipy.constants.find()
ModuleNotFoundError: No module named 'scipy'
Save the constant that you need for later use
import scipy.constants
res = scipy.constants.physical_constants["alpha particle mass"]
res
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 import scipy.constants
2 res = scipy.constants.physical_constants["alpha particle mass"]
3 res
ModuleNotFoundError: No module named 'scipy'
Integration#
SciPy can be used to solve integrals, which will be very useful to you in the future when solving chemical reactions engineering exercises, so remember these functions!
Here, we will use the quad function, but you can check out the full guide on integration here.
import scipy.integrate
# you can use the ? to get info
scipy.integrate?
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 import scipy.integrate
3 # you can use the ? to get info
4 get_ipython().run_line_magic('pinfo', 'scipy.integrate')
ModuleNotFoundError: No module named 'scipy'
Let’s define an example function that we will integrate
from math import cos, exp, pi, sin
from scipy.integrate import quad
# define a function to integrate
def f(x):
return (2 * x * pi) + 3.2
# call quad to integrate f from -2 to 2
res, err = quad(f, -2, 2)
print(f"The numerical result is {res} (+-{err})")
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[4], line 2
1 from math import cos, exp, pi, sin
----> 2 from scipy.integrate import quad
4 # define a function to integrate
5 def f(x):
ModuleNotFoundError: No module named 'scipy'
print(sin(0))
0.0
As you can see, quad returns two values, the first number is the value of the integral and the second is the estimate of the absolute error in the value of the integral.
Exercise: Change the function above to return the sin of the function.
Level: Easy.
# Your code here
# define a function to integrate
def f(x):
return sin((2 * x * pi) + 3.2)
# call quad to integrate f from -2 to 2
res, err = quad(f, -2, 2)
print(f"The numerical result is {res} (+-{err})")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 7
4 return sin((2 * x * pi) + 3.2)
6 # call quad to integrate f from -2 to 2
----> 7 res, err = quad(f, -2, 2)
9 print(f"The numerical result is {res} (+-{err})")
NameError: name 'quad' is not defined
Exercise: Change the lower and upper limits of the integral. Try out a few different values (at least 5).
Level: Easy.
# Your code here
Exercise: Store the numerical values you calculated before into a data structure that you prefer. Plot the values using matplotlib.
Level: Medium.
# Your code here
You can also pass the function f using a lambda function.
f = lambda x: exp(-x**2)
integral = scipy.integrate.quad(f, 0, 1)
integral
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 2
1 f = lambda x: exp(-x**2)
----> 2 integral = scipy.integrate.quad(f, 0, 1)
3 integral
NameError: name 'scipy' is not defined
Equations#
Another package, SymPy can be used also to solve equations, let’s see how that works in a practical example.
In this example, we will calculate the experimental gas (air) hold-up and the theoretical gas hold-up for the conditions applicable in a bubble column experiment.
The gas hold-up is calculated from the equation given below. \begin{equation} \frac{ɛ}{(1-ɛ)^4}= 0.2 \cdot Bo^{1/8}\cdot Ga^{1/12}\cdot Fr \end{equation} In which \(Bo = gD^2\rho/\sigma\), \(Ga = gD^3\rho^2/\eta^2\), and \(Fr = v_G / \sqrt{gD}\).
import numpy as np
from sympy import symbols, solve
g = 9.81 # m/s^2 - gravitational constant
rho = 1000 # kg/m^3 - liquid density
sigma = 0.072 # kg/s^2 - liquid surface tension
eta = 0.00100160 # kg/s*m - dynamic viscosity of water at 20 degC
D = 0.592 # m - Internal diameter of column tube
S = np.pi * (D/2)**2 # m^2 - Cross-sectional area of column
Q = np.array([300, 600]) * 0.001 * 1/60 # m^3/s - volume flowrate
def superficialVelocity(Q, S):
return Q/S
for air_flow in Q:
vG = superficialVelocity(air_flow, S) # m/s
print("vG (m/s)", vG)
Bo = g*D**2*rho/sigma
print("Bo:",Bo)
Ga = g*D**3*rho**2/eta**2
print("Ga", Ga)
Fr = vG/np.sqrt(g*D)
print("Fr", Fr)
epsilon = symbols('epsilon')
expr = 0.2*Bo**(1/8)*Ga**(1/12)*Fr-epsilon/((1-epsilon)**4)
sol = solve(expr) # this is how to use the sympy solve package
print(sol)
print('Simple correlation', 0.6*vG**0.7 )
print()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[10], line 1
----> 1 import numpy as np
2 from sympy import symbols, solve
4 g = 9.81 # m/s^2 - gravitational constant
ModuleNotFoundError: No module named 'numpy'
Solving Systems of Differential Equations#
We will use scipy.integrate.solve_ivp.
Please take some time to read the documentation and get familiar with it.
Note that SciPy has also another package for solving ODEs: scipy.integrate.odeint.
You can read more on the comparison of the two packages in this blog post here. In this course, it is recommended to use scipy. integrate.solve_ivp.
from scipy.integrate import solve_ivp, odeint
def exponential_decay(t, y):
return -0.5 * y
sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
print(sol.t)
print()
print(sol.y)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[11], line 1
----> 1 from scipy.integrate import solve_ivp, odeint
2 def exponential_decay(t, y):
3 return -0.5 * y
ModuleNotFoundError: No module named 'scipy'
import numpy as np
import matplotlib.pyplot as plt
def lotkavolterra(t, z, a, b, c, d):
x, y = z
return [a*x - b*x*y, -c*y + d*x*y]
sol = solve_ivp(lotkavolterra, [0, 15], [10, 5],
args=(1.5, 1, 3, 1), # parameter values a=1.5, b=1, c=3 and d=1 are passed with the args argument
dense_output=True) # dense_output=True computes a continuous solution
#generating t and solving the system of ODEs
t = np.linspace(0, 15, 300)
z1 = sol.sol(t)
# plotting the solution
plt.plot(t, z1.T)
plt.xlabel('t')
plt.legend(['x', 'y'], shadow=True)
plt.title('Lotka-Volterra System')
plt.show()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[12], line 1
----> 1 import numpy as np
2 import matplotlib.pyplot as plt
4 def lotkavolterra(t, z, a, b, c, d):
ModuleNotFoundError: No module named 'numpy'
round(467293847.67899, 2)
467293847.68
Exercise: Write your own function and use the solve_ivp package. Plot the values using matplotlib.
Level: Medium.
# Your code here
Linear Algebra#
For more details, check the documentation page here. We will follow the example of the least-squares solution to equation Ax = b.
To get more information of Least squares.
from scipy.linalg import lstsq
import matplotlib.pyplot as plt
x = np.array([1, 2.5, 3.5, 4, 5, 7, 8.5])
y = np.array([0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6])
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[15], line 1
----> 1 from scipy.linalg import lstsq
2 import matplotlib.pyplot as plt
4 x = np.array([1, 2.5, 3.5, 4, 5, 7, 8.5])
ModuleNotFoundError: No module named 'scipy'
# We want to fit a quadratic polynomial of the form y = a + b*x**2 to this data.
# We first form the “design matrix” M, with a constant column of 1s and a column containing x**2
M = x[:, np.newaxis]**[0, 2]
M
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 3
1 # We want to fit a quadratic polynomial of the form y = a + b*x**2 to this data.
2 # We first form the “design matrix” M, with a constant column of 1s and a column containing x**2
----> 3 M = x[:, np.newaxis]**[0, 2]
4 M
NameError: name 'x' is not defined
Exercise: What does np.newaxis do?
Level: Easy.
# Explanation
We want to find the least-squares solution to M.dot(p) = y, where p is a vector with length 2 that holds the parameters a and b.
p, res, rnk, s = lstsq(M, y)
print('Least-squares solution', p)
print('Residues', res)
print('Effective rank of M', rnk)
print('Singular values of M', s)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 p, res, rnk, s = lstsq(M, y)
2 print('Least-squares solution', p)
3 print('Residues', res)
NameError: name 'lstsq' is not defined
plt.plot(x, y, 'o', label='data')
xx = np.linspace(0, 9, 101)
yy = p[0] + p[1]*xx**2
plt.plot(xx, yy, label='least squares fit, $y = a + bx^2$')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(framealpha=1, shadow=True)
plt.grid(alpha=0.25)
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 plt.plot(x, y, 'o', label='data')
2 xx = np.linspace(0, 9, 101)
3 yy = p[0] + p[1]*xx**2
NameError: name 'plt' is not defined