Chemical & Biological Engineering CBE Computational Turbulent Reacting Flow CTRFlow

Publications Calendar Research Teaching Zoom LDS

Publications Calendar Research Teaching Zoom LDS

Python Arrays

Import Libraries

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Array basics were covered earlier

• Review those if needed
• See this page for more information on built-in Python Lists.

Numpy arrays

• Python lists are nice, but for numerical work, numpy arrays are much better.
• 90% of the work we will do in science, math, and engineering, is best done using numpy arrays.
• All of the above regarding Python lists is the same, except
  • Each element type needs to be the same.
  • The way we create arrays is different.

Some good links:

• Numpy for Matlab users
  • This too
• Tentative numpy tutorial
• Numpy quickstart tutorial
• Linear algebra
• Cookbook with good recipes

Length of a list

len function

a = np.array([1.1, 2.2, 3.3, 4.4])

len(a)

$\displaystyle 4$

Creating arrays

• array (we've already seen this)
• linspace, arange
• empty, zeros, ones, full

# linspace

np.linspace(0,5,11)

array([0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. ])

# arange

np.arange(0,20,2.2)

array([ 0. ,  2.2,  4.4,  6.6,  8.8, 11. , 13.2, 15.4, 17.6, 19.8])

# empty

np.empty(5)         # creates space, but arbitrary values (not set, faster)

array([ 1. ,  2. ,  2.5,  5. , 10. ])

# zeros

np.zeros(5)         # array of 5 zeros

array([0., 0., 0., 0., 0.])

# ones

np.ones(5)           # array of 5 ones

array([1., 1., 1., 1., 1.])

# full

np.full(5, 2.2)      # array of 5 2.2's

array([2.2, 2.2, 2.2, 2.2, 2.2])

Access elements

• use square brackets with index number
• array indexing starts at zero
• for 2D arrays (matrix) x[i,j] is row i, column j
• can access backwards0
  • x[-1] is the last element
  • x[-2] is the second to the last element, etc.
  • See the figure below. If we wrap the array into a circle. Then element 0 is at the top. We index by +1 in moving clockwise, and index by -1 in moving counter-clockwise.

x = np.arange(3,21,3)               # x is         [3 6 9 12 15 18]
print(x)

# Access elements with index values

print( x[[1,3,5]] )

# or use an "index array"

i = [1,3,5]
print( x[i] )

[ 3  6  9 12 15 18]
[ 6 12 18]
[ 6 12 18]

Array slices

• Can access multiple elements using syntax: a[s:e]
  • s is the start (optional)
  • e is the end (optional)
  • (stops short of the end)

x = np.arange(0,11,1.5)       

print( x[:] )                      # get all elements
print( x[2:] )                     # get indices 2 to the end
print( x[:-2] )                    # get elements start to one before 2nd to last (stops short)
print( x[2:5] )                    # get indices 2 to 4 (stops short of 5; prints 5-2=3 elements)

[ 0.   1.5  3.   4.5  6.   7.5  9.  10.5]
[ 3.   4.5  6.   7.5  9.  10.5]
[0.  1.5 3.  4.5 6.  7.5]
[3.  4.5 6. ]

Example

We have an array of times $t$, and we want the array of differences $dt$ between times.

• How big is $t$?
• How big will $dt$ be?
• Create $dt$ using a loop.
• Create $dt$ using slice-indexing.
  • Note, for arrays a and b, a-b is done element-wise.

t = np.array([ 
    0.01136646,  0.01655847,  0.02548247,  0.07698612,  0.17964401,
    0.19813924,  0.22419994,  0.28622096,  0.28884997,  0.3627257 ,
    0.42555072,  0.44965767,  0.47257237,  0.59335946,  0.73842352,
    0.85900369,  0.88851226,  0.89153229,  0.89355469,  0.97705246,
    ])

# loop

dt = np.empty(len(t)-1)
for i in range(len(dt)):
    dt[i] = t[i+1] - t[i]
print(dt)

# slice

dt = t[1:] - t[:-1]
print(dt)

[0.00519201 0.008924   0.05150365 0.10265789 0.01849523 0.0260607
 0.06202102 0.00262901 0.07387573 0.06282502 0.02410695 0.0229147
 0.12078709 0.14506406 0.12058017 0.02950857 0.00302003 0.0020224
 0.08349777]
[0.00519201 0.008924   0.05150365 0.10265789 0.01849523 0.0260607
 0.06202102 0.00262901 0.07387573 0.06282502 0.02410695 0.0229147
 0.12078709 0.14506406 0.12058017 0.02950857 0.00302003 0.0020224
 0.08349777]

Array properties

• size and shape

x = np.array([[1,2,3,4],[7,8,9,10]])
print(x)

np.shape(x)            # rows x columns

[[ 1  2  3  4]
 [ 7  8  9 10]]

$\displaystyle \left( 2, \ 4\right)$

np.size(x)             # total number of elements

$\displaystyle 8$

len(x)                 # length of first dimension (rows)

$\displaystyle 2$

x = np.array([1,2,3,4])
len(x)

$\displaystyle 4$

Arithmetic

a = np.arange(9)
a

array([0, 1, 2, 3, 4, 5, 6, 7, 8])

#------- THIS

b = 0.1 * a                 # multiplies each element of a by 0.1 and assigns whole resulting array to b

#------- NOT THIS

b = np.empty(len(a))
for i in range(len(b)):
    b[i] = 0.1*a[i]
    
#-------

print(b)

[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8]

a+b                 # adds element-wise

array([0. , 1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8])

a**2                # squares element-wise

array([ 0,  1,  4,  9, 16, 25, 36, 49, 64])

BE CAREFUL about assignment

x = np.array([1,2,3,4]) 
y = x                    # note, y is now just a reference to x.
y[3] = 44                # if we change an element of y, it will change in x too!
print(x)                 # this is not usually what we want
print(y)

[ 1  2  3 44]
[ 1  2  3 44]

x = np.array([1,2,3,4])
y = x.copy()                 # now x and y are two different objects
y[3] = 44                    #   could also do y = np.array(x)
print(x)
print(y)

[1 2 3 4]
[ 1  2  3 44]

Basic Plots

x = np.linspace(0,20,1000)
y = np.exp(-x/5)*np.sin(x)

plt.plot(x,y);

Problem 1

• Make an array of $x$ from 0 to 5 with 1000 points.
• create the following arrays
  • $2x$
  • $x^2 -3x + 4$
  • $\sin(5x)$
• Plot each array versus $x$ on a single plot.

x = np.linspace(0,5,1000)
y = 2*x
z = x**2 - 3*x + 4
w = np.sin(5*x)

plt.plot(x,y)
plt.plot(x,z)
plt.plot(x,w);

Problem 2

• Ideal gas law $$\rho = \frac{MP}{RT}$$
• M=29 kg/kmol
• R=8314 J/kmol*K
• P=101325 Pa
• Find $\rho$ for $T=300,400,\ldots1000$ K
  • Store $\rho$ and $T$ as arrays
  • Plot $\rho$ versus $T$

M = 29     # kg/kmol
R = 8314   # J/kmol*K
P = 101325 # Pa

# define T array using a numpy array building function

# compute rho using the above formula (don't use a for loop) 

# plot the rho versus T

Problem 3

• Redo Problem 2 by writing a function for $\rho(T)$ as you would have in a previous lecture.
  • Call the function with an array of $T$ and store the result as an array of $\rho$.
  • Plot $\rho$ versus $T$.

Example

def f(T):                 # if T is an array, then 
    return 5*T            # return 5 times the T array, which is an array
    
T = np.array([300, 400, 500])

fiveTimesT = f(T)    # fiveTimesT is np.array([1500, 2000, 2500])

Problem 4

• You have an array of $x$ points $x=0.0, 0.1, \ldots, 50$.
• Create $x$
• Create $y=\sin(x)$
• At interior points, use a central difference approximation to find the derivative $$\left(\frac{dy}{dx}\right)_i \approx \frac{y_{i+1}-y_{i-1}}{2\Delta x},$$ where $\Delta x = x_1-x_0$.
• At edge point use a one-sided difference approximation, $$\left(\frac{dy}{dx}\right)_0 \approx \frac{y_{1}-y_{0}}{\Delta x},$$ $$\left(\frac{dy}{dx}\right)_{n-1} \approx \frac{y_{n-1}-y_{n-2}}{\Delta x}.$$
• Plot $dy/dx$ versus $x$ , and $\cos(x)$ versus $x$.

# y: * (*  *  *  *  *  *) *        # boundary, (interior), boundary
#    0  1  2  3  4  5  6  7
#                     n-2 n-1

#---------- create array x

#---------- create array y

#---------- create array dydx
# First make an empty or unfilled array to store the values: np.empty, or np.zeros
# (the size is the same as x and y since you want the derivative at every point)
# Fill in dydx[0] using the formula given
# Fill in the last value of dydx using the formula given (what "index" for dydx is best to use?)
# Fill in the interior values using an array slice

#---------- plot results

Problem 5

• Solve the rate equation below to $t=5$ s using the Explicit Euler method $y_{i+1} = y_i + \Delta tf(y_i,t_i).$

$$\frac{dy}{dt} = -ky,$$ $$y_0=1,$$ $$k=2.$$

• Use $\Delta t = 0.5$.
• Store $t$ and $y$ in arrays.
• Plot $y$ versus $t$, and $y_0e^{-kt}$ versus $t$ on the same plot.

#---------- think through the approach:
# what variables should be defined?
# what arrays should be defined
# how to fill in dydt (recommend using a for loop)

Problem 6

• The Haaland equation relates the friction factor $f$ in turbulent pipe flow to the Reynolds number $Re$. $$\frac{1}{\sqrt{f}} = -1.8\log_{10}\left[\left(\frac{\epsilon/D}{3.7}\right)^{1.11}+ \frac{6.9}{Re}\right].$$

• Write a function to compute $f(Re, \epsilon/D)$.

• Create an array of 100 points that is uniformly spaced on a log scale, for $1000\le Re\le 1\times 10^8$

• Plot $f$ versus $\log_{10}(Re)$.

• Use $\epsilon/D=0.001$.