Numpy

The NumPy library

NumPy is mainly used to create and edit arrays. An array is a data structure similar to a list, with the difference that it can contain only one type of object. For example you can have an array of integers, an array of floats, an array of strings etc, however you can't have an array that contains two datatypes at the same time. But then why use arrays instead of lists? Because they are quicker, they have more functionalities, and they better represent mathematical concepts like vectors and matrices.

 

So essentially what NumPy does is provide the user with a powerful object called ndarray (but generally referred to simply as array) that expands the functionalities of Python quite a bit.

Loading the library

In Python, when you want to use an external library like NumPy, you have to load it before you can access its functionalities. Usually you would also need to install it first, but luckily Jupyter has all the libraries that we need pre-installed. Open Jupyter, write the following code in an empty cell and run it:

 

import numpy as np
a=np.zeros((4))
print(a)

[0. 0. 0. 0.]

 

To load a library, we use the keyword import. we can then specify an abbreviation for that library that will be used throughout the code, using the keyword as. In line 1 we import numpy, and we assign the nickname "np" to it. This way to access the library's functionalities it will be enough to write "np" followed by whatever function you need. In this example we create an array of zeros using NumPy's built in method zeros, that takes as input a tuple that specifies the dimensions of the array to be created. We want a one-dimensional array of four elements, so the tuple will just contain 4.

 

From this point on, in the code examples written in this tutorial the line import numpy as np will be omitted. They will all work, as long as you have already loaded NumPy in another cell.

Creating arrays

You can create NumPy arrays from lists and tuples with the function array():

 

a=np.array([1,2,3,4])
b=np.array([[1,2],[3.5, 4.2]])
print("a={0}".format(a))
print("b={0}".format(b))

a=[1 2 3 4]
b=[[1. 2. ]
[3.5 4.2]]

 In line 1 we pass the array function a lists of integers, so NumPy creates an array of integers. In line 2 we pass a multidimensional list that contains both integers and floats, and as you can see NumPy automatically converts the int values to float, because an array can contain only one type of data.

 

If you know the dimensions of an array that you are going to need, but you want to populate its positions later, you can create arrays of zeros and ones:

a=np.zeros((3))
b=np.ones((2,3), dtype=int)
print("a={0}".format(a))
print("b={0}".format(b))

a=[0. 0. 0.]
b=[[1 1 1]
[1 1 1]]

 The functions ones() and zeros() take as inputs a tuple that specifies the dimensions of the axes of the array. In line 2 we pass the tuple (2,3), in order to crate an array of ones with two lines and three columns. Both zeros() and ones() can take the additional parameter dtype, which is used to specify the dtype of the elements in the array. When we create b we specify int as the dtype, so b will be composed of integer values.

 

To create an array that contains a sequence of numbers that goes from a starting value to a finishing value, you can use the linspace() function:

a=np.linspace(1,10,6)
print("a={0}".format(a))

a=[ 1. 2.8 4.6 6.4 8.2 10. ]

linspace takes as arguments the starting point, the ending point and the number of elements of the sequence that you want to create. In the previous code cell we create an array that goes from 1 to 10, composed by 6 evenly spaced values.

Retrieving an array's dimension and dtype

When you need to retrieve information about an existing array, NumPy offers the following methods:

  • ndarray.ndim: returns the number of axes of the array
  • ndarray.shape: returns the dimensions of the array as a tuple. Essentially it returns the number of elements in each axis
  • ndarray.size: the total number of elements in the array
  • ndarray.dtype: returns the dtype of the elements of the array

Let's see these methods in action:

a=np.array([[1,2],[3,4]])
print(a.ndim)
print(a.shape)
print(a.size)
print(a.dtype)

2
(2, 2)
4
int32

As you can see, a is composed by int32-type objects. This is one of the various int-type objects provided by NumPy, and it represents a 32-bit long integer. You generally don't need to worry about which type of integer or float NumPy is using, as long as you recognize that it is  indeed an int or float.

Access and modify the elements of an array

NumPy arrays are indexed starting from zero like lists, however the syntax used to access a certain value within an array is a little bit different. Let's take a look at code snippet:

a=np.array([2.45, 3.75])
a[0]=16
b=np.array([[0,1],[2,3]])
print(a[0])
print(b[1,1])

16.0
3

To access an element of an array we use the square brackets like we did for lists. Notice that when we need to retrieve a value from a two-dimensional array like b, we specify the position of the element using only one set of square brackets, in this manner:

[i,j]

If you remember, for lists the same operation was made using this syntax instead:

[i][j]

Be careful not to confuse the two!

As we did for lists, we can change the value stored at a certain position of the array, like we see in line 2.

Mathematical operations between arrays

NumPy performs mathematical operations between arrays in a element-wise fashion. This means that the arrays must all have the same number of elements in them, otherwise you will get an error.Let's see some of the possible operations that can be performed:

a=np.array([1,2,3,4]) b=np.array([2,4,8,16]) print(a+b)
print(a-b)
print(a*b)
print(a/b)
print(a**2)
print(np.sqrt(b))

[ 3 6 11 20]\
[ -1 -2 -5 -12]
[ 2 8 24 64]
[0.5 0.5 0.375 0.25 ]
[ 1 4 9 16]
[1.41421356 2. 2.82842712 4. ]

 If you have used Excel before, you may have noticed that for example the operation in line 3 is akin to summing two columns. These kinds of operations between arrays are very common, and indeed very useful. but what if we wanted to perform the matrix product instead of the element-wise multiplication? Then we would use @ instead of * as the multiplication symbol:

a=np.array([[1,2],[3,4]])
b=np.array([[2,4],[8,16]])
print(a@b)

[[18 36]
[38 76]]

Array slicing

Array slicing is used to perform operations only on certain parts of a given array. The next code cell shows the most commonly used slicing operations for one-dimensional arrays: 

a=np.linspace(0,10,11) print(a)
print(a[2:5])
print(a[1:8:2])

[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
[2. 3. 4.]
[1. 3. 5. 7.]

 As you can see we can use semicolons to print only the elements between a starting point and an ending point an an array. Moreover, in line 3 we see how to select elements at certain intervals: the syntax 1:8:2 means"take every second element of a between the positions 1 and 8".

 

For multidimensional arrays the syntax is the same, except each axis of the array can be sliced independently. To understand what this means, let's go over code the next piece of code:

mat=np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
print(mat[:,2])
print(mat[3,:])
print(mat[:,2:4])

[ 3 7 11 15]
[13 14 15 16]
[[} 3 4]
[ 7 8]
[11 12]
[15 16]]

 In line 2 we select all the lines by using the semicolon, and only the third column by passing 2 as the positional index for the second axis. Line 3 is similar, except it selects the fourth line of the matrix. The print order in line 3 displays the third and fourth columns instead.

 

NOTE:  You may have noticed that when slicing using the semicolon, NumPy cuts the array from the starting position to the ending position excluded. This is why line 4 does not give an error, even tough the index 4 is out of bounds for an array with 4 columns. (remember that arrays are indexed from 0!)

Array union and concatenation

The functions hstack and vstack can be used to merge two or more arrays:

a=np.array([1,1])
b=np.array([2,2])
c=np.array([3,3])
print(np.hstack((a,b,c)))
print(np.vstack((a,b,c)))

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

 while append and insert can be used to add new elements inside the array:

a=np.array([1,2,3])
print(a)
a=np.append(a,4)
print(a)
a=np.insert(a,2,7)
print(a)

[1 2 3]
[1 2 3 4]
[1 2 7 3 4]

As you can see, append adds an element at the end of an array, and insert adds an element at a certain position. The second argument of insert can be an iterable object (like a list or another array) that specifies more than one position where to add the element specified by the third argument. Both append and insert take the optional argument axis, used to specify the axis along which the operation must be performed.

 

NumPy offers many other functionalities on top of the ones explained in the previous sections. So many, in fact, that they would require a separate book to be explained in their entirety. In the next chapters we will introduce some of them, however we will focus our attention only on the ones that are most useful and relevant to the practical goals that we have.

EXAMPLE: calculating the internal forces of a beam using NumPy

We can finally move on to the first practical example of this book. In this section many of the concepts that have been explained before will come together, hopefully giving you a hint of their possible practical applications. Since this is the first time we use Python (and Jupyter) in a more structured way, the following example will be very simple from an engineering point of view, and yet very useful to understand the functionalities explained in the previous chapters.

Problem definition

The goal of this example is to calculate the bending moment and shear of simply supported beam using NumPy. The beam is shown in the next figure:

Where: 

$$q=20\;kN/m$$

$$l=5\;m$$

Input the problem's quantities in Python

The first step to solve many problems using a computer is simply to store the various given quantities in variables, so that they can be used throughout the code. But before we can do that, we have to load NumPy, so we can use its functionalities: 

import numpy as np

 It's common practice to put all your library-loading code lines in the first cell of the notebook, to keep it separated form the rest. Remember to run this cell at least once before you run other cells that contain NumPy functionalities, otherwise you will get an error!

 

In the next cell, we store the length of the beam and the value of the distributed load in two variables:

l=5 #m q=20 #kN/m

 Notice how each quantity has its own unit of measure written as a comment next to it. This may seem a bit superfluous given the simplicity of the problem, but it will most certainly become a necessary habit when tackling more difficult applications. It's common practice to keep code inputs in their own cell as well, so that you know where to find them in case you need to modify them. From this point on, however, it's up to you to decide how to organize your code. Keep in mind that the easier it is to understand your code, the easier it is to find bugs.

Calculating the bending moment and shear

The beam's ground reactions and coordinate system are shown in the next figure:

The equilibrium equations are:

$$H_A=0$$

$$V_A+V_b-ql=0$$

$$V_Bl-\cfrac{ql^2}{2}=0$$

solving the system we obtain

$H_A=0$        $V_A=\cfrac{ql}{2}$        $V_B=\cfrac{ql}{2}$

so the bending moment \textbf{M} and the shear \textbf{V} of the beam will be equal to

$$\mathbf{M}=V_Ax-\cfrac{qx^2}{2}=\cfrac{q}{2}\left(lx-x^2\right)$$

$$\mathbf{V}=V_A-qx=q\left(\cfrac{l}{2}-x\right)$$

Our goal is to translate all of these equations in Python language, and obtain two vectors named M and V that contain the value of the bending moment and the shear of the beam for x that goes from A to B. All of this is done in the next code snippet:

x=np.linspace(0,l,20) M=q/2*(l*x-x**2) V=q*(l/2-x) print(M)

[ 0.         12.46537396 23.54570637 33.24099723 41.55124654 48.47645429  54.0166205  58.17174515 60.94182825 62.32686981 62.32686981 60.94182825  58.17174515 54.0166205  48.47645429 41.55124654 33.24099723 23.54570637  12.46537396  0.        ]

 In line 1 we create the x coordinate using the NumPy function linspace. As you can see we pass l as the ending value and 20 as the number of steps, meaning that NumPy will create a vector of evenly spaced numbers that go from 0 to l, which in this case is equal to 5. From this point on it's just a matter of copying the formulas: NumPy will perform all the operations element-wise, so M and V will be arrays composed of 20 elements each, corresponding to every coordinate value contained in x.

 

If you are familiar with Excel, you can think of x as a column of values to which you apply the formulas in lines 3 and 4, in order to obtain two new columns containing the values of M and V.

 

In future we will learn how to use matplotlib to plot various graphs and diagrams. For now, just copy and paste the following code in an empty cell to display M and V:

which is still pretty basic, but it's enough to understand that the values of M and V are correct.

Go to the next link:

Sympy tutorial