Variables, loops, lists, and arrays

Getting access to Python

Simple mathematical calculations can be done in plain Python, but for more advanced scientific computing, as we do here, several add-on packages are needed. Getting all software correctly installed used to be quite challenging, but today there are several easy-to-use approaches to get access to Python.

Mac and Windows

We recommend to download and install Anaconda, which is a free distribution of Python that comes with most of the packages you need for advanced scientific computing.

Ubuntu

Debian-based Linux systems, such as Ubuntu, can also use Anaconda, but it is more common to use the apt-get command-line tool or the Ubuntu installer to install a set of Debian packages. Here is a list of the packages you need to install for this introduction:

Terminal> sudo apt-get install python-pip python-setuptools \ 
  python-scipy python-sympy python-cython python-matplotlib \ 
  python-dev python-profiler pydb spyder imagemagick gedit vim \ 
  emacs python-mode git mercurial lib-avtools gnome-terminal

In addition, run

Terminal> pip install nose
Terminal> pip install pytest
Terminal> pip install ipython --upgrade
Terminal> pip install tornado --upgrade
Terminal> pip install pyzmq --upgrade

Web access

You can also access Python directly through a web browser without having it installed on your local computer. We refer to the document How to access Python for doing scientific computing for more details on such tools and also installations based on Anaconda and Ubuntu.

Mathematical example

We shall use a famous mathematical formula in our forthcoming programming examples. The formula tells how long distance \( s \) an object has moved in a time interval \( [0,t] \) if it starts out with a velocity \( v_0 \) and undergoes constant acceleration \( a \): $$ \begin{equation} s = v_0t + \frac{1}{2}at^2\thinspace . \tag{1} \end{equation} $$ We ma y view \( s \) as a function of \( t \): \( s(t) \), and also include the parameters in the notation: \( s(t;v_0,a) \).

A program for evaluating a formula

Here is a Python program for computing \( s \), given \( t=0.5 \), \( v_0=2 \), and \( a=0.2 \):

t = 0.5
v0 = 2
a = 0.2
s = v0*t + 0.5*a*t**2
print(s)
print('s=%g' % s)
print('s\t = \t %.3f' % s)

The program is pure text and must be placed in a pure text file using a text editor. Popular text editors are Gedit, Nano, Emacs, and Vim on Linux, TextWrangler on Mac OS X, and Notepad++ on Windows. Save the text to a program file whose name ends in .py, say distance.py.

The program is run in a terminal window (Command Prompt on Windows, Terminal application on Mac OS X, xterm or gnome-terminal on Linux):

Terminal> python distance.py
1.025

The result of the print statement is the number 1.025 in the terminal window.

As an alternative to writing programs in a text editor and executing them in a terminal window, you may use the Spyder graphical interface which gives a more Matlab-style environment to work in. Anaconda installations come with Spyder.

The code snippet above first contains four assignment statements where we assign numbers or the results of arithmetic expressions to variables. The variables have names coinciding with the mathematical notation used in the formula: t, v0, a, and s. You may think of variables in this programs just as variables in mathematics.

Formatted output with text and numbers

For any object s in Python, print s will (normally) print the contents of s. However, sometimes we want to combine the content of an object with some text. Say we want to print s=1.025 rather than just the number. This is easily accomplished using printf syntax:

print('s=%g' % s)

The output is specified as a string, enclosed in single or double quotes. Inside the string, there can be "slots" for numbers (or other objects), indicated by a percentage sign followed by a specification of kind of data that will be inserted at this place. In the string above, there is one such slot, %g, where the g implies a real number written as compactly as possible.

It is easy to control the number of decimals using printf syntax. Printing out s=1.03, i.e., s with two decimals, is done by

print('s=%.2f' % s)

where the f signifies a decimal number and the preceding .2 means 2 decimals. Scientific notation, as in s=1.03E+00 (\( 1.03\cdot 10^0 \)), is specified as %.2E (2 decimals).

The printf syntax is available in numerous programming languages. Python also offers a related variant, called format string syntax:

print('s={s:.2f}'.format(s=s))

While loops

Suppose we want to make a table with two columns, one with \( t \) values and one with the corresponding \( s \) values. Now we have to repeat a lot of calculations with the formula (1). This is easily done with a loop. There are two types of loops in Python: while loops and for loops.

Let the \( t \) values go from 0 to 2 in increments of 0.1. The following program applies a while loop:

v0 = 2
a = 0.2
dt = 0.1  # Increment
t = 0     # Start value
while t <= 2:
    s = v0*t + 0.5*a*t**2
    print(t, s)
    t = t + dt

The result of running this program in a terminal window is

Terminal> python while.py
0 0.0
0.1 0.201
0.2 0.404
0.3 0.609
0.4 0.816
0.5 1.025
0.6 1.236
0.7 1.449
0.8 1.664
0.9 1.881
1.0 2.1
1.1 2.321
1.2 2.544
1.3 2.769
1.4 2.996
1.5 3.225
1.6 3.456
1.7 3.689
1.8 3.924
1.9 4.161

So, how do we interpret the contents of this program? First we initialize four variables: v0, a, dt, and t. Everything after # on a line is a comment and does not affect what happens in the program, but is meant to be of help for a human reading the program. Then comes the while loop:

while condition:
    <intented statement>
    <intented statement>
    <intented statement>

Observe the colon at the end of the while line. The set of indented statements are repeated as long as the expression or variable condition evaluates to True. In the present case, the condition is the boolean expression t <= 2, so as long as the value is less than or equal to 2, t <= 2 evaluates to True, otherwise it evaluates to False.

The Python Online Tutor is a great tool to examine the program flow. Consider this little loop run in the Python Online Tutor (view in Chrome):

The Python Online Tutor executes each statement and displays the contents of variables such that you can track the program flow and the evolution of variables.

So, why is not a printed for 1.2? That is, why does a == 1.2 evaluate to True when a is (expected to be) 2? We must look at the accuracy of a to investigate this question and write it out with 16 decimals in scientific notation (printf specification %.16E):

The problem is that da is not exactly 0.4, but contains a small round-off error. Doing a = a + da then results in a slightly inaccurate a. When the exact a should reach 1.2, the a in the program has a an error and equals 1.2000000000000002. Obviously the test for exact equality, a == 1.2, becomes False, and the loop is terminated.

a = 1.2
b = 0.4 + 0.4 + 0.4
boolean_condition1 = a == b              # False
print(boolean_condition1)
tol = 1E-14
boolean_condition2 = abs(a - b) < tol    # True
print(boolean_condition2)

The Python Online Tutor is ideal for demonstrating program flow and contents of variables in small and simple programs. However, you should use a real debugger instead when searching for errors in real programs.

Lists

The table created in the previous section has two columns of data. We could store all the numbers in each column in a list object. A list is just a collection of objects, here numbers, in a given sequence. For example,

L = [-1, 1, 8.0]

is a list of three numbers. Lists are enclosed in square brackets and may contain any type of objects separated by commas. Here we mix a filename (string), a real number, and an integer:

L = ['mydata.txt', 3.14, 10]

The different list elements can be reached via indexing: L[0] is the first element, L[1] is the second, and L[-1] is the last element. Here is an interactive Python shell where we can write Python statements and examine the contents of variables as we perform various operations:

L = ['mydata.txt', 3.14, 10]
print(L[0])
print(L[1])
del(L[0])  # delete the first element
print(L)
print(len(L))  # length of L
L.append(-1)  # add -1 at the end of the list
print(L)

Let us store the numbers in the previous table in lists, one for each column. We can start with empty lists [] and use append to add a new element to the lists in each pass of the while loop. Thereafter, we can run a new while loop and print the contents of the lists in a nice, tabular fashion:

v0 = 2
a = 0.2
dt = 0.1  # Increment
t = 0
t_values = []
s_values = []
while t <= 2.1:
    s = v0*t + 0.5*a*t**2
    t_values.append(t)
    s_values.append(s)
    t = t + dt
print(s_values)  # Just take a look at a created list
print(t_values)

# Print a nicely formatted table
i = 0
while i <= len(t_values)-1:
    print('%.18f  %.4f' % (t_values[i], s_values[i]))
    i += 1   # Compact form for i = i + 1

The output looks like

[0.0, 0.201, 0.404, 0.6090000000000001, 0.8160000000000001,
 1.025, 1.236, 1.4489999999999998, 1.664, 1.8809999999999998,
 2.0999999999999996, 2.3209999999999997, ...]
0.00  0.0000
0.10  0.2010
0.20  0.4040
0.30  0.6090
0.40  0.8160
0.50  1.0250
0.60  1.2360
0.70  1.4490
...

Note that print s_values here leads to output with many decimals and small round-off errors. To get complete control of the formatting of real numbers in the table, we use the printf syntax.

Lists come with a lot of functionality. See the Python Tutorial for many more examples.

For loops

A for loop is used for visiting elements in a list, one by one:

>>> L = [1, 4, 8, 9]
>>> for e in L:
...     print(e)
...
1
4
8
9

The variable e is successively set equal to the elements in the list, in the right order. Note the colon at the end of the for line. The statements to be executed in each pass of the loop (here only print e) must be indented. When e is set equal to the last element and all indented statements are executed, the loop is over, and the program flow continues with the next statement that is not indented. Try the following code out in the Python Online Tutor:

A for loop over the valid indices in a list is created by

for i in range(len(somelist)):
    # Work with somelist[i]

The range function returns a list of integers: range(a, b, s) returns the integers a, a+s, a+2*s, ... up to but not including b. Just writing range(b) implies a=0 and s=1, so range(len(somelist)) returns [0, 1, 2].

# Need loop over 0, 0.1, 0.2, ..., 1
values = []
for i in range(11):
    values.append(i*0.1)

# Shorter version using list comprehension (same as the loop above)
values = [i*0.1 for i in range(11)]

for value in values:
    # work with value

We can now rewrite our program that used lists and while loops to use for loops instead:

v0 = 2
a = 0.2
dt = 0.1  # Increment
t_values = []
s_values = []
n = int(round(2/dt)) + 1  # No of t values
for i in range(n):
    t = i*dt
    s = v0*t + 0.5*a*t**2
    t_values.append(t)
    s_values.append(s)
print(s_values)  # Just take a look at a created list

# Make nicely formatted table
for t, s in zip(t_values, s_values):
    print('%.2f  %.4f' % (t, s))

# Alternative implementation
for i in range(len(t_values)):
    print('%.2f  %.4f' % (t_values[i], s_values[i]))

Observe that we have here used a slightly different technique for computing the \( t \) values inside the first loop: now we set \( t \) as \( i\Delta t \), where \( \Delta t \) (dt in the code) is the increment (0.1) between each \( t \) value. The computation of n, the number of \( t \) values, makes use of round to make a correct mathematical rounding to the nearest integer (and int makes an integer object out of the rounded real number). (In an interval \( [a,b] \) divided into subintervals of equal length \( \Delta t \), there will be \( 1 + (b-a)/\Delta t \) points in total.)

Running through multiple lists simultaneously is done with the zip construction:

for e1, e2, e3, ... in zip(list1, list2, list3, ...):

One may instead create a for loop over all the legal index values instead and index each array,

for i in range(len(list1)):
    e1 = list1[i]
    e2 = list2[i]
    ...

Arrays

Lists are useful for collecting a set of numbers or other objects in a single variable. Arrays are much like lists, but tailored for collection of numbers. The primary advantage of arrays is that you can use them very efficiently and conveniently in mathematical computations, but the downside is that an array has (in practice) a fixed length and all elements must be of the same type. This is usually no important restriction in scientific computations.

Much of Python's functionality are coded in modules. To use such functionality, one must import the relevant modules. Arrays, for example, are available in the numpy module, which must be imported. Functions or variables in the module must be prefixed by numpy. This will be demonstrated below. Here is an example involving basic operations with arrays:

L = [1, 4, 10.0]    # List of numbers
import numpy
a = numpy.array(L)  # Make corresponding array
print(a)
print(a[1])
print(a.dtype)       # Data type of an element
b = 2*a + 1
print(b)

Note that all elements in the a array are of float type (because one element in L was float). Arithmetic expressions such as 2*a+1 work with a as array, but not as list. In fact, we can pass arrays to mathematical functions:

>>> c = numpy.log(a)
>>> print(c)
[ 0.          1.38629436  2.30258509]

The numpy module has a lot of very useful utilities. To create \( n+1 \) uniformly distributed coordinates in an interval \( [a,b] \), stored in an array, one can use linspace:

t = numpy.linspace(a, b, n+1)

This construction makes it easy to create arrays for the \( t \) and \( s \) values in our tables:

import numpy
v0 = 2
a = 0.2
dt = 0.1  # Increment
n = int(round(2/dt)) + 1  # No of t values

t_values = numpy.linspace(0, 2, n+1)
s_values = v0*t_values + 0.5*a*t_values**2

# Make nicely formatted table
for t, s in zip(t_values, s_values):
    print('%.2f  %.4f' % (t, s))

Mathematical functions

Python offers access to all standard mathematical functions such as \( \sin x \), \( \cos x \), \( \tan x \), \( \sinh x \), \( \cosh x \), \( \tanh x \), all their inverses (called asin(x), asinh(x), and so forth), \( e^x \) (exp(x)), \( \ln x \) (log(x)), and \( x! \) (factorial(x)). However, one has to import a module to get access to these functions. For scalars (single numbers) the relevant module is math:

>>> import math
>>> print(math.sin(math.pi))
1.2246467991473532e-16

which shows that the sine function is only approximate (to 16 digits). Many prefer to write mathematical expressions without the math prefix:

from math import sin, pi
print(sin(pi))

# Or import everything from math
from math import *
print(sin(pi), log(e), tanh(0.5))

The numpy module contains sine, cosine, and other mathematical functions that work on scalars as well as arrays.

from numpy import *
x = linspace(0, 1, 101)
y = exp(-x)*sin(pi*x)

import numpy as np
x = np.linspace(0, 1, 101)
y = np.exp(-x)*np.sin(np.pi*x)

import numpy as np
from numpy import sin, exp, pi
x = np.linspace(0, 1, 101)
y = exp(-x)*sin(pi*x)

Plotting

We can easily make a graph of a function \( s(t) \) using the module matplotlib. The technique is to compute an array of \( t \) values and a corresponding array of function values \( s(t) \). Plotting programs will draw straight lines between the points on the curve, so a sufficient number of points are needed to give the impression of a smooth curve. Our \( s(t) \) function is plotted by the following code:

import numpy as np
import matplotlib.pyplot as plt

v0 = 0.2
a = 2
n = 21  # No of t values for plotting

t = np.linspace(0, 2, n+1)
s = v0*t + 0.5*a*t**2

plt.plot(t, s)
plt.savefig('myplot.png')
plt.show()

The plotfile myplot.png looks like

Matlab users may prefer to do

from numpy import *
from matplotlib.pyplot import *

such that they can use linspace and plot without any prefix, just as in Matlab.

Two curves can easily be plotted, this time also with labels on the axis and a box with legends such that we can distinguish the two curves:

import numpy as np
import matplotlib.pyplot as plt

v0 = 0.2
a = 2
n = 21  # No of t values for plotting

t = np.linspace(0, 2, n+1)
s = v0*t + 0.5*a*t**2

plt.plot(t, s)
plt.savefig('myplot.png')
plt.show()

• «
• 1
• 2
• 3
• 4
• 5
• 6
• »