Python Man
Python Man
Python Man
Science
Release 0.9.23
David Pine
CONTENTS
Introduction
1.1 Introduction to Python and its use in science . . . . . . . . . . . . . . .
Launching Python
2.1 Installing Python on your computer . . .
2.2 The Canopy window . . . . . . . . . . .
2.3 The Interactive Python Pane . . . . . . .
2.4 Interactive Python as a calculator . . . .
2.5 Python Modules . . . . . . . . . . . . .
2.6 Python functions: a first look . . . . . .
2.7 Variables . . . . . . . . . . . . . . . . .
2.8 Script files and programs . . . . . . . . .
2.9 Importing Modules . . . . . . . . . . . .
2.10 Getting help: documentation in IPython .
2.11 Programming is a detail-oriented activity
2.12 Exercises . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5
5
5
6
11
14
15
17
20
24
25
26
27
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
29
30
31
35
46
48
52
55
55
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3
3
4.2
4.3
4.4
4.5
5
Screen output
File input . .
File output .
Exercises . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
57
62
66
70
Plotting
5.1 An interactive session with pyplot
5.2 Basic plotting . . . . . . . . . . . . .
5.3 Logarithmic plots . . . . . . . . . . .
5.4 More advanced graphical output . . .
5.5 Exercises . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
73
74
75
89
91
96
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
99
. 99
. 105
. 110
. 112
Functions
7.1 User-defined functions . . . . . . .
7.2 Methods and attributes . . . . . . .
7.3 Example: linear least squares fitting
7.4 Exercises . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
115
115
129
130
136
Curve Fitting
141
8.1 Using linear regression for fitting non-linear functions . . . . . . . . . . 141
8.2 Nonlinear fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
8.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
157
158
161
165
170
173
178
A Installing Python
179
A.1 Installing Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
A.2 Testing your installation of Python . . . . . . . . . . . . . . . . . . . . . 180
B IPython Notebooks
183
B.1 An interactive interface . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
ii
B.2
B.3
B.4
B.5
B.6
B.7
B.8
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
183
185
187
187
190
191
191
C Python Resources
193
C.1 Web resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
C.2 Books . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
Index
197
iii
iv
Contents:
CONTENTS
CONTENTS
CHAPTER
ONE
INTRODUCTION
1.1 Introduction to Python and its use in science
This manual is meant to serve as an introduction to the Python programming language
and its use for scientific computing. Its ok if you have never programmed a computer
before. This manual will teach you how to do it from the ground up.
The Python programming language is useful for all kinds of scientific and engineering
tasks. You can use it to analyze and plot data. You can also use it to numerically solve
science and engineering problems that are difficult or even impossible to solve analytically.
While we want to marshall Pythons powers to address scientific problems, you should
know that Python is a general purpose computer language that is widely used to address
all kinds of computing tasks, from web applications to processing financial data on Wall
Street and various scripting tasks for computer system management. Over the past decade
it has been increasingly used by scientists and engineers for numerical computations,
graphics, and as a wrapper for numerical software originally written in other languages,
like Fortran and C.
Python is similar to Matlab and IDL, two other computer languages that are frequently
used in science and engineering applications. Like Matlab and IDL, Python is an interpreted language, meaning you can run your code without having to go through an extra
step of compiling, as required for the C and Fortran programming languages. It is also
a dynamically typed language, meaning you dont have to declare variables and set aside
memory before using them. Dont worry if you dont know exactly what these terms
mean. Their primary significance for you is that you can write Python code, test, and use
it quickly with a minimum of fuss.
One advantage of Python over similar languages like Matlab and IDL is that it is free. It
can be downloaded from the web and is available on all the standard computer platforms,
including Windows, MacOS, and Linux. This also means that you can use Python without
being tethered to the internet, as required for commercial software that is tied to a remote
license server.
Another advantage is Pythons clean and simple syntax, including its implementation of
object oriented programming (which we do not emphasize in this introduction).
An important disadvantage is that Python programs can be slower than compiled languages like C. For large scale simulations and other demanding applications, there can
be a considerable speed penalty in using Python. In these cases, C, C++, or Fortran is
recommended, although intelligent use of Pythons array processing tools contained in
the NumPy module can greatly speed up Python code. Another disadvantage is that compared to Matlab and IDL, Python is less well documented. This stems from the fact that it
is public open source software and thus is dependent on volunteers from the community
of developers and users for documentation. The documentation is freely available on the
web but is scattered among a number of different sites and can be terse. This manual will
acquaint you with the most commonly-used web sites. Search engines like Google can
help you find others.
You are not assumed to have had any previous programming experience. However, the
purpose of this manual isnt to teach you the principles of computer programming; its
to provide a very practical guide to getting started with Python for scientific computing.
Perhaps once you see some of the powerful tasks that you can accomplish with Python,
you will be inspired to study computational science and engineering, as well as computer
programming, in greater depth.
Chapter 1. Introduction
CHAPTER
TWO
LAUNCHING PYTHON
2.1 Installing Python on your computer
If you havent already installed Python on your computer, see Installing Python, which
includes instructions for installing Python on Macs running under MacOSX and on PCs
running under Windows.
Once you have installed Python, find the Canopy icon on your computer and launch
the application. Wait for the Canopy welcome screen to appear, and then click on the
Editor icon. The Canopy window should appear, like the one shown below. This is the
window you will generally use to work with Python.
This means that Canopy is running a particular version or shell of interactive Python
called IPython. The IPython shell has been specifically designed for scientific and engineering use. The standard Python interactive shell uses the prompt >>>. You can pretty
much do everything you want to do with either shell, but we will be using the IPython
shell as we want to take advantage of some of its special features for scientific computing.
By typing short commands at the prompt, IPython can be used to perform various system
tasks, such as running programs and creating and moving files around on your computer.
This is a different kind of computer interface than the icon-based interface (or graphical
user interface GUI) that you usually use to communicate with your computer. While
it may seem more cumbersome for some tasks, it can be more powerful for other tasks,
particularly those associated with programming.
Before getting started, we point out that like most modern computer languages, Python is
case sensitive. That is, Python distinguishes between upper and lower case letters. Thus,
two words spelled the same but having different letters capitalized are treated as different
names in Python. Keep that in mind as we introduce different commands.
or this on a PC:
In [3]: pwd
Out[3]: C:\\Users\\pine
Typing cd .. (cd space two periods) moves the IPython shell up one directory
in the directory tree, as illustrated by the set of commands below.
In [4]: cd ..
/Users
In [5]: pwd
Out[5]: u/Users
The directory moved up one from /Users/pine to /Users. Now type ls (list) and
press RETURN. The console should list the names of the files and subdirectories in the
current directory.
In [6]: ls
Shared/
pine/
In this case, there are only two directories (indicated by the slash) and not files. Type cd
~ again to return to your home directory and then type pwd to verify where you are in
your directory tree. [Technically, ls isnt a magic command, but typing it without the %
sign lists the contents of the current directory, irrespective of whether Automagic is ON
or OFF.]
Lets create a directory within your documents directory that you can use to store your
Python programs. We will call it PyProgs. First, return to your home directory by typing
cd ~. To create directory PyProgs, type mkdir PyProgs (make directory). Type
ls to confirm that you have created PyProgs and then type cd PyProgs to switch to
that directory.
Now lets say you want to return to the previous subdirectory, Documents or My
Documents, which should be one up in the directory tree if you have followed along.
Type cd .. and then type pwd. You should find that you are back in the previous directory, Documents or My Documents. If you type ls, you should see the new directory
PyProgs that you just created.
More Magic Commands
The most important magic command is %run filename where filename is the name of a
Python program you have created. We havent done this yet but include it here just for
reference. We will come back to this later.
Some other useful magic commands include %hist, which lists the recent commands
issued to the IPython terminal, and %edit, which opens a new empty file in the code
editor window. Typing %edit filename, will open the file filename if it exists in the
current directory, or it will create a new file by that name if it does not, and will open it as
a blank file in the code editor window.
There are a number of other magic commands. You can get a list of them by typing
lsmagic.
In [7]: lsmagic
Available line magics:
%alias %alias_magic %autocall %automagic %bookmark %cd
%clear %colors %config %connect_info %debug %dhist %dirs
%doctest_mode %ed %edit %env %gui %guiref %hist %history
%install_default_config %install_ext %install_profiles
%killbgscripts %less %load %load_ext %loadpy %logoff %logon
%logstart %logstate %logstop %lsmagic %macro %magic %man
%more %notebook %page %pastebin %pdb %pdef %pdoc %pfile
%pinfo %pinfo2 %popd %pprint %precision %profile %prun
%psearch %psource %pushd %pwd %pycat %pylab %qtconsole
%quickref %recall %rehashx %reload_ext %rep %rerun %reset
%reset_selective %run %save %sc %store %sx %system %tb
%time %timeit %unalias %unload_ext %who %who_ls %whos
%xdel %xmode
Available cell magics:
%%! %%bash %%capture %%file %%javascript %%latex %%perl
%%prun %%pypy %%python %%python3 %%ruby %%script %%sh %%svg
%%sx %%system %%timeit
Automagic is ON, % prefix IS NOT needed for line magics.
There are a lot of magic commands, most of which we dont need right now. We will
introduce them in the text as needed.
10
Python returns the correct product, as expected. You can do more complicated calculations:
In [2]: 6+21/3
Out[2]: 13.0
Notice that the effect of the parentheses in In [3]: (6+21)/3 is to cause the addition to be performed first and then the division. Without the parentheses, Python will
always perform the multiplication and division operations before performing the addition
and subtraction operations. The order in which arithmetic operations are performed is the
same as for most calculators: exponentiation first, then multiplication or division, then
addition or subtraction, then left to right.
Symbol
+
*
/
//
%
**
Example
12+7
12-7
12*7
12/7
12//7
12%7
12**7
Output
19
5
84
1.714285
1
5
35831808
Floor division, designated by the symbols //, means divide and keep only the integer
part without rounding. Remainder, designated by the symbols %, gives the remainder of
after a floor division.
2.4. Interactive Python as a calculator
11
For the binary operators +, -, *, and //, the output is an integer if the inputs are integers.
The only exception is if the result of the calculation is out of the bounds of Python integers,
in which case Python automatically converts the result to a long integer. The output of
the division operator / is a floating point as of version 3 of Python. If an integer output is
desired when two integers are divided, the floor division operator // must be used. See
Important note on integer division in Python.
Floating point numbers are essentially rational numbers and can have a fractional part;
integers, by their very nature, have no fractional part. In most versions of Python running
on PCs or Macs, floating point numbers go between approximately 2 10308 and
2 10308 . Here are some examples of integer arithmetic:
In [9]: 12.*3.
Out[9]: 36.0
12
In [10]: 123.4*(-53.9)/sqrt(5.)
Out[10]: -2974.5338992050501
In [11]: 11./5.
Out[11]: 2.2
In [12]: 11.//5.
Out[12]: 2.0
In [13]: 11.%5.
Out[13]: 1.0
In [14]: 6.022e23*300.
Out[14]: 1.8066e+26
Note that the result of any operation involving only floating point numbers as inputs is a
real number, even in the cases where the floor division // or remainder % operators are
used. The last output also illustrates an alternative way of writing floating point numbers
as a mantissa followed by and e or E followed by a power of 10: so 1.23e-12 is equivalent
to 1.23 1012 .
We also sneaked into our calculations sqrt, the square root function. We will have more
to say about functions in a few pages.
Complex numbers are written in Python as a sum of a real and imaginary part. For
example, the complex number 32i is represented as 3-2j in Python where j represents
Notice that to obtain the expected result or 3, you must write the argument of the
square root function as a complex number. Otherwise, Python returns nan (not a number).
If you multiply an integer by a floating point number, the result is a floating point number.
Similarly, if you multiply a floating point number by a complex number, the result is a
13
complex number. Python always promotes the result to the most complex of the inputs.
14
15
The argument of the function can be a number or any kind of expression whose output
produces a number. For example, the function log(x) calculates the nature logarithm.
All of the following expressions are legal and produce the expected output:
In [2]: log(sin(0.5))
Out[2]: -0.73516668638531424
In [3]: log(sin(0.5)+1.0)
Out[3]: 0.39165386283471759
In [4]: log(5.5/1.2)
Out[4]: 1.5224265354444708
Description
Square root of x
Exponential of x, i.e. ex
Natural log of x, i.e. ln x
Base 10 log of x
Converts x from radians to degrees
Converts x from degrees to radians
Sine of x (x in radians)
Cosine x (x in radians)
Tangent x (x in radians)
Arc sine (in radians) of x
Arc cosine (in radians) of x
Arc tangent (in radians) of x
Absolute value of x
Rounds a float to nearest integer
Rounds a float down to nearest integer
Rounds a float up to nearest integer
-1 if x < 0, +1 if x > 0, 0 if x = 0
The functions discussed here all have one input and one output. Python functions can,
in general, have multiple inputs and multiple outputs. We will discuss these and other
features of functions later when we take up functions in the context of user-defined functions.
16
2.7 Variables
2.7.1 Names and the assignment operator
A variable is a name that is used to store data. It can be used to store different kinds of
data, but here we consider the simplest case where the data is a single numerical value.
Here are a few examples:
In [1]: a = 23
In [2]: p, q = 83.4, sqrt(2)
The equal sign = is the assignment operator. In the first statement, it assigns the value
of 23 to the variable a. In the second statement it assigns a value of 83.4 to p and a value
of 1.4142135623730951 to q. To be more precise, the name of a variable, such as a, is
associated with a memory location in your computer; the assignment variable tells the
computer to put a particular piece of data, in this case a numerical value, in that memory
location. Note that Python stores the numerical value, not the expression used to generate
it. Thus, q is assigned the 17-digit number
1.4142135623730951 generated by evaluating
the expression sqrt(2), not with 2. (Actually the value of q is stored as a binary,
base 2, number using scientific notation with a mantissa and an exponent.)
Suppose we write
In [3]: b = a
In this case Python associates a new memory location with the name b, distinct from the
one associated with a, and sets the value stored at that memory location to 23, the value
of a. The following sequence of statements demonstrate that fact. Can you see how?
Notice that simply typing a variable name and pressing Return prints out the value of
the variable.
In [4]: a=23
2.7. Variables
17
In [5]: b=a
In [6]: a
Out[6]: 23
In [7]: b
Out[7]: 23
In [8]: a=12
In [9]: a
Out[9]: 12
In [10]: b
Out[10]: 23
The assignment variable works from right to left; that is, it assigns the value of the number
on the right to the variable name on the left. Therefore, the statement 5=a makes no
sense in Python. The assignment operator = in Python is not equivalent to the equals
sign = we are accustomed to in algebra.
The assignment operator can be used to increment or change the value of a variable
In [11]: b = b+1
In [12]: b
Out[12]: 24
The statement, b = b+1 makes no sense in algebra, but in Python (and most computer
languages), it makes perfect sense: it means add 1 to the current value of b and assign
the result to b. This construction appears so often in computer programming that there
is a special set of operators to perform such changes to a variable: +=, -=, *=, and /=.
Here are some examples of how they work:
In [13]: c , d = 4, 7.92
In [14]: c += 2
In [15]: c
Out[15]: 6
In [16]: c *= 3
In [16]: c
Out[16]: 18
18
In [17]: d /= -2
In [17]: d
Out[17]: -3.96
In [18]: d -= 4
In [19]: d
Out[19]: -7.96
The variable names distance, time_traveled, and velocity immediately remind you of what is being calculated here. This is good practice. But so is keeping
variable names reasonably short, so dont go nuts!
2.7. Variables
19
and
as
assert
break
class
continue
def
del
elif
else
except
exec
finally
for
from
global
if
import
in
is
lambda
not
or
pass
print
raise
return
try
while
with
yield
In addition, you should not use function names, like sin, cos, and sqrt, defined in the
SciPy, NumPy, or any other library that you are using.
20
The number (or hash) symbol # is the comment character in Python; anything on a line
following # is ignored when the code is executed. Judicious use of comments in your
code will make your code much easier to understand days, weeks, or months after the
time you wrote it. Use comments generously.
Now you are ready to run the code. Before doing so, you first need to use the IPython
console to move to the PyProgs directory where the file containing the code resides.
From the IPython console, use the cd command to move to the PyProgs directory. For
example, you might type
In [1]: cd ~/Documents/PyProgs/
To run or execute a script, simply type run filename, which in this case means type run
myTrip.py. When you run a script, Python simply executes the sequence of commands
in the order they appear.
In [2]: run myTrip.py
Once you have run the script, you can see the values of the variables calculated in the
2.8. Script files and programs
21
script simply by typing the name of the variable. IPython responds with the value of that
variable.
In [3]: time
Out[3]: 6.666666666666667
In [4]: gallons
Out[4]: 13.333333333333334
In [5]: cost
Out[5]: 54.666666666666664
You can change the number of digits IPython displays using the command %precision:
In [6]: %precision 2
Out[6]: u%.2f
In [7]: time
Out[7]: 6.67
In [8]: gallons
Out[8]: 13.33
In [9]: cost
Out[9]: 54.67
the script will return the values of the variables time, gallons, and cost that the
script calculated. We will discuss the print function in much greater detail, as well as
other methods for data output, in Chapter 4 on Input and Output.
22
We have introduced extra spaces into some of the expressions to improve readability.
They are not necessary; where and whether you include them is largely a matter of taste.
There are two important differences between the code above and the commands we would
have written into the IPython console to execute the same set of commands. The first is
the statement on the second line
...
import numpy as np
...
and the second is the np. in front of the sqrt function on the last line. If you leave out
the import numpy as np line and remove the np. in front of the sqrt function,
you will get the following error message
----> 7 dr = sqrt( (x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2 )
NameError: name sqrt is not defined
The reason for the error is that the sqrt function is not a part of core Python. But it is a
part of the NumPy module discussed earlier. To make the NumPy library available to the
script, you need to add the statement import numpy as np. Then, when you call a
NumPy function, you need to write the function with the np. prefix. Failure to do either
will result in a error message. Now we can run the script.
In [10]: run twoPointDistance.py
In [11]: dr
Out[11]: 34.48
23
These statements import the entire library named in the import statement and associate
a prefix with the imported library: np and plt in the above examples. Functions from
within these libraries are then called by attaching the appropriate prefix with a period
before the function name. Thus, the functions sqrt or sin from the NumPy library are
called using the syntax np.sqrt or np.sin; the functions plot or xlabel from the
maplotlib.pyplot would be called using plt.plot or plt.xlabel.
Alternatively, the NumPy and MatPlotLib libraries can be called simply by writing
import numpy
import maplotlib.pyplot
When loaded this way, the sqrt function would be called as numpy.sqrt and the
plot function would be called as MatPlotLib.pyplot.plot. The import as
syntax allows you to define nicknames for numpy and maplotlib.pyplot. Nearly
any nickname can be chosen, but the Python community has settled on the nicknames np
and plt for numpy and maplotlib.pyplot, so you are advised to stick with those.
Using the standard nicknames makes your code more readable.
You can also import a single functions or subset of functions from a module without
importing the entire module. For example, suppose you wanted to import just the natural
24
which would assign the value 1.6094379124341003 to the variable a. If you wanted
to import the three functions, log, sin, and cos, you would write
from numpy import log, sin, cos
and would similarly use them without an np. prefix. In general, we do not recommend
using the the from module import ... way of importing functions. When reading
code, t makes it harder to determine from which modules functions are imported, and can
lead to clashes between similarly named functions from different modules. Nevertheless,
you will see the form used in programs you encounter on the web and elsewhere so it is
important to understand the syntax.
Often, the information provided can be quite extensive and you might find it useful to
clear the IPython window with the clear command so you can easily scroll back to find
the beginning of the documentation. You may have also noticed that when you type the
name of a function plus the opening parenthesis, IPython displays a window showing the
first dozen lines or so of the documentation on that function.
25
26
2.12 Exercises
1. A ball is thrown vertically up in the air from a height h0 above the ground at an
initial velocity v0 . Its subsequent height h and velocity v are given by the equations
h = h0 + v0 t 21 gt2
v = v0 gt
where g = 9.8 is the acceleration due to gravity in m/s2 . Write a script that finds
the height h and velocity v at a time t after the ball is thrown. Start the script by
setting h0 = 1.2 (meters) and v0 = 5.4 (m/s) and have your script print out the
values of height and velocity (see Note about printing). Then use the script to find
the height and velocity after 0.5 seconds. Then modify your script to find them after
2.0 seconds.
2. Write a script that defines the variables V0 = 10, a = 2.5, and z = 4 31 , and then
evaluates the expression
z
V = V0 1
.
a2 + z 2
Then find V for z = 8 32 and print it out (see Note about printing). Then find V for
z = 13 by changing the value of z in your script.
3. Write a single Python script that calculates the following expressions:
2 + e2.8
(a)
13 2
1 (1 + ln 2)3.5
1+ 5
!
2 2
(c) sin
2+ 2
(b)
After running your script in the IPython shell, typing a, b, or c at the IPython
prompt should yield the value of the expressions in (a), (b), or (c), respectively.
4. A quadratic equation with the general form
ax2 + bx + c = 0
has two solutions given by the quadratic formula
b b2 4ac
x=
.
2a
2.12. Exercises
27
(a) Given a, b, and c as inputs, write a script that gives the numerical values of
the two solutions. Write the constants a, b, and c as floats, and show that your
script gives the correct solutions for a few test cases when the solutions are real
numbers, that is, when the discriminant b2 4ac 0. Use the print function
in your script, discussed at the end of Section 2.8.1 Scripting Example 1, to
print out your two solutions.
(b) Written this way, however, your script gives an error message when the solutions are complex. For example, see what happens when a = 1, b = 2, and
c = 3. You can fix this using statements in your script like a = a+0j after
setting a to some float value. Thus, you can make the script work for any set
of real inputs for a, b, and c. Again, use the print function to print out your
two solutions.
28
CHAPTER
THREE
29
way than for lists. The elements of lists and arrays are numbered consecutively, and to
access an element of a list or an array, you simply refer to the number corresponding
to its position in the sequence. The elements of dictionaries are accessed by keys,
which can be either strings or (arbitrary) integers (in no particular order). Dictionaries
are an important part of core Python. However, we do not make much use of them in this
introduction to scientific Python, so our discussion of them is limited.
3.1 Strings
Strings are lists of characters. Any character that you can type from a computer keyboard,
plus a variety of other characters, can be elements in a string. Strings are created by
enclosing a sequence of characters within a pair of single or double quotes. Examples of
strings include "Marylyn", "omg", "good_bad_#5f>!", "{0:0.8g}", and "We
hold these truths ...". Caution: the defining quotes must both be single or
both be double quotes.
Strings can be assigned variable names
In [1]: a = "My dogs name is"
In [2]: b = "Bingo"
In forming the string c, we concatenated three strings, a, b, and a string literal, in this
case a space " ", which is needed to provide a space to separate string a from b.
You will use strings for different purposes: labeling data in data files, labeling axes in
plots, formatting numerical output, requesting input for your programs, as arguments in
functions, etc.
Because numbersdigitsare also alpha numeric characters, strings can be made up of
numbers:
In [5]: d = "927"
In [6]: e = 927
The variable d is a string while the variable e is an integer. What do you think happens if
you try to add them by writing d+e? Try it out and see if you understand the result.
30
3.2 Lists
Python has two data structures, lists and tuples, that consist of a list of one or more
elements. The elements of lists or tuples can be numbers or strings, or both. Lists (we will
discuss tuples later) are defined by a pair of square brackets on either end with individual
elements separated by commas. Here are two examples of lists:
In [1]: a = [0, 1, 1, 2, 3, 5, 8, 13]
In [2]: b = [5., "girl", 2+0j, "horse", 21]
We can access individual elements of a list using the variable name for the list with square
brackets:
In [3]: b[0]
Out[3]: 5.0
In [4]: b[1]
Out[4]: girl
In [5]: b[2]
Out[5]: (2+0j)
The first element of b is b[0], the second is b[1], the third is b[2], and so on. Some
computer languages index lists starting with 0, like Python and C, while others index
lists (or things more-or-less equivalent) starting with 1 (like Fortran and Matlab). Its
important to keep in mind that Python uses the former convention: lists are zero-indexed.
The last element of this array is b[4], because b has 5 elements. The last element can
also be accessed as b[-1], no matter how many elements b has, and the next-to-last
element of the list is b[-2], etc. Try it out:
In [6]: b[4]
Out[6]: 21
In [7]: b[-1]
Out[7]: 21
In [8]: b[-2]
Out[8]: horse
3.2. Lists
31
Here we see that 2 was added to the previous value of b[0] and the string horse
was replaced by the floating point number 3.14159. We can also manipulate individual
elements that are strings:
In [13]: b[1] = b[1] + "s & boys"
In [14]: b
Out[14]: [10.0, girls & boys, (2+0j), 3.14159, 21]
You can also add lists, but the result might surprise you:
In [15]: a
Out[15]: [0, 1, 1, 2, 3, 5, 8, 13]
In [16]: a+a
Out[16]: [0, 1, 1, 2, 3, 5, 8, 13, 0, 1, 1, 2, 3, 5, 8, 13]
In [17]: a+b
Out[17]: [0, 1, 1, 2, 3, 5, 8, 13, 10.0, girls & boys, (2+0j),
3.14159, 21]
You access a subset of a list by specifying two indices separated by a colon :. This
32
is a powerful feature of lists that we will use often. Here are a few other useful slicing
shortcuts:
In [21]: b[2:]
Out[21]: [(2+0j), 3.14159, 21]
In [22]: b[:3]
Out[22]: [10.0, girls & boys, (2+0j)]
In [23]: b[:]
Out[23]: [10.0, girls & boys, (2+0j), 3.14159, 21]
Thus, if the left slice index is 0, you can leave it out; similarly, if the right slice index is
the length of the list, you can leave it out also.
What does the following slice of an array give you?
In [24]: b[1:-1]
You can get the length of a list using Pythons len function:
In [25]: len(b)
Out[25]: 5
You can add one or more elements to the beginning or end of a list using the + operator:
In [29]: a = range(1,10,3)
3.2. Lists
33
In [30]: a
Out[30]: [1, 4, 7]
In [31]: a += [16, 31, 64, 127]
In [32]: a
Out[32]: [1, 4, 7, 16, 31, 64, 127]
In [33]: a = [0, 0] + a
In [34]: a
Out[34]: [0, 0, 1, 4, 7, 16, 31, 64, 127]
3.2.3 Tuples
Finally, a word about tuples: tuples are lists that are immutable. That is, once defined, the
individual elements of a tuple cannot be changed. Whereas a list is written as a sequence
of numbers enclosed in square brackets, a tuple is written as a sequence of numbers
enclosed in round parentheses. Individual elements of a tuple are addressed in the same
way as individual elements of lists are addressed, but those individual elements cannot be
changed. All of this illustrated by this simple example:
In [37]: c = (1, 1, 2, 3, 5, 8, 13)
In [37]: c[4]
Out[38]: 5
In [39]: c[4] = 7
--------------------------------------------------------------------------TypeError: tuple object does not support item assignment
When we tried to change c[4], the system returned an error because we are prohibited
from changing an element of a tuple. Tuples offer some degree of safety when we want
to define lists of immutable constants.
34
Here we have a three-element list where each element consists of a two-element list. Such
constructs can be useful in making tables and other structures. They also become relevant
later on in our discussion of NumPy arrays and matrices, which we introduce below.
We can access the various elements of a list with a straightforward extension of the indexing scheme we have been using. The first element of the list a above is a[0], which is
[3, 9]; the second is a[1], which is [8, 5]. The first element of a[1] is accessed
as a[1][0], which is 8, as illustrated below:
In [41]: a[0]
Out[41]: [3, 9]
In [42]: a[1]
Out[42]: [8, 5]
In [43]: a[1][0]
Out[43]: 8
In [44]: a[2][1]
Out[44]: 1
Multidimensional tuples work exactly like multidimensional lists, except they are immutable.
35
0,
0,
1,
4,
Notice that b is an integer array, as it was created from a list of integers. On the other hand,
c is a floating point array even though only one of the elements of the list from which it
was made was a floating point number. The array function automatically promotes all
of the numbers to the type of the most general entry in the list, which in this case is a
floating point number. In the case that elements of the list is made up of numbers and
strings, all the elements become strings when an array is formed from a list.
The second way arrays can be created is using the NumPy linspace or logspace
functions. The linspace function creates an array of N evenly spaced points between
a starting point and an ending point. The form of the function is linspace(start,
stop, N). If the third argument N is omitted, then N=50.
In [6]: linspace(0, 10, 5)
Out[6]: array([ 0. , 2.5,
5. ,
7.5, 10. ])
The linspace function produced 5 evenly spaced points between 0 and 10 inclusive.
NumPy also has a closely related function logspace that produces evenly spaced points
on a logarithmically spaced scale. The arguments are the same as those for linspace
except that start and stop refer to a power of 10. That is, the array starts at 10start
and ends at 10stop .
In [7]: %precision 1
Out[7]: u%.1f
In [8]: logspace(1, 3, 5)
Out[8]: array([
10. ,
31.6,
36
100. ,
316.2, 1000. ])
The logspace function created an array with 5 points evenly spaced on a logarithmic
axis starting at 101 and ending at 103 . The logspace function is particularly useful
when you want to create a log-log plot.
The third way arrays can be created is using the NumPy arange function, which is
similar to the Python range function for creating lists. The form of the function is
arange(start, stop, step). If the third argument is omitted step=1. If the
first and third arguments are omitted, then start=0 and step=1.
In [9]: arange(0, 10, 2)
Out[9]: array([0, 2, 4, 6, 8])
In [10]: arange(0., 10, 2)
Out[10]: array([ 0., 2., 4., 6., 8.])
In [11]: arange(0, 10, 1.5)
Out[11]: array([ 0. , 1.5, 3. , 4.5, 6. , 7.5, 9. ])
The arange function produces points evenly spaced between 0 and 10 exclusive of the
final point. Notice that arange produces an integer array in the first case but a floating
point array in the other two cases. In general arange produces an integer array if the
arguments are all integers; making any one of the arguments a float causes the array that
is created to be a float.
A fourth way to create an array is with the zeros and ones functions. As their names
imply, they create arrays where all the elements are either zeros or ones. They each take
on mandatory argument, the number of elements in the array, and one optional argument
that specifies the data type of the array. Left unspecified, the data type is a float. Here are
three examples
In [12]: zeros(6)
Out[12]: array([ 0., 0., 0., 0., 0., 0.])
In [13]ones(8)
Out[13]: array([ 1., 1., 1., 1., 1., 1., 1., 1.])
In [14]ones(8, dtype=int)
Out[14]: array([1, 1, 1, 1, 1, 1, 1, 1])
37
0.,
Here we can see that each element of the array has been multiplied by 6. This works
not only for multiplication, but for any other mathematical operation you can imagine:
division, exponentiation, etc.
In [18]: a/5
Out[18]: array([-0.2, 0. , 0.2, 0.4, 0.6, 0.8, 1. ])
In [19]: a**3
Out[19]: array([
-1.,
0.,
1.,
8.,
27.,
64., 125.])
In [20]: a+4
Out[20]: array([ 3., 4., 5., 6., 7., 8., 9.])
In [21]: a-10
38
4.,
6.,
In [23]: sin(a)
Out[23]: array([-0.84147098, 0.
, 0.84147098, 0.90929743,
0.14112001, -0.7568025 , -0.95892427])
In [24]: exp(-a)
Out[24]: array([ 2.71828183, 1.
, 0.36787944, 0.13533528,
0.04978707, 0.01831564, 0.00673795])
In [25]: 1. + exp(-a)
Out[25]: array([ 3.71828183, 2.
, 1.36787944, 1.13533528,
1.04978707, 1.01831564, 1.00673795])
In [26]: b = 5*ones(8)
In [27]: b
Out[27]: array([ 5., 5., 5., 5., 5., 5., 5., 5.])
In [28] b += 4
In [29] b
Out[29]: array([ 9., 9., 9., 9., 9., 9., 9., 9.])
In each case, you can see that the same mathematical operations are performed individually on each element of each array. Even fairly complex algebraic computations can be
carried out this way.
Lets say you want to create an x-y data set of y = cos x vs. x over the interval from -3.14
to 3.14. Here is how you might do it.
In [30]: x = linspace(-3.14, 3.14, 21)
In [31]: y = cos(x)
In [32]: x
Out[32]: array([-3.14 ,
-1.256,
0.628,
2.512,
In [33]: y
39
You can use arrays as inputs for any of the functions introduced in the section on Python
functions: a first look. You might well wonder what happens if Python encounters an
illegal operation. Here is one example.
In [34]: a
Out[34]: array([-1., 0., 1., 2., 3., 4., 5.])
In [35]: log(a)
-c:1: RuntimeWarning: divide by zero encountered in log
-c:1: RuntimeWarning: invalid value encountered in log
Out[35]: array([
nan, -inf, 0.
, 0.693, 1.099, 1.386,
1.609])
We see that NumPy calculates the logarithm where it can, and returns nan (not a number)
for an illegal operation, taking the logarithm of a negative number, and -inf, or for
the logarithm of zero. The other values in the array are correctly reported. NumPy also
prints out a warning message to let you know that something untoward has occurred.
Arrays can also be added, subtracted, multiplied, and divided by each other on an elementby-element basis, provided the two arrays have the same size. Consider adding the two
arrays a and b defined below:
In [36]: a = array([34., -12, 5.])
In [37]: b = array([68., 5.0, 20.])
In [38]: a+b
Out[38]: array([ 102.,
-7.,
25.])
The result is that each element of the two arrays are added. Similar results are obtained
for subtraction, multiplication, and division:
In [39]: a-b
Out[39]: array([-34., -17., -15.])
In [40]: a*b
Out[40]: array([ 2312.,
40
-60.,
100.])
In [41]: a/b
Out[41]: array([ 0.5 , -2.4 ,
0.25])
These kinds of operations with arrays are called vectorized operations because the entire
array, or vector, is processed as a unit. Vectorized operations are much faster than processing each element of arrays one by one. Writing code that takes advantage of these
kinds of vectorized operations is almost always to be preferred to other means of accomplishing the same task, both because it is faster and because it is usually syntactically
simpler. You will see examples of this later on when we discuss loops in Chapter 6.
We can get find the average velocity for time interval i by the formula
vi =
yi yi1
ti ti1
We can easily calculate the entire array of of velocities using the slicing and vectorized
subtraction properties of NumPy arrays by noting that we can create two y arrays displaced by one index
In [44]: y[:-1]
Out[44]: array([
0. ,
1.3,
5. ,
10.9,
18.9,
28.7])
In [45]: y[1:]
Out[45]: array([
1.3,
5. ,
10.9,
18.9,
28.7,
40. ])
41
In [46]: y[1:]-y[:-1]
Out[46]: array([ 1.3,
3.7,
5.9,
8. ,
9.8,
11.3])
Of course, these are the average velocities over each interval so the times best associated
with each interval are the times halfway in between the original time array, which we can
calculate using a similar trick of slicing:
In [49]: tv = (t[1:]+t[:-1])/2.
In [50]: tv
Out[50]: array([ 0.245,
0.745,
1.25 ,
1.79 ,
2.315,
2.875])
42
Notice the syntax used above in which two one-dimensional lists [1., 4, 5] and [9,
7, 4] are enclosed in square brackets to make a two-dimensional list. The array
function converts the two-dimensional list, a structure we introduced earlier, to a twodimensional array. When it makes the conversion from a list to an array, the array function
makes all the elements have the same data type as the most complex entry, in this case
a float. This points out an important difference between NumPy arrays and lists: all
elements of a NumPy array must be of the same data type: floats, or integers, or complex
numbers, etc.
There are a number of other functions for creating arrays. For example, a 3 row by 4
column array or 3 4 array with all the elements filled with 1 can be created using the
ones function introduced earlier.
In [53]: a = ones((3,4), dtype=float)
In [54]: a
Out[54]: array([[ 1., 1., 1., 1.],
[ 1., 1., 1., 1.],
[ 1., 1., 1., 1.]])
Using a tuple to specify the size of the array in the first argument of the ones function
creates a multidimensional array, in this case a two-dimensional array with the two elements of the tuple specifying the number of rows and columns, respectively. The zeros
function can be used in the same way to create a matrix or other multidimensional array
of zeros.
The eye(N) function creates an N N two-dimensional identity matrix with ones along
the diagonal:
In [55]: eye(4)
Out[55]: array([[
[
[
[
1.,
0.,
0.,
0.,
0.,
1.,
0.,
0.,
0.,
0.,
1.,
0.,
0.],
0.],
0.],
1.]])
Multidimensional arrays can also be created from one-dimensional arrays using the
reshape function. For example, a 2 3 array can be created as follows:
In [56]: c = arange(6)
In [57]: c
Out[57]: array([0, 1, 2, 3, 4, 5])
In [58]: c = reshape(c, (2,3))
43
In [59]: c
Out[59]: array([[0, 1, 2],
[3, 4, 5]])
which means the same thing. Caution: both the b[0][2] and the b[0,2] syntax work
for NumPy arrays and mean exactly the same thing; for lists only the b[0][2] syntax
works.
Matrix operations
Addition, subtraction, multiplication, division, and exponentiation all work with multidimensional arrays the same way they work with one dimensional arrays, on an elementby-element basis, as illustrated below:
In [62]: b
Out[62]: array([[ 1., 4., 5.],
[ 9., 7., 4.]])
In [63]: 2*b
Out[63]: array([[ 2., 8., 10.],
[ 18., 14., 8.]])
In [64]: b/4.
Out[64]: array([[ 0.25, 1. , 1.25],
[ 2.25, 1.75, 1. ]])
In [65]: b**2
Out[65]: array([[ 1., 16., 25.],
[ 81., 49., 16.]])
In [66]: b-2
44
Multiplying two arrays together is done on an element-by-element basis. Using the matrices b and c defined above, multiplying them together gives
In [68]: b
Out[68]: array([[ 1., 4., 5.],
[ 9., 7., 4.]])
In [69]: c
Out[69]: array([[0, 1, 2],
[3, 4, 5]])
In [70]: b*c
Out[70]: array([[ 0., 4., 10.],
[ 27., 28., 20.]])
Of course, this requires that both arrays have the same shape. Beware: array multiplication, done on an element-by-element basis, is not the same as matrix multiplication
as defined in linear algebra. Therefore, we distinguish between array multiplication and
matrix multiplication in Python.
Normal matrix multiplication is done with NumPys dot function. For example, defining
d to be the transpose of c using using the array functions T transpose method creates
an array with the correct dimensions that we can use to find the matrix product of b and
d:
In [71]: d = c.T
In [72]: d
Out[72]: array([[0, 3],
[1, 4],
[2, 5]])
In [73]: dot(b,d)
Out[73]: array([[ 14., 44.],
[ 15., 75.]])
45
3.4 Dictionaries
A Python list is a collection of Python objects indexed by an ordered sequence of integers
starting from zero. A dictionary is also collection of Python objects, just like a list,
but one that is indexed by strings or numbers (not necessarily integers and not in any
particular order) or even tuples! For example, suppose we want to make a dictionary of
room numbers indexed by the name of the person who occupies each room. We create
our dictionary using curly brackets {...}.
46
The dictionary above has three entries separated by commas, each entry consisting of a
key, which in this case is a string, and a value, which in this case is a room number. Each
key and its value are separated by a colon. The syntax for accessing the various entries is
similar to a that of a list, with the key replacing the index number. For example, to find
out the room number of Olivia, we type
In [2]: room["Olivia"]
Out[2]: 764
The key need not be a string; it can be any immutable Python object. So a key can be a
string, an integer, or even a tuple, but it cant be a list. And the elements accessed by their
keys need not be a string, but can be almost any legitimate Python object, just as for lists.
Here is a weird example
In [3]: weird = {"tank":52, 846:"horse", "bones":[23,
...: "fox", "grass"], "phrase":"I am here"}
In [4]: weird["tank"]
Out[4]: 52
In [5]: weird[846]
Out[5]: horse
In [6]: weird["bones"]
Out[6]: [23, fox, grass]
In [7]: weird["phrase"]
Out[7]: I am here
3.4. Dictionaries
47
You can get a list of all the keys or values of a dictionary by typing the dictionary name
followed by .keys() or .values().
In [13]: d.keys()
Out[13]: [last name, first name, birthday]
In [14]: d.values()
Out[14]: [Alberts, Marie, January 27]
In other languages, data types similar to Python dictionaries may be called hashmaps
or associative arrays, so you may see such term used if you read on the web about
dictionaries.
If rand has no argument, a single random number is generated. Otherwise, the argument
specifies the number of size of the array of random numbers that is created.
If you want random numbers uniformly distributed over some other interval, say from a
to b, then you can do that simply by stretching the interval so that it has a width of b a
48
and displacing the lower limit from 0 to a. The following statements produce random
numbers uniformly distributed from 10 to 20:
In [3]: a, b = 10, 20
In [4]: (b-a)*rand(20) + a
Out[4]: array([ 10.99031149,
18.25559651,
17.84258224,
12.05787676,
19.15872073,
10.16094915,
19.61360422,
18.11685555, 11.48302458,
17.55568817, 11.86290145,
12.1309852 , 14.30479884,
19.63135536, 16.58552886,
17.59104303, 11.48499468,
13.95534353, 18.21502143,
19.21058726])
The figure below shows histograms for the distributions of 10,000 random numbers generated by the np.random.rand (blue) and np.random.randn (green) functions.
As advertised, the np.random.rand function produces an array of random numbers
that is uniformly distributed between the values of 0 and 1, while np.random.randn
function produces an array of random numbers that follows a distribution of mean 0 and
standard deviation 1.
If we want a random numbers with a Gaussian distribution of width centered about x0 ,
we stretch the interval by a factor of and displace it by x0 . The following code produces
20 random numbers normally distributed around 15 with a width of 10:
In [5]: x0, sigma = 15, 10
In [6]: sigma*randn(20) + x0
Out[6]: array([ 9.36069244,
18.50471781,
12.45076637,
18.2756275 ,
20.8457924 ,
15.44237582,
31.19120157,
13.49260733, 6.12550102,
9.89499319, 14.09576728,
17.83073628, 2.95085564,
14.781659 , 31.80264078,
13.87890601, 25.41433678,
21.2385386 , -3.91668973,
26.24254326])
49
1.2
1.0
P(x)
0.8
0.6
0.4
0.2
0.0 4
50
the
see
51
3.6 Exercises
1. Create at array of 9 evenly spaced numbers going from 0 to 29 (inclusive) and give
it the variable name r. Find the square of each element of the array (as simply as
possible). Find twice the value of each element of the array in two different ways:
(i) using addition and (ii) using multiplication.
2. Create the following arrays:
(a) an array of 100 elements all equal to e, the base of the natural logarithm;
(b) an array in 1-degree increments of all the angles in degrees from 0 to 360
degrees;
(c) an array in 1-degree increments of all the angles in radians from 0 to 360
degrees;
(d) an array from 12 to 17, not including 17, in 0.2 increments;
(e) an array from 12 to 17, including 17, in 0.2 increments.
3. The position of a ball at time t dropped with zero initial velocity from a height h0
is given by
y = h0 12 gt2
2
where g = 9.8 m/s . Suppose h0 = 10 m. Find the sequence of times when the
ball passes each half meter assuming the ball is dropped at t = 0. Hint: Create a
NumPy array for y that goes from 10 to 0 in increments of -0.5 using the arange
function. Solving the above equation for t, show that
s
2(h0 y)
t=
.
g
Using this equation and the array you created, find the sequence of times when the
ball passes each half meter. Save your code as a Python script. It should yield the
following results for the y and t arrays:
In [2]: y
Out[2]: array([10. , 9.5, 9. , 8.5, 8. , 7.5, 7. , 6.5,
6. , 5.5, 5. , 4.5, 4. , 3.5, 3. , 2.5,
2. , 1.5, 1. , 0.5])
In [3]: t
Out[3]: array([ 0.
52
, 0.31943828, 0.45175395,
0.55328334,
0.7824608 ,
0.95831485,
1.10656667,
1.23717915,
1.35526185,
0.63887656, 0.71428571,
0.84515425, 0.9035079 ,
1.01015254, 1.05945693,
1.15175111, 1.19522861,
1.27775313, 1.31707778,
1.39239919])
3.6. Exercises
-4.9246827 ,
-7.3340579 ,
-9.12293148,
-10.61351563,
-11.91879801,
-13.09446421,
53
54
CHAPTER
FOUR
When the raw_input statement is executed, it prints to the computer screen the text in
the quotes and waits for input from the user. The user types a string of characters and
presses the return key. The raw_input function then assigns that string to the variable
name on the right of the assignment operator =.
55
Python prints out the string argument of the raw_input function and waits for a response from you. Lets go ahead and type 450 for 450 miles and press return. Now
type the variable name distance to see its value
In [2]: distance
Out[2]: u450
The value of the distance is 450 as expected, but it is a string (the u stands for unicode which refers to the string coding system Python uses). Because we want to use
450 as a number and not a distance, we need to convert it from a string to a number. We
can do that with the eval function by writing
In [3]: distance = eval(distance)
In [4]: distance
Out[4]: 450
The eval function has converted distance to an integer. This is fine and we are ready to
move on. However, we might prefer that distance be a float instead of an integer.
There are two ways to do this. We could assume the user is very smart and will type
450. instead of 450, which will cause distance to be a float when eval does the
conversion. That is, the number 450 is dynamically typed to be a float or an integer
depending on whether or not the user uses a decimal point. Alternatively, we could use
the function float in place of eval, which would ensure that distance is a floating
point variable. Thus, our code would look like this (including the user response):
In [5]: distance = raw_input("Input distance of trip in miles: ")
Input distance of trip in miles: 450
In [5]: distance
Out[5]: u450
In [7]: distance = float(distance)
In [8]: distance
Out[8]: 450.0
Now lets incorporate what we have learned into the code we wrote for Scripting Example
1
56
1
2
3
4
5
6
7
8
9
mpg = 30.
speed = 60.
costPerGallon = 4.10
# car mileage
# average speed
# price of gas
10
11
12
13
time = distance/speed
gallons = distance/mpg
cost = gallons*costPerGallon
Lines 4 and 5 can be combined into a single line, which is a little more efficient:
distance = float(raw_input("Input distance of trip in miles: "))
Whether you use float or int or eval depends on whether you want a float, an integer,
or a dynamically typed variable. In this program, it doesnt matter.
Now you can simply run the program and then type time, gallons, and cost to view
the results of the calculations done by the program.
Before moving on to output, we note that sometimes you may want string input rather
that numerical input. For example, you might want the user to input their name, in which
case you would simply use the raw_input function without converting its output.
57
The program prints out the results as a tuple of time (in hours), gasoline used (in gallons),
and cost (in dollars). Of course, the program doesnt give the user a clue as to which
quantity is which. The user has to know.
3
4
5
6
7
8
9
10
11
time = distance/speed
gallons = distance/mpg
cost = gallons*costPerGallon
12
13
14
15
16
17
Running this program, with the distance provided by the user, gives
In [9]: run myTripNiceIO.py
What is the trip distance in miles? 450
Duration of trip = 7.5 hours
Gasoline used = 15.0 gallons (@ 30 mpg)
Cost of gasoline = $61.50 (@ $4.10/gallon)
Now the output is presented in a way that is immediately understandable to the user.
Moreover, the numerical output is formatted with an appropriate number of digits to the
right of the decimal point. For good measure, we also included the assumed mileage
(30 mpg) and the cost of the gasoline. All of this is controlled by the str.format()
function within the print function.
The argument of the print function is of the form str.format() where str is a
58
string that contains text that is written to be the screen, as well as certain format specifiers
contained in curly braces {}. The format function contains the list of variables that are
to be printed.
The \n at the start of the string in the print statement on line 12 in the newline
character. It creates the blank line before the output is printed.
The positions of the curly braces determine where the variables in the format
function at the end of the statement are printed.
The format string inside the curly braces specifies how each variable in the format
function is printed.
The number before the colon in the format string specifies which variable in the list
in the format function is printed. Remember, Python is zero-indexed, so 0 means
the first variable is printed, 1 means the second variable, etc.
The zero after the colon specifies the minimum number of spaces reserved for printing out the variable in the format function. A zero means that only as many spaces
as needed will be used.
The number after the period specifies the number of digits to the right of the decimal
point that will be printed: 1 for time and gallons and 2 for cost.
The f specifies that a number with a fixed number of decimal points. If the f format
specifier is replaced with e, then the number is printed out in exponential format
(scientific notation).
In addition to f and e format types, there are two more that are commonly used: d for
integers (digits) and s for strings. There are, in fact, many more formatting possibilities. Python has a whole Format Specification Mini-Language that is documented at
http://docs.python.org/library/string.html#formatspec. Its very flexible but arcane. You
might find it simplest to look at the Format examples section further down the same
web page.
The program below illustrates most of the formatting you will need for writing a few
variables, be they strings, integers, or floats, to screen or to data files (which we discuss
in the next section).
string1 = "How"
string2 = "are you my friend?"
int1 = 34
int2 = 942885
float1 = -3.0
float2 = 3.141592653589793e-14
print( ***)
59
print(string1)
print(string1 + + string2)
print( 1. {} {}.format(string1, string2))
print( 2. {0:s} {1:s}.format(string1, string2))
print( 3. {0:s} {0:s} {1:s} - {0:s} {1:s}
.format(string1, string2))
print( 4. {0:10s}{1:5s}
.format(string1, string2))
print( ***)
print(int1, int2)
print( 6. {0:d} {1:d}.format(int1, int2))
print( 7. {0:8d} {1:10d}.format(int1, int2))
print( ***)
print( 8. {0:0.3f}.format(float1))
print( 9. {0:6.3f}.format(float1))
print(10. {0:8.3f}.format(float1))
print(2*11. {0:8.3f}.format(float1))
print( ***)
print(12. {0:0.3e}.format(float2))
print(13. {0:10.3e}.format(float2))
print(14. {0:10.3f}.format(float2))
print( ***)
print(15. 12345678901234567890)
print(16. {0:s}--{1:8d},{2:10.3e}
.format(string2, int1, float2))
60
12. 3.142e-14
13. 3.142e-14
14.
0.000
***
15. 12345678901234567890
16. are you my friend?--
34, 3.142e-14
Successive empty brackets {} like those that appear in the statement above print(
1. {} {}.format(string1, string2)) are numbered consecutively starting at 0 and will print out whatever variables appear inside the format() method using
their default format.
Finally, note that the code starting on lines 14 and 16 each are split into two lines. We
have done this so that the lines fit on the page without running off the edge. Python allows
you to break lines up like this to improve readability.
Simply using the print function does print out the array, but perhaps not in the
format you desire. To control the output format, you use the NumPy function
set_printoptions. For example, suppose you want to see no more than two digits
to the right of the decimal point. Then you simply write
In [12]: set_printoptions(precision=2)
In [13]: print(a)
[ 3.
5.67
8.33 11.
13.67 16.33
19.
If you want to change the number of digits to the right of the decimal point to 4, you set
the keyword argument precision to 4
In [14]: set_printoptions(precision=4)
In [15]: print(a)
[ 3.
5.6667
8.3333 11.
13.6667
16.3333
19.
Suppose you want to use scientific notation. The method for doing it is somewhat arcane,
using something called a lambda function. For now, you dont need to understand how
4.2. Screen output
61
it works to use it. Just follow the examples shown below, which illustrate several different
output formats using the print function with NumPy arrays.
In [16]: set_printoptions(
...: formatter={float: lambda x: format(x, 6.2e)})
In [17]: print(a)
[3.00e+00 5.67e+00 8.33e+00 1.10e+01 1.37e+01 1.63e+01 1.90e+01]
To specify the format of the output, you use the formatter keyword argument. The
first entry to the right of the curly bracket is a string that can be float, as it is above,
or int, or str, or a number of other data types that you can look up in the online
NumPy documentation. The only other thing you should change is the format specifier
string. In the above example, it is 6.2e, specifying that Python should allocate at least
6 spaces, with 2 digits to the right of the decimal point in scientific (exponential) notation.
For fixed width floats with 3 digits to the right of the decimal point, use the f in place of
the e format specifier, as follows
In [18]:
...:
In [19]:
[ 3.000
set_printoptions(
formatter={float: lambda x: format(x, 6.3f)})
print(a)
5.667 8.333 11.000 13.667 16.333 19.000]
62
data point
0
1
2
3
4
5
6
7
8
9
10
11
12
time (sec)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
height (mm)
180
182
178
165
160
148
136
120
99
83
55
35
5
uncertainty (mm)
3.5
4.5
4.0
5.5
2.5
3.0
2.5
3.0
4.0
2.5
3.6
1.75
0.75
We would like to read these data into a Python program, associating the data in each
column with an appropriately named array. While there are a multitude of ways to do this
in Python, the simplest by far is to use the NumPy loadtxt function, whose use we
illustrate here. Suppose that the name of the text file is MyData.txt. Then we can read
the data into four different arrays with the following statement
In [1]: dataPt, time, height, error = np.loadtxt("MyData.txt",
skiprows=5 , unpack=True)
In this case, the loadtxt function takes three arguments: the first is a string that is
the name of the file to be read, the second tells loadtxt to skip the first 5 lines at the
top of file, sometimes called the header, and the third tells loadtxt to output the data
(unpack the data) so that it can be directly read into arrays. loadtxt reads however
many columns of data are present in the text file to the array names listed to the left of the
= sign. The names labeling the columns in the text file are not used, but you are free to
choose the same or similar names, of course, as long as they are legal array names. By the
way, for the above loadtxt call to work, the file MyData.txt should be in the current
working directory of the IPython shell. Otherwise, you need to specify the directory path
with the file name.
It is critically important that the data file be a text file. It cannot be a MSWord file,
for example, or an Excel file, or anything other than a plain text file. Such files can be
created by a text editor programs like Notepad and Notepad++ (for a PC) or TextEdit and
TextWrangler (for a Mac). They can also be created by MSWord and Excel provided you
explicitly save the files as text files. Beware: You should exit any text file you make and
save it with a program that allows you to save the text file using UNIX-type formatting,
which uses a line feed (LF) to end a line. Some programs, like MSWord under Windows,
may include a carriage return (CR) character, which can confuse loadtxt. Note that
63
we give the file name a .txt extension, which indicates to most operating systems that
this is a text file, as opposed to an Excel file, for example, which might have a .xlsx or
.xls extension.
If you dont want to read in all the columns of data, you can specify which columns to
read in using the usecols key word. For example, the call
In [2]: time, height = loadtxt(MyData.txt, skiprows=5 ,
usecols = (1,2), unpack=True)
reads in only columns 1 and 2; columns 0 and 3 are skipped. As a consequence, only
two array names are included to the left of the = sign, corresponding to the two column
that are read. Writing usecols = (0,2,3) would skip column 1 and read in only the
data in colums 0, 2, and 3. In this case, 3 array names would need to be provided on the
left hand side of the = sign.
One convenient feature of the loadtxt function is that it recognizes any white space as
a column separator: spaces, tabs, etc.
Finally you should remember that loadtxt is a NumPy function. So if you are using it in
a Python module, you must be sure to include an import numpy as np statement
before calling np.loadtxt.
64
65
5,2.5,148,3
6,3,136,2.5
7,3.5,120,3
8,4,99,4
9,4.5,83,2.5
10,5,55,3.6
11,5.5,35,1.75
12,6,5,0.75
As its name suggests, the CSV file is simply a text file with the data that was formerly in
spreadsheet columns now separated by commas. We can read the data in this file into a
Python program using the loadtxt NumPy function once again. Here is the code
In [3]: dataPt, time, height, error = loadtxt("MyData.csv",
skiprows=5 , unpack=True, delimiter=,)
The form of the function is exactly the same as before except we have added the argument
delimiter=, that tells loadtxt that the columns are separated by commas instead
of white space (spaces or tabs), which is the default. Once again, we set the skiprows
argument to skip the header at the beginning of the file and to start reading at the first row
of data. The data are output to the arrays to the right of the assignment operator = exactly
as in the previous example.
We illustrate savetext below with a script that first creates four arrays by reading in the
data file MyData.txt, as discussed in the previous section, and then writes that same
data set to another file MyDataOut.txt.
1
import numpy as np
2
3
66
skiprows=5, unpack=True)
4
5
6
7
np.savetxt(MyDataOut.txt,
zip(dataPt, time, height, error), fmt="%12.1f")
The first argument of of savetxt is a string, the name of the data file to be created. Here
we have chosen the name MyDataOut.txt, inserted with quotes, which designates it
as a string literal. Beware, if there is already a file of that name on your computer, it will
be overwrittenthe old file will be destroyed and a new one will be created.
The second argument is the data array the is to be written to the data file. Because we
want to write not one but four data arrays to the file, we have to package the four data
arrays as one, which we do using the zip function, a Python function that combines
returns a list of tuples, where the ith tuple contains the ith element from each of the
arrays (or lists, or tuples) listed as its arguments. Since there are four arrays, each row
will be a tuple with four entries, producing a table with four columns. Note that the first
two arguments, the filename and data array, are regular arguments and thus must appear
as the first and second arguments in the correct order. The remaining arguments are all
keyword arguments, meaning that they are optional and can appear in any order, provided
you use the keyword.
The next argument is a format string that determines how the elements of the array are
displayed in the data file. The argument is optional and, if left out, is the format 0.18e,
which displays numbers as 18 digit floats in exponential (scientific) notation. Here we
choose a different format, 12.1f, which is a float displayed with 1 digit to the right of
the decimal point and a minimum width of 12. By choosing 12, which is more digits than
any of the numbers in the various arrays have, we ensure that all the columns will have
the same width. It also ensures that the decimal points in column of numbers are aligned.
This is evident in the data file below, MyDataOut.txt, which was produced by the above
script.
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
180.0
182.0
178.0
165.0
160.0
148.0
136.0
120.0
99.0
83.0
55.0
35.0
5.0
3.5
4.5
4.0
5.5
2.5
3.0
2.5
3.0
4.0
2.5
3.6
1.8
0.8
67
We omitted the optional delimiter keyword argument, which leaves the delimiter as
the default space.
We also omitted the optional header keyword argument, which is a string variable that
allows you to write header text above the data. For example, you might want to label the
data columns and also include the information that was in the header of the original data
file. To do so, you just need to create a string with the information you want to include
and then use the header keyword argument. The code below illustrates how to do this.
1
import numpy as np
2
3
4
5
6
7
8
9
10
info
info
info
info
info
11
12
13
np.savetxt(MyDataOut.txt,
zip(dataPt, time, height, error), header=info, fmt="%12.1f")
Now the data file produces has a header preceding the data. Notice that the header rows
all start with a # comment character, which is the default setting for the savetxt function. This can be changed using the keyword argument comments. You can find more
information about savetxt using the IPython help function or from the online NumPy
documentation.
# Data for falling mass experiment
# Date: 16-Aug-2013
# Data taken by Lauren and John
#
#
data point
time (sec) height (mm)
0.0
0.0
180.0
1.0
0.5
182.0
2.0
1.0
178.0
3.0
1.5
165.0
4.0
2.0
160.0
5.0
2.5
148.0
6.0
3.0
136.0
7.0
3.5
120.0
8.0
4.0
99.0
9.0
4.5
83.0
10.0
5.0
55.0
68
uncertainty (mm)
3.5
4.5
4.0
5.5
2.5
3.0
2.5
3.0
4.0
2.5
3.6
11.0
12.0
5.5
6.0
35.0
5.0
1.8
0.8
This data file, with a csv extension, can be directly read into a spreadsheet program like
Excel.
69
4.5 Exercises
1. Write a Python program that calculates how much money you can spend each day
for lunch for the rest of the month based on todays date and how much money you
currently have in your lunch account. The program should ask you: (1) how much
money you have in your account, (2) what todays date is, and (3) how many days
there are in month. The program should return your daily allowance. The results of
running your program should look like this:
How much money (in dollars) in your lunch account? 118.39
What day of the month is today? 17
How many days in this month? 30
You can spend $8.46 each day for the rest of the month.
Extra: Create a dictionary (see Dictionaries) that stores the number of days in each
month (forget about leap years) and have your program ask what month it is rather
than the number of days in the month.
2. From the IPython terminal, create the following three NumPy arrays:
a = array([1, 3, 5, 7])
b = array([8, 7, 5, 4])
c = array([0, 9,-6,-8])
Print d out at the terminal prompt. What kind of object is d? Hint: It is not a NumPy
array. Convert d into a NumPy array and call that array e. Type e at the terminal
prompt so that e is printed out on the IPython terminal. One of the elements of e is
-8. Show how to address and print out just that element of e. Show how to address
that same element of d. What has become of the three original arrays a, b, and c,
that is, how do they appear in e?
3. Create the following data file and then write a Python script to read it into a three
NumPy arrays with the variable names f, a, da for the frequency, amplitude, and
amplitude error.
Date: 2013-09-16
Data taken by Liam and Selena
frequency (Hz) amplitude (mm)
70
0.7500
1.7885
2.8269
3.8654
4.9038
5.9423
6.9808
8.0192
9.0577
10.0962
11.1346
12.1731
13.2115
14.2500
13.52
12.11
14.27
16.60
22.91
35.28
60.99
33.38
17.78
10.99
7.47
6.72
4.40
4.07
0.32
0.92
0.73
2.06
1.75
0.91
0.99
0.36
2.32
0.21
0.48
0.51
0.58
0.63
Show that you have correctly read in the data by having your script print out to your
computer screen the three arrays. Format the printing so that it produces output like
this:
f =
[ 0.75
1.7885
2.8269
3.8654
4.9038
5.9423
6.9808
8.0192
9.0577 10.0962 11.1346 12.1731
13.2115 14.25 ]
a =
[ 13.52 12.11 14.27 16.6
22.91 35.28 60.99 33.38
17.78 10.99
7.47
6.72
4.4
4.07]
da =
[ 0.32 0.92 0.73 2.06 1.75 0.91 0.99 0.36 2.32
0.21 0.48 0.51 0.58 0.63]
Note that the array f is displayed with four digits to the right of the decimal point
while the arrays a and da are displayed with only two. The columns of the displayed arrays need not line up as they do above.
4. Write a script to read the data from the previous problem into three NumPy arrays
with the variable names f, a, da for the frequency, amplitude, and amplitude error
and then, in the same script, write the data out to a data file, including the header,
with the data displayed in three columns, just as its displayed in the problem above.
Its ok if the header lines begin with the # comment character. Your data file should
have the extension .txt.
5. Write a script to read the data from the previous problem into three NumPy arrays
with the variable names f, a, da for the frequency, amplitude, and amplitude error
and then, in the same script, write the data out to a csv data file, without the header,
to a data file with the data displayed in three columns. Use a single format specifier
4.5. Exercises
71
and set it to "%0.16e". If you have access the spreadsheet program (like MS
Excel), try opening the file you have created with your Python script and verify that
the arrays are displayed in three columns. Note that your csv file should have the
extension .csv.
72
CHAPTER
FIVE
PLOTTING
The graphical representation of dataplottingis one of the most important tools for
evaluating and understanding scientific data and theoretical predictions. However, plotting is not a part of core Python but is provided through one of several possible library
modules. The most highly developed and widely used plotting package for Python is
MatPlotLib (http://MatPlotLib.sourceforge.net/). It is a powerful and flexible program
that has become the de facto standard for 2-d plotting with Python.
Because MatPlotLib is an external libraryin fact its a collection of librariesit must
be imported into any routine that uses it. MatPlotLib makes extensive use of NumPy so
the two should be imported together. Therefore, for any program for which you would
like to produce 2-d plots, you should include the lines
import numpy as np
import matplotlib.pyplot as plt
There are other MatPlotLib sub-libraries, but the pyplot library provides nearly everything that you need for 2-d plotting. The standard prefix for it is plt. MatPlotLib is automatically loaded with the IPython shell so you do not need to use import matplotlib.pyplot
nor do you need to use the plt prefix when working in the IPython shell.
One final word before we get started: We only scratch the surface of what is possible using
MatPlotLib and as you become familiar with it, you will surely want to do more than this
manual describes. In that case, you need to go the the web to get more information. A
good place to start is http://matplotlib.org/api/pyplot_summary.html. Another interesting
web page is http://matplotlib.org/gallery.html.
73
5.0
4.5
4.0
3.5
3.0
2.5
2.0
1.5
1.00
74
Chapter 5. Plotting
plot(x, y)
where x and y are arrays (or lists) that have the same size. If the x array is missing, that
is, if there is only a single array, as in our example above, the plot function uses 0, 1,
..., N-1 for the x array, where N is the size of the y array. Thus, the plot function
provides a quick graphical way of examining a data set.
More typically, you supply both an x and a y data set to plot. Taking things a bit further,
you may also want to plot several data sets on the same graph, use symbols as well as
lines, label the axes, create a title and a legend, and control the color of symbols and lines.
All of this is possible but requires calling a number of plotting functions. For this reason,
plotting is usually done using a Python script or program.
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 4.*np.pi, 33)
y = np.sin(x)
plt.plot(x, y)
plt.show()
75
Chapter 5. Plotting
1
2
import numpy as np
import matplotlib.pyplot as plt
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# create plot
plt.figure(1, figsize = (6,4) )
plt.plot(x, y, b-, label=theory)
plt.plot(xdata, ydata, ro, label="data")
plt.xlabel(x)
plt.ylabel(transverse displacement)
plt.legend(loc=upper right)
plt.axhline(color = gray, zorder=-1)
plt.axvline(color = gray, zorder=-1)
20
21
22
23
24
77
plt.show()
1.0
transverse displacement
25
theory
data
0.5
0.0
0.5
1.010
0
x
10
78
Chapter 5. Plotting
79
plot function draws between data points will be visible. For plotting a typical function,
something on the order of 100-200 data points usually produces a smooth curve, depending on just how curvy the function is. On the other hand, only two points are required to
draw a smooth straight line.
Detailed information about the MatPlotLib plotting functions are available online, starting
with the site http://matplotlib.org/api/pyplot_summary.html. The main MatPlotLib site is
http://matplotlib.org/.
description
solid line style
dashed line style
dash-dot line style
dotted line style
point marker
pixel marker
circle marker
triangle_down marker
triangle_up marker
triangle_left marker
triangle_right marker
tri_down marker
tri_up marker
character
3
4
s
p
*
h
H
+
x
D
d
|
_
description
tri_left marker
tri_right marker
square marker
pentagon marker
star marker
hexagon1 marker
hexagon2 marker
plus marker
x marker
diamond marker
thin_diamond marker
vline marker
hline marker
This second table gives the character codes for eight different colors. Many more are
possible but the color specification becomes more complex. You can consult the webbased MatPlotLib documentation for further details.
80
Chapter 5. Plotting
character
b
g
r
c
m
y
k
w
color
blue
green
red
cyan
magenta
yellow
black
white
Here are some examples of how these format specifiers can be used:
plot(x, y, ro)
plot(x, y, ks-)
plot(x, y, g^)
plot(x, y, k-)
plot(x, y, ms)
You can also make two calls sequentially for added versatility. For example, by sequentially calling the last two plot calls, the plot produces magenta squares on top of black
lines connecting the data points.
These format specifiers give rudimentary control of the plotting symbols and lines. MatPlotLib provides much more precise and detailed control of the plotting symbol size, line
types, and colors using optional keyword arguments instead of the plotting format strings
introduced above. For example, the following command creates a plot of large yellow
diamond symbols with blue edges connected by a green dashed line:
plot(x, y, color=green, linestyle=dashed, marker=d,
markerfacecolor=yellow, markersize=12,
markeredgecolor=blue)
Try it out! The online MatPlotLib documentation provides all the plotting format keyword
arguments and their possible values.
81
transverse displacement
20
theory
data
15
10
5
0
50
10
15
20
25
30
35
40
45
import numpy as np
import matplotlib.pyplot as plt
3
4
5
82
Chapter 5. Plotting
6
7
8
9
10
11
12
13
14
15
16
17
18
# create plot
plt.figure(1, figsize = (6,4) )
plt.plot(x, y, b-, label="theory")
plt.errorbar(xdata, ydata, fmt=ro, label="data",
xerr=0.75, yerr=yerror, ecolor=black)
plt.xlabel(x)
plt.ylabel(transverse displacement)
plt.legend(loc=upper right)
19
20
21
22
23
24
We have more to say about the errorbar function in the sections on logarithmic plots.
But the brief introduction given here should suffice for making most plots not involving
logarithmic axes.
83
plt.figure()
plt.plot(theta, ytan)
plt.show()
1400
1200
1000
800
600
400
200
0
2000
10
The resulting plot, shown above, doesnt quite look like what you might have expected
for tan vs . The problem is that tan diverges at = /2, 3/2, 5/2, ..., which
leads to large spikes in the plots as values in the theta array come near those values.
Of course, we dont want the plot to extend all the way out to in the y direction, nor
can it. Instead, we would like the plot to extend far enough that we get the idea of what is
going on as y , but we would still like to see the behavior of the graph near y = 0.
We can restrict the range of ytan values that are plotted using the MatPlotLib function
ylim, as we demonstrate in the script below.
import numpy as np
import matplotlib.pyplot as plt
theta = np.arange(0.01, 10., 0.04)
ytan = np.tan(theta)
plt.figure()
plt.plot(theta, ytan)
plt.ylim(-8, 8)
# restricts range of y axis from -8 to +8
plt.axhline(color="gray", zorder=-1)
84
Chapter 5. Plotting
plt.show()
The figure produced by this script is shown below. The plot now looks much more like
the familiar tan function we know. We have also include a call to the axline function
to create an x axis.
8
6
4
2
0
2
4
6
80
10
85
8
6
4
2
0
2
4
6
80
10
import numpy as np
import matplotlib.pyplot as plt
3
4
5
6
7
8
9
10
11
plt.figure()
plt.plot(theta, ytanM)
plt.ylim(-8, 8)
plt.axhline(color="gray", zorder=-1)
12
13
plt.show()
86
Chapter 5. Plotting
5.2.4 Subplots
cot(theta)
tan(theta)
Often you want to create two or more graphs and place them next to one another, generally
because they are related to each other in some
p way. The plot below shows an example of
2 1 vs are plotted. The two curves
such a plot. In the top graph, tan and (8/)p
(8/)2 1. In the bottom cot and
cross
p each other at the points where tan =
(8/)2p 1 vs are plotted. These two curves cross each other at the points where
cot = (8/)2 1.
8
6
4
2
0
2
4
6
80
8
6
4
2
0
2
4
6
80
4
theta
4
theta
import numpy as np
import matplotlib.pyplot as plt
3
4
5
6
7
8
87
10
11
plt.figure(1)
12
13
14
15
16
17
18
19
20
21
22
plt.subplot(2, 1, 1)
plt.plot(theta, y)
plt.plot(theta, ytan)
plt.ylim(-8, 8)
plt.axhline(color="gray", zorder=-1)
plt.axvline(x=np.pi/2., color="gray", linestyle=--, zorder=-1)
plt.axvline(x=3.*np.pi/2., color="gray", linestyle=--, zorder=-1)
plt.axvline(x=5.*np.pi/2., color="gray", linestyle=--, zorder=-1)
plt.xlabel("theta")
plt.ylabel("tan(theta)")
23
24
25
26
27
28
29
30
31
32
plt.subplot(2, 1, 2)
plt.plot(theta, -y)
plt.plot(theta, ycot)
plt.ylim(-8, 8)
plt.axhline(color="gray", zorder=-1)
plt.axvline(x=np.pi, color="gray", linestyle=--, zorder=-1)
plt.axvline(x=2.*np.pi, color="gray", linestyle=--, zorder=-1)
plt.xlabel("theta")
plt.ylabel("cot(theta)")
33
34
plt.show()
The function subplot, called on lines 13 and 24, creates the two subplots in the above
figure. subplot has three arguments. The first specifies the number of rows that the
figure space is to be divided into; on line 13, its two. The second specifies the number
of columns that the figure space is to be divided into; on line 13, its one. The third
argument specifies which rectangle the will contain the plot specified by the following
function calls. Line 13 specifies that the plotting commands that follow will be act on
the first box. Line 24 specifies that the plotting commands that follow will be act on the
second box.
We have also labeled the axes and included dashed vertical lines at the values of where
tan and cot diverge.
88
Chapter 5. Plotting
9000
8000
7000
6000
5000
4000
3000
2000
1000
00
104
theory
data
counts per second
For data sets that vary exponentially in the independent variable, it is often useful to use
one or more logarithmic axes. Radioactive decay of unstable nuclei, for example, exhibits
an exponential decrease in the number of particles emitted from the nuclei as a function
of time. In the plot below, for example, we show the decay of the radioactive isotope
Phosphorus-32 over a period of 6 months, where the radioactivity is measured once each
week. Starting at a decay rate of nearly 104 electrons (counts) per second, the decay rate
diminishes to only about 1 count per second after about 6 months or 180 days. If we plot
counts per second as a function of time on a normal plot, as we have done in the plot on
the left below, then the count rate is indistinguishable from zero after about 100 days. On
the other hand, if we use a logarithmic axis for the count rate, as we have done in the plot
on the right below, then we can follow the count rate well past 100 days and can readily
distinguish it from zero. Moreover, if the data vary exponentially in time, then the data
will fall along a straight line, as they do for the case of radioactive decay.
theory
data
103
102
101
100 0
89
1
2
import numpy as np
import matplotlib.pyplot as plt
3
4
5
6
7
8
9
10
11
12
13
14
# create plot
plt.figure(1, figsize = (10,4) )
15
16
17
18
19
20
21
plt.subplot(1, 2, 1)
plt.plot(t, N, b-, label="theory")
plt.plot(time, counts, ro, label="data")
plt.xlabel(time (days))
plt.ylabel(counts per second)
plt.legend(loc=upper right)
22
23
24
25
26
27
28
plt.subplot(1, 2, 2)
plt.semilogy(t, N, b-, label="theory")
plt.semilogy(time, counts, ro, label="data")
plt.xlabel(time (days))
plt.ylabel(counts per second)
plt.legend(loc=upper right)
29
30
plt.tight_layout()
31
32
33
The semilogx and semilogy functions work the same way as the plot function.
You just use one or the other depending on which axis you want to be logarithmic.
The tight_layout() function
You may have noticed the tight_layout() function, called without arguments on
line 30 of the program. This is a convenience function that adjusts the sizes of the plots to
make room for the axes labels. If it is not called, the y-axis label of the right plot runs into
the left plot. The tight_layout() function can also be useful in graphics windows
90
Chapter 5. Plotting
5
6
7
import numpy as np
import matplotlib.pyplot as plt
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
91
exponential
20000
104
distance (mm)
distance (mm)
15000
10000
5000
00
exponential
105
103
102
101
4
6
time (ms)
1.0
10
100 0
sinc function
4
6
time (ms)
10
0.8
electric field
0.6
0.4
0.2
0.0
0.2
0.410
0
angle (deg)
10
92
Chapter 5. Plotting
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
After defining several arrays for plotting, the above program opens a figure window in
line 23 with the statement
fig = plt.figure(figsize=(9,8))
The MatPlotLib statement above creates a Figure object, assigns it the name fig, and
opens a blank figure window. Thus, just as we give lists, arrays, and numbers variable
names (e.g. a = [1, 2, 5, 7], dd = np.array([2.3, 5.1, 3.9]), or st
= 4.3), we can give a figure object and the window in creates a name: here it is fig.
In fact we can use the figure function to open up multiple figure objects with different
figure windows. The statements
fig1 = plt.figure()
fig2 = plt.figure()
93
open up two separate windows, one named fig1 and the other fig2. We can then use the
names fig1 and fig2 to plot things in either window. The figure function need not
take any arguments if you are satisfied with the default settings such as the figure size and
the background color. On the other hane, by supplying one or more keyword arguments,
you can customize the figure size, the background color, and a few other properties. For
example, in the program listing (line 23), the keyword argument figsize sets the width
and height of the figure window; the default size is (8, 6); in our program we set it to
(9, 8), which is a bit wider and higher than the default size. In the example above, we
also choose to open only a single window, hence the single figure call.
The fig.add_subplot(2,2,1) in line 30 is a MatPlotLib function that divides the
figure window into 2 rows (the first argument) and 2 columns (the second argument). The
third argument creates a subplot in the first of the 4 subregions (i.e. of the 2 rows 2
columns) created by the fig.add_subplot(2,2,1) call. To see how this works,
type the following code into a Python module and run it:
1
2
import numpy as np
import matplotlib.pyplot as plt
3
4
5
fig = plt.figure(figsize=(9,8))
ax1 = fig.add_subplot(2,2,1)
6
7
plt.show()
You should get a figure window with axes drawn in the upper left quadrant. The
fig. prefix used with the add_subplot(2,2,1) function directs Python to draw
these axes in the figure window named fig. If we had opened two figure windows,
changing the prefix to correspond to the name of one or the other of the figure windows would direct the axes to be drawn in the appropriate window. Writing ax1 =
fig.add_subplot(2,2,1) assigns the name ax1 to the axes in the upper left quadrant of the figure window.
The ax1.plot(x, y) in line 27 directs Python to plot the previously-defined x
and y arrays onto the axes named ax1. The ax2 = fig.add_subplot(2,2,2)
draws axes in the second, or upper right, quadrant of the figure window. The ax3 =
fig.add_subplot(2,1,2) divides the figure window into 2 rows (first argument)
and 1 column (second argument), creates axes in the second or these two sections, and
assigns those axes (i.e. that subplot) the name ax3. That is, it divides the figure window into 2 halves, top and bottom, and then draws axes in the half number 2 (the third
argument), or lower half of the figure window.
You may have noticed in above code that some of the function calls are a bit different
from those used before: xlabel(time (ms)) becomes set_xlabel(time
(ms)), title(exponential) becomes set_title(exponential),
94
Chapter 5. Plotting
etc.
The call ax2.set_yscale(log) sets the y-axes in the second plot to be logarithmic, thus creating a semi-log plot. Creating properly-labeled logarthmic axes like this is
more straightforward with the advanced syntax illustrated in the above example.
Using the prefixes ax1, ax2, or ax3, direct graphical instructions to their respective
subplots. By creating and specifying names for the different figure windows and subplots
within them, you access the different plot windows more efficiently. For example, the
following code makes four identical subplots in a single figure window using a for loop.
In [1]: fig = figure()
In [2]: ax1 = fig.add_subplot(221)
In [3]: ax2 = fig.add_subplot(222)
In [4]: ax3 = fig.add_subplot(223)
In [5]: ax4 = fig.add_subplot(224)
In [6]: for ax in [ax1, ax2, ax3, ax4]:
...:
ax.plot([3,5,8],[6,3,1])
In [7]: show()
95
5.5 Exercises
1. Plot the function y = 3x2 for 1 x 3 as a continuous line. Include enough
points so that the curve you plot appears smooth. Label the axes x and y.
2. Plot the following function for 15 x 15:
y=
cos x
1 + 15 x2
Include enough points so that the curve you plot appears smooth. Label the axes x
and y.
3. Plot the functions sin x and cos x vs x on the same plot with x going from to
. Make sure the limits of x-axis do not extend beyond the limits of the data. Plot
sin x in the color green and cos x in the color black and include a legend to label the
two curves. Place the legend within the plot, but such that it does not cover either
of the sine or cosine traces.
4. Create a data file with the data shown below.
(a) Read the data into Python program and plot t vs y using circles for data points
with error bars. Use the data in the dy column as the error estimates for the y
data. Label the horizontal and vertical axes time (s) and position (cm).
(b) On the same graph, plot the function below as a smooth line. Make the line
pass behind the data points.
1
t
y(t) = 3 + sin
t et/10
2
5
96
d
2.94
8.29
9.36
11.60
9.32
7.75
8.06
5.60
dy
0.7
1.2
1.2
1.4
1.3
1.1
1.2
1.0
Chapter 5. Plotting
29.0
32.5
36.0
39.5
43.0
4.50
4.01
2.62
1.70
2.03
0.8
0.8
0.7
0.6
0.6
force N
Your plot should also include a visual straight best fit to the data as well as visual
fits that give the smallest and largest slopes consistent with the data. Note, you
only need two points to define a straight line so the straight lines you draw on the
plot should be arrays of length 2 and no longer. All of your fitted lines should lie
behind the data. Try to make your plot look like the one below. In addition, add
a legend to your plot the gives the slope with its uncertainty obtained from your
visual fits to the data.
distance cm
The web page http://matplotlib.org/api/pyplot_summary.html gives a summary of
the main plotting commands available in MatPlotLib. The two important ones here
5.5. Exercises
97
are plot and errorbar, which make regular plots and plots with error bars,
respectively. You will find the following keyword arguments useful: yerr, ls,
marker, mfc, mec, ms, and ecolor, which you can find described by clicking
on the errorbar function link on the web page cited above.
7. The data file below shows data obtained for the displacement (position) vs time of
a falling object, together with the estimated uncertainty in the displacement.
Measurements of fall velocity vs time
Taken by A.P. Crawford and S.M. Torres
19-Sep-13
time (s)
position (m)
uncertainty (m)
0.0
0.0
0.04
0.5
1.3
0.12
1.0
5.1
0.2
1.5
10.9
0.3
2.0
18.9
0.4
2.5
28.7
0.4
3.0
40.3
0.5
3.5
53.1
0.6
4.0
67.5
0.6
4.5
82.3
0.6
5.0
97.6
0.7
5.5
113.8
0.7
6.0
131.2
0.7
6.5
148.5
0.7
7.0
166.2
0.7
7.5
184.2
0.7
8.0
201.6
0.7
8.5
220.1
0.7
9.0
238.3
0.7
9.5
256.5
0.7
10.0
275.6
0.8
(a) Use these data to calculate the velocity and acceleration (in a Python program
.py file), together with their uncertainties propagated from the displacement
vs time uncertainties. Be sure to calculate time arrays corresponding the midpoint in time between the two displacements or velocities for the velocity and
acceleration arrays, respectively.
(b) In a single window frame, make three vertically stacked plots of the displacement, velocity, and acceleration vs time. Show the error bars on the different
plots. Make sure that the time axes of all three plots cover the same range of
times. Why do the relative sizes of the error bars grow progressively greater
as one progresses from displacement to velocity to acceleration?
98
Chapter 5. Plotting
CHAPTER
SIX
6.1 Conditionals
Conditional statements allow a computer program to take different actions based on
whether some condition, or set of conditions is true or false. In this way, the programmer
can control the flow of a program.
99
4
5
d = b*b - 4.*a*c
6
7
8
9
10
11
12
if d >= 0.0:
print("Solutions are real")
elif b == 0.0:
print("Solutions are imaginary")
else:
print("Solutions are complex")
# block 1
# block 2
# block 3
13
14
print("Finished!")
After getting the inputs of from the user, the program evaluates the discriminant d. The
code d >= 0.0 has a boolean truth value of either True or False depending on
whether or not d 0. You can check this out in the interactive IPython shell by typing the following set of commands
In [2]: d = 5
In [3]: d >= 2
Out[3]: True
In [4]: d >= 7
Out[4]: False
100
Therefore, the if statement in line 7 is simply testing to see if the statement d >= 0.0
if True or False. If the statement is True, Python executes the indented block of
statements following the if statement. In this case, there is only one line in indented
block. Once it executes this statement, Python skips past the elif and else blocks and
executes the print("Finished!") statement.
If the if statement in line 7 is False, Python skips the indented block directly below the if statement and executes the elif statement. If the condition b == 0.0
is True, it executes the indented block immediately below the elif statement and
then skips the else statement and the indented block below it. It then executes the
print("Finished!") statement.
Finally, if the elif statement is False, Python skips to the else statement and executes the block immediately below the else statement. Once finished with that indented
block, it then executes the print("Finished!") statement.
As you can see, each time a False result is obtained in an if or elif statement, Python
skips the indented code block associated with that statement and drops down to the next
conditional statement, that is, the next elif or else. A flowchart of the if-elif-else code
is shown below.
start
if
d0?
block 1
True
False
elif
b=0?
True
block 2
False
block 3
else
finish
101
however, we only check to see if b = 0. The reason that we dont have to check if d < 0
is that the elif statement is executed only if the condition if d >= 0.0 on line 7
is False. Similarly, we dont have to check if if b = 0 and d < 0 for the final else
statement because this part of the if, elif, and else block will only be executed if
the preceding if and elif statements are False. This illustrates a key feature of the
if, elif, and else statements: these statements are executed sequentially until one of
the if or elif statements is found to be True. Therefore, Python reaches an elif or
else statement only if all the preceding if and elif statements are False.
The if-elif-else logical structure can accomodate as many elif blocks as desired.
This allows you to set up logic with more than the three possible outcomes illustrated in
the example above. When designing the logical structure you should keep in mind that
once Python finds a true condition, it skips all subsequent elif and else statements in
a given if, elif, and else block, irrespective of their truth values.
if-else example
You will often run into situations where you simply want the program to execute one of
two possible blocks based on the outcome of an if statement. In this case, the elif
block is omitted and you simply use an if-else structure. The following program testing whether an integer is even or odd provides a simple example.
a = int(raw_input("Please input an integer: "))
if a%2 == 0:
print("{0:0d} is an even number.".format(a))
else:
print("{0:0d} is an odd number.".format(a))
102
start
if
d0?
True
block 1
False
else
block 2
finish
This works exactly as the preceding code. Note, however, that if the block of code associated with an if or elif statement is more than one line long, the entire block of code
should be written as indented text below the if or elif statement.
The flowchart below shows the logical structure of a simple if structure.
start
if
d0?
True
block 1
False
finish
6.1. Conditionals
103
function
less than
less than or equal to
greater than
greater than or equal to
equal
not equal
both must be true
one or both must be true
reverses the truth value
# a is not equal to 5
104
In [11]: a>2
Out[11]: True
In [12]: not a>2
Out[12]: False
Logical statements like those above can be used in if, elif, and, as we shall see below,
while statements, according to your needs.
6.2 Loops
In computer programming a loop is statement or block of statements that is executed
repeatedly. Python has two kinds of loops, a for loop and a while loop. We first
introduce the for loop and illustrate its use for a variety of tasks. We then introduce the
while loop and, after a few illustrative examples, compare the two kinds of loops and
discuss when to use one or the other.
6.2. Loops
105
Buster
Arf, arf!
Maggie
Arf, arf!
Lucy
Arf, arf!
All done.
The for loop works as follows: the iteration variable or loop index dogname is set
equal to the first element in the list, "Max", and then the two lines in the indented body
are executed. Then dogname is set equal to second element in the list, "Molly", and the
two lines in the indented body are executed. The loop cycles through all the elements of
the list, and then moves on to the code that follows the for loop and prints All done.
When indenting a block of code in a Python for loop, it is critical that every line be
indented by the same amount. Using the <tab> key causes the Code Editor to indent 4
spaces. Any amount of indentation works, as long as it is the same for all lines in a for
loop. While code editors designed to work with Python (including Canopy and Spyder)
translate the <tab> key to 4 spaces, not all text editors do. In those cases, 4 spaces are
not equivalent to a <tab> character even if they appear the same on the display. Indenting
some lines by 4 spaces and other lines by a <tab> character will produce an error. So
beware!
The figure below shows the flowchart for a for loop. It starts with an implicit conditional
asking if there are any more elements in the sequence. If there are, it sets the iteration
variable equal to the next element in the sequence and then executes the bodythe indented textusing that value of the iteration variable. It then returns to the beginning to
see if there are more elements in the sequence and continues the loop until there is none
remaining.
6.2.2 Accumulators
Lets look at another application of Pythons for loop. Suppose you want to calculate
the sum of all the odd numbers between 1 and 100. Before writing a computer program
to do this, lets think about how you would do it by hand. You might start by adding
1+3=4. Then take the result 4 and add the next odd integer, 5, to get 4+5=9; then 9+7=16,
then 16+9=25, and so forth. You are doing repeated additions, starting with 1+3, while
keeping track of the running sum, until you reach the last number 99.
In developing an algorithm for having the computer sum the series of numbers, we are
going to do the same thing: add the numbers one at a time while keeping track of the
running sum, until we reach the last number. We will keep track of the running sum with
the variable s, which is called the accumulator. Initially s=0, since we havent adding
106
start
more items
in seq?
No
Yes
<intervar> =
next item in seq
<body>
finish
s = 0
for i in range(1, 100, 2):
s = s+i
print(s)
The range function defines the list [1, 3, 5, ..., 97, 99]. The for loop
successively adds each number in the list to the running sum until it reaches the last
element in the list and the sum is complete. Once the for loop finishes, the program
exits the loop and the final value of s, which is the sum of the odd numbers from 1 to 99,
is printed out. Copy the above program and run it. You should get an answer of 2500.
6.2. Loops
107
while <condition>:
<body>
where <condition> is a statement that can be either True or False and <body> is
a series of Python commands that is executed repeatedly until <condition> becomes
false. This means that somewhere in <body>, the truth value of <condition> must be
changed so that it becomes false after a finite number of iterations. Consider the following
example.
Suppose you want to calculate all the Fibonacci numbers smaller than 1000. The Fibonacci numbers are determined by starting with the integers 0 and 1. The next number
in the sequence is the sum of the previous two. So, starting with 0 and 1, the next Fibonacci number is 0 + 1 = 1, giving the sequence 0, 1, 1. Continuing this process, we
obtain 0, 1, 1, 2, 3, 5, 8, ... where each element in the list is the sum of the previous two.
Using a for loop to calculate the Fibonacci numbers is impractical because we do not
know ahead of time how many Fibonacci numbers there are smaller than 1000. By contrast a while loop is perfect for calculating all the Fibonacci numbers because it keeps
calculating Fibonacci numbers until it reaches the desired goal, in this case 1000. Here is
the code using a while loop.
x, y = 0, 1
while x < 1000:
print(x)
x, y = y, x+y
We have used the multiple assignment feature of Python in this code. Recall that all the
values on the left are assigned using the original values of x and y.
The figure below shows the flowchart for the while loop. The loop starts with the evaluation of a condition. If the condition is False, the code in the body is skipped, the flow
exits the loop, and then continues with the rest of the program. If the condition is True,
the code in the bodythe indented textis executed. Once the body is finished, the flow
returns to the condition and proceeds along the True or False branches depending on
the truth value of the condition. Implicit in this loop is the idea that somewhere during
the execution of the body of the while loop, the variable that is evaluated in the condition
is changed in some way. Eventually that change will cause the condition to return a value
of False so that the loop will end.
One danger of a while loop is that it entirely possible to write a loop that never
terminatesan infinite loop. For example, if we had written while y > 0:, in place
of while x < 1000:, the loop would never end. If you execute code that has an infinite loop, you can often terminate the program from the keyboard by typing ctrl-C a
couple of times. If that doesnt work, you may have to terminate and then restart Python.
108
start
<condition>
False
True
<body>
(change condition)
finish
Running this on my computer returns the result in about 8 secondsnot bad for having
performed 10 million multiplications. Of course we could have performed the same calculation using the array multiplication we learned in Chapter 3 (Strings, Lists, Arrays,
and Dictionaries). Here is the code.
6.2. Loops
109
import numpy as np
a = np.linspace(0, 32, 1e7)
print(a)
a = a*a
print(a)
Running this on my computer returns the results faster than I can discern, but certainly
much less than a second. This illustrates an important point: for loops are slow. Array
operations run much faster and are therefore to be preferred in any case where you have
a choice. Sometimes finding an array operation that is equivalent to a loop can be difficult, especially for a novice. Nevertheless, doing so pays rich rewards in execution time.
Moreover, the array notation is usually simpler and clearer, providing further reasons to
prefer array operations over loops.
Suppose we want to construct a vector from the diagonal elements of this matrix. We
could do so with a for loop with an accumulator as follows
In [2]: diag = []
...: for i in [0, 1, 2]:
...:
diag.append(x[i][i])
...:
In [3]: diag
Out[3]: [1, 5, 9]
List comprehensions provide a simpler, cleaner, and faster way of doing the same thing
In [4]: diagLC = [x[i][i] for i in [0, 1, 2]]
In [5]: diagLC
Out[5]: [1, 5, 9]
110
Notice here how y serves as a dummy variable accessing the various elements of the list
diagLC.
Extracting a column from a 2-dimensaional array such as x is quite easy. For example the
second row is obtained quite simply in the following fashion
In [7]: x[1]
Out[7]: [4, 5, 6]
Obtaining a column is not as simple, but a list comprehension makes it quite straightforward:
In [7]: c1 = [a[1] for a in x]
In [8]: c1
Out[8]: [2, 5, 8]
Suppose you have a list of numbers and you want to extract all the elements of the list that
are divisible by three. A slightly fancier list comprehension accomplishes the task quite
simply and demonstrates a new feature:
In [10]: y = [-5, -3, 1, 7, 4, 23, 27, -9, 11, 41]
In [14]: [a for a in y if a%3==0]
Out[14]: [-3, 27, -9]
111
6.4 Exercises
1. Write a program to calculate the factorial of a positive integer input by the user.
Recall that the factorial function is given by x! = x(x 1)(x 2)...(2)(1) so that
1! = 1, 2! = 2, 3! = 6, 4! = 24, ...
(a) Write the factorial function using a Python while loop.
(b) Write the factorial function using a Python for loop.
Check your programs to make sure they work for 1, 2, 3, 5, and beyond, but especially for the first 5 integers.
2. The following Python program finds the smallest non-trivial (not 1) prime factor of
a positive integer.
n = int(raw_input("Input an integer > 1: "))
i = 2
while (n % i) != 0:
i += 1
print("The smallest factor of n is:", i )
(a) Type this program into your computer and verify that it works as advertised.
Then briefly explain how it works and why the while loop always terminates.
(b) Modify the program so that it tells you if the integer input is a prime number
or not. If it is not a prime number, write your program so that it prints out the
smallest prime factor. Using your program verify that the following integers
are prime numbers: 101, 8191, 947431.
3. Consider the matrix list x = [[1, 2, 3], [4, 5, 6], [7, 8, 9]].
Write a list comprehension to extract the last column of the matrix [3, 6, 9]. Write
another list comprehension to create a vector of twice the square of the middle
column [8, 50, 128].
4. Write a program that calculates the value of an investment after some number of
years specified by the user if
(a) the principal is compounded annually
(b) the principle is compounded monthly
(c) the principle is compounded daily
112
Your program should ask the user for the initial investment (principal), the interest
rate in percent, and the number of years the money will be invested (allow for
fractional years). For an initial investment of $1000 at an interest rate of 6%, after
10 years I get $1790.85 when compounded annually, $1819.40 when compounded
monthly, and $1822.03 when compounded daily, assuming 12 months in a year and
365.24 days in a year, where the monthly interest rate is the annual rate divided
by 12 and the daily rate is the annual rate divided by 365 (dont worry about leap
years).
5. Write a program that determines the day of the week for any given calendar date
after January 1, 1900, which was a Monday. Your program will need to take into
account leap years, which occur in every year that is divisible by 4, except for years
that are divisible by 100 but are not divisible by 400. For example, 1900 was not
a leap year, but 2000 was a leap year. Test that your program gives the following
answers: Monday 1900 January 1, Tuesday 1933 December 5, Wednesday 1993
June 23, Thursday 1953 January 15, Friday 1963 November 22, Saturday 1919
June 28, Sunday 2005 August 28.
6.4. Exercises
113
114
CHAPTER
SEVEN
FUNCTIONS
As you develop more complex computer code, it becomes increasingly important to organize your code into modular blocks. One important means for doing so is user-defined
Python functions. User-defined functions are a lot like built-in functions that we have
encountered in core Python as well as in NumPy and Matplotlib. The main difference is
that user-defined functions are written by you. The idea is to define functions to simplify
your code and to allow you to reuse the same code in different contexts.
The number of ways that functions are used in programming is so varied that we cannot
possibly enumerate all the possibilities. As our use of Python functions in scientific program is somewhat specialized, we introduce only a few of the possible uses of Python
functions, ones that are the most common in scientific programming.
sin x
.
x
Lets write a Python function for the sinc function. Here is our first attempt:
115
def sinc(x):
y = np.sin(x)/x
return y
Every function definition begins with the word def followed by the name you want to
give to the function, sinc in this case, then a list of arguments enclosed in parentheses,
and finally terminated with a colon. In this case there is only one argument, x, but in
general there can be as many arguments as you want, including no arguments at all. For
the moment, we will consider just the case of a single argument.
The indented block of code following the first line defines what the function does. In
this case, the first line calculates sinc x = sin x/x and sets it equal to y. The return
statement of the last line tells Python to return the value of y to the user.
We can try it out in the IPython shell. First we type in the function definition.
In [1]: def sinc(x):
...:
y = sin(x)/x
...:
return y
Because we are doing this from the IPython shell, we dont need to import NumPy; its
preloaded. Now the function sinc x is available to be used from the IPython shell
In [2]: sinc(4)
Out[2]: -0.18920062382698205
In [3]: a = sinc(1.2)
In [4]: a
Out[4]: 0.77669923830602194
In [5]: sin(1.2)/1.2
Out[5]: 0.77669923830602194
Inputs and outputs 4 and 5 verify that the function does indeed give the same result as an
explicit calculation of sin x/x.
You may have noticed that there is a problem with our definition of sinc x when x=0.0.
Lets try it out and see what happens
In [6]: sinc(0.0)
Out[6]: nan
IPython returns nan or not a number, which occurs when Python attempts a division by
zero, which is not defined. This is not the desired response as sinc x is, in fact, perfectly
well defined for x = 0. You can verify this using LHopitals rule, which you may have
116
Chapter 7. Functions
learned in your study of calculus, or you can ascertain the correct answer by calculating
the Taylor series for sinc x. Here is what we get
sinc x =
x
sin x
=
x
x3
3!
+
x
x5
5!
+ ...
=1
x2
x4
+
+ ... .
3!
5!
From the Taylor series, it is clear that sinc x is well-defined at and near x = 0 and that, in
fact, sinc(0) = 1. Lets modify our function so that it gives the correct value for x=0.
In [7]: def sinc(x):
...:
if x==0.0:
...:
y = 1.0
...:
else:
...:
y = sin(x)/x
...:
return y
In [8]: sinc(0)
Out[8]: 1.0
In [9]: sinc(1.2)
Out[9]: 0.77669923830602194
Now our function gives the correct value for x=0 as well as for values different from zero.
0.5, 1. ,
4.5])
1.5,
2. ,
2.5,
3. ,
3.5,
In [12]: sinc(x)
---------------------------------------------------------ValueError
Traceback (most recent call last)
----> 1 sinc(x)
1 def sinc(x):
----> 2
if x==0.0:
3
y = 1.0
117
4
5
else:
y = np.sin(x)/x
The if statement in Python is set up to evaluate the truth value of a single variable, not of
multielement arrays. When Python is asked to evaluate the truth value for a multi-element
array, it doesnt know what to do and therefore returns an error.
An obvious way to handle this problem is to write the code so that it processes the array
one element at a time, which you could do using a for loop, as illustrated below.
1
2
3
4
5
6
7
8
def sinc(x):
y = []
# creates an empty list to store results
for xx in x:
# loops over all elements in x array
if xx==0.0:
# adds result of 1.0 to y list if
y += [1.0] # xx is zero
else:
# adds result of sin(xx)/xx to y list if
y += [np.sin(xx)/xx] # xx is not zero
return np.array(y) # converts y to array and returns array
9
10
11
import numpy as np
import matplotlib.pyplot as plt
12
13
14
15
16
17
18
19
plt.plot(x, y)
plt.axhline(color="gray", zorder=-1)
plt.axvline(color="gray", zorder=-1)
plt.show()
The for loop evaluates the elements of the x array one by one and appends the results to
the list y one by one. When it is finished, it converts the list to an array and returns the
array. The code following the function definition plots sinc x as a function of x.
In the program above, you may have noticed that the NumPy library is imported after
the sinc(x) function definition. As the function uses the NumPy functions sin and
array, you may wonder how this program can work. Doesnt the import numpy
statement have to be called before any NumPy functions are used? The answer it an emphatic YES. What you need to understand is that the function definition is not executed
when it is defined, nor can it be as it has no input x data to process. That part of the code
is just a definition. The first time the code for the sinc(x) function is actually executed
is when it is called on line 14 of the program, which occurs after the NumPy library is
118
Chapter 7. Functions
imported in line 10. The figure below shows the plot of the sinc x function generated by
the above code.
1.0
0.8
0.6
0.4
0.2
0.0
0.2
0.410
10
The first argument of the where function is a conditional statement involving an array.
The where function applies the condition to the array element by element, and returns the
second argument for those array elements for which the condition is True, and returns the
third argument for those array elements that are False. We can apply it to the sinc(x)
function as follows
7.1. User-defined functions
119
def sinc(x):
z = np.where(x==0.0, 1.0, np.sin(x)/x)
return z
The where function creates the array y and sets the elements of y equal to 1.0 where the
corresponding elements of x are zero, and otherwise sets the corresponding elements to
sin(x)/x. This code executes much faster, 25 to 100 times, depending on the size of
the array, than the code using a for loop. Moreover, the new code is much simpler to
write and read. An additional benefit of the where function is that it can handle single
variables and arrays equally well. The code we wrote for the sinc function with the for
loop cannot handle single variables. Of course we could rewrite the code so that it did,
but the code becomes even more clunky. Its better just to use NumPys where function.
The moral of the story
The moral of the story is that you should avoid using for and while loops to process
arrays in Python programs whenever an array-processing method is available. As a beginning Python programmer, you may not always see how to avoid loops, and indeed,
avoiding them is not always possible, but you should look for ways to avoid loops, especially loops that iterate a large number of times. As you become more experienced,
you will find that using array-processing methods in Python becomes more natural. Using them can greatly speed up the execution of your code, especially when working with
large arrays.
7.1.3 Functions with more (or less) than one input or output
Python functions can have any number of input arguments and can return any number of
variables. For example, suppose you want a function that outputs n (x, y) coordinates
around a circle of radius r centered at the point (x0 , y0 ). The inputs to the function would
be r, x0 , y0 , and n. The outputs would be the n (x, y) coordinates. The following code
implements this function.
def circle(r, x0, y0, n):
theta = np.linspace(0., 2.*np.pi, n, endpoint=False)
x = r * np.cos(theta)
y = r * np.sin(theta)
return x0+x, y0+y
This function has four inputs and two outputs. In this case, the four inputs are simple
numeric variables and the two outputs are NumPy arrays. In general, the inputs and
120
Chapter 7. Functions
outputs can be any combination of data types: arrays, lists, strings, etc. Of course, the
body of the function must be written to be consistent with the prescribed data types.
Functions can also return nothing to the calling program but just perform some task. For
example, here is a program that clears the terminal screen
import subprocess
import platform
def clear():
subprocess.Popen( "cls" if platform.system() ==
"Windows" else "clear", shell=True)
The function is invoked by typing clear(). It has no inputs and no outputs but it
performs a useful task. This function uses two standard Python libraries, subprocess
and platform that are useful for performing computer system tasks. Its not important
that you know anything about them at this point. We simply use them here to demonstrate
a useful cross-platform function that has no inputs and returns no values.
The default values of the arguments x0, y0, and n are specified in the argument of the
function definition in the def line. Arguments whose default values are specified in this
manner are called keyword arguments, and they can be omitted from the function call if
the user is content using those values. For example, writing circle(4) is now a perfectly legal way to call the circle function and it would produce 12 (x, y) coordinates
centered about the origin (x, y) = (0, 0). On the other hand, if you want the values of x0,
121
y0, and n to be something different from the default values, you can specify their values
as you would have before.
If you want to change only some of the keyword arguments, you can do so by using the
keywords in the function call. For example, suppose you are content with have the circle
centered on (x, y) = (0, 0) but you want only 6 points around the circle rather than 12.
Then you would call the circle function as follows:
circle(2, n=6)
The unspecified keyword arguments keep their default values of zero but the number of
points n around the circle is now 6 instead of the default value of 12.
The normal arguments without keywords are called positional arguments; they have to
appear before any keyword arguments and, when the function is called, must be supplied
values in the same order as specified in the function definition. The keyword arguments,
if supplied, can be supplied in any order providing they are supplied with their keywords.
If supplied without their keywords, they too must be supplied in the order they appear
in the function definition. The following function calls to circle both give the same
output.
In [13]: circle(3, n=3, y0=4, x0=-2)
Out[13]: (array([ 1. , -3.5, -3.5]),
array([ 4.
, 6.59807621,
1.40192379]))
By now you probably have noticed that we used the keyword argument endpoint in
calling linspace in our definition of the circle function. The default value of
endpoint is True, meaning that linspace includes the endpoint specified in the
second argument of linspace. We set it equal to False so that the last point was not
included. Do you see why?
122
Chapter 7. Functions
p = 1
for num in args:
p *= num
return p
In [15]: product(11., -2, 3)
args = (11.0, -2, 3)
Out[15]: -66.0
In [16]: product(2.31, 7)
args = (2.31, 7)
Out[16]: 16.17
Now lets find the derivative of f0 (x) = 4x5 at x = 3 using the function deriv:
In [17]: deriv(f0, 3)
Out[17]: 1620.0001482502557
The exact result is 1620, so our function to numerically calculate the derivative works
pretty well.
123
Suppose we want to calculate the derivative of this function for a particular set of parameters a and p. Now we face a problem, because it might seem that there is no way to pass
the parameters a and p to the deriv function. Moreover, this is a generic problem for
functions such as deriv that use a function as an input, because different functions you
want to use as inputs generally come with different parameters. Therefore, we would like
to write our program deriv so that it works, irrespective of how many parameters are
needed to specify a particular function.
This is what the optional positional argument *params defined in deriv is for: to pass
parameters of f1, like a and b, through deriv. To see how this works, lets set a and b
to be 4 and 5, respectively, the same values we used in the definition of f0, so that we can
compare the results:
In [16]: deriv(f1, 3, 1.e-9, 4, 5)
Out[16]: 1620.0001482502557
We get the same answer as before, but this time we have used deriv with a more general
form of the function f1 (x) = axp .
The order of the parameters is important. The function deriv uses x, the first argument
of f1, as its principal argument, and then uses a and p, in the same order that they are
defined in the function f1, to fill in the additional argumentsthe parametersof the
function f1.
Optional arguments must appear after the regular positional and keyword arguments in a
function call. The order of the arguments must adhere to the following convention:
def func(pos1, pos2, ..., keywd1, keywd2, ..., *args, **kwargs):
That is, the order of arguments is: positional arguments first, then keyword arguments, then optional positional arguments (*args), then optional keyword arguments
(**kwargs). Note that to use the *params argument, we had to explicitly include the
keyword argument h even though we didnt need to change it from its default value.
Python also allows for a variable number of keyword arguments**kwargsin a function call. While *args is a tuple, kwargs is a dictionary, so that the value of an optional
keyword argument is accessed through its dictionary key.
124
Chapter 7. Functions
def sinc(x):
z = np.where(x==0.0, 1.0, np.sin(x)/x)
return z
4
5
6
import numpy as np
import matplotlib.pyplot as plt
7
8
9
10
11
12
13
14
plt.plot(x, y)
plt.axhline(color="gray", zorder=-1)
plt.axvline(color="gray", zorder=-1)
plt.show()
Running this program produces a plot like the plot of sinc shown in the previous section.
Notice that the array variable z is only defined within the function definition of sinc. If
we run the program from the IPython terminal, it produces the plot, of course. Then if we
ask IPython to print out the arrays, x, y, and z, we get some interesting and informative
results, as shown below.
7.1. User-defined functions
125
, -9.99969482,
9.99969482, 10.
-9.99938964, ...,
])
In [17]: y
Out[17]: array([-0.05440211, -0.05437816, -0.0543542 , ...,
-0.0543542 , -0.05437816, -0.05440211])
In [18]: z
--------------------------------------------------------NameError
Traceback (most recent call last)
NameError: name z is not defined
When we type in x at the In [16]: prompt, IPython prints out the array x (some of the
output is suppressed because the array x has many elements); similarly for y. But when
we type z at the In [18]: prompt, IPython returns a NameError because z is not
defined. The IPython terminal is working in the same namespace as the program. But the
namespace of the sinc function is isolated from the namespace of the program that calls
it, and therefore isolated from IPython. This also means that when the sinc function ends
with return z, it doesnt return the name z, but instead assigns the values in the array
z to the array y, as directed by the main program in line 9.
Passing variables and arrays to functions: mutable and immutable objects
What happens to a variable or an array passed to a function when the variable or array
is changed within the function? It turns out that the answers are different depending on
whether the variable passed is a simple numeric variable, string, or tuple, or whether it
is an array or list. The program below illustrates the different ways that Python handles
single variables vs the way it handles lists and arrays.
1
2
3
4
5
6
7
8
9
import numpy as np
10
126
Chapter 7. Functions
11
12
13
14
15
s
v
t
l
a
=
=
=
=
=
16
17
18
19
20
21
22
23
24
25
print(*************)
print("s = {0:s}".format(s))
print("v = {0:5.2f}".format(v))
print("t = {0:s}".format(t))
print("l = {0:s}".format(l))
print("a = "),
# comma suppresses line feed
print(a)
print(*************)
print(*call "test"*)
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
print(*************)
print("s1 = {0:s}".format(s1))
print("v1 = {0:5.2f}".format(v1))
print("t1 = {0:s}".format(t1))
print("l1 = {0:s}".format(l1))
print("a1 = "),
print(a1)
print(*************)
print("s = {0:s}".format(s))
print("v = {0:5.2f}".format(v))
print("t = {0:s}".format(t))
print("l = {0:s}".format(l))
print("a = "),
# comma suppresses line feed
print(a)
print(*************)
The function test has five arguments, a string s, a numerical variable v, a tuple t, a list
l, and a NumPy array a. test modifies each of these arguments and then returns the
modified s, v, t, l, a. Running the program produces the following output.
In [17]: run passingVars.py
*************
s = How do you do?
v = 5.00
t = (97.5, 82.9, 66.7)
l = [3.9, 5.7, 7.5, 9.3]
a = [ 3.9 5.7 7.5 9.3]
127
*************
*call "test"*
*************
s1 = I am doing fine
v1 = 9.87
t1 = (1.1, 2.9)
l1 = [3.9, 5.7, 7.5, end]
a1 = [ 963.2
5.7
7.5
*************
s = How do you do?
v = 5.00
t = (97.5, 82.9, 66.7)
l = [3.9, 5.7, 7.5, end]
a = [ 963.2
5.7
7.5
*************
9.3]
9.3]
The program prints out three blocks of variables separated by asterisks. The first block
merely verifies that the contents of s, v, t, l, and a are those assigned in lines 10-13.
Then the function test is called. The next block prints the output of the call to the
function test, namely the variables s1, v1, t1, l1, and a1. The results verify that the
function modified the inputs as directed by the test function.
The third block prints out the variables s, v, t, l, and a from the calling program after
the function test was called. These variables served as the inputs to the function test.
Examining the output from the third printing block, we see that the values of the string s,
the numeric variable v, and the contents of t are unchanged after the function call. This
is probably what you would expect. On the other hand, we see that the list l and the array
a are changed after the function call. This might surprise you! But these are important
points to remember, so important that we summarize them in two bullet points here:
Changes to string, variable, and tuple arguments of a function within the function
do not affect their values in the calling program.
Changes to values of elements in list and array arguments of a function within the
function are reflected in the values of the same list and array elements in the calling
function.
The point is that simple numerics, strings and tuples are immutable while lists and arrays
are mutable. Because immutable objects cant be changed, changing them within a function creates new objects with the same name inside of the function, but the old immutable
objects that were used as arguments in the function call remain unchanged in the calling
program. On the other hand, if elements of mutable objects like those in lists or arrays are
changed, then those elements that are changed inside the function are also changed in the
calling program.
128
Chapter 7. Functions
Any object in Python can and in general does have a number of attributes that are accessed
in just the way demonstrated above, with a period and the attribute name following the
name of the particular object. In general, attributes involve properties of the object that
are stored by Python with the object and require no computation. Python just looks up the
attribute and returns its value.
Objects in Python also have associated with them a number of specialized functions called
methods that act on the object. In contrast to attributes, methods generally involve Python
performing some kind of computation. Methods are accessed in a fashion similar to
attributes, by appending a period followed the methods name, which is followed by a
pair of open-close parentheses, consistent with methods being a kind of function that acts
on the object. Often methods are used with no arguments, as methods by default act on the
object whose name they follow. In some cases. however, methods can take arguments.
Examples of methods for NumPy arrays are sorting, calculating the mean, or standard
deviation of the array. The code below illustrates a few array methods.
129
In [21]: a
Out[21]:
array([ 0.859057 ,
0.05067356,
0.31527767,
0.27228037, 0.87780026,
0.83490135, 0.54844515,
0.15868803])
0.14341207,
0.33583966,
In [22]: a.sum()
Out[22]: 4.3963751104791005
# sum
In [23]: a.mean()
Out[23]: 0.43963751104791005
# mean or average
In [24]: a.var()
Out[24]: 0.090819477333711512
# variance
In [25]: a.std()
Out[25]: 0.30136270063448711
# standard deviation
In [26]: a.sort()
In [27]: a
Out[27]:
array([ 0.05067356,
0.31527767,
0.859057 ,
In [28]: a.clip(0.3,
Out[29]:
array([ 0.3
,
0.31527767,
0.8
,
0.14341207, 0.15868803,
0.33583966, 0.54844515,
0.87780026])
0.27228037,
0.83490135,
0.8)
0.3
, 0.3
,
0.33583966, 0.54844515,
0.8
])
0.3
0.8
,
,
The clip() method provides an example of a method that takes an argument, in this
case the arguments are the lower and upper values to which array elements are cutoff if
their values are outside the range set by these values.
Chapter 7. Functions
f (x; a, b, c, ...), where x is the independent variable and a, b, c, ... are parameters to be
adjusted so that the function f (x; a, b, c, ...) best fits the experimental data. For example,
suppose we had some data of the velocity vs time for a falling mass. If the mass falls
only a short distance such that its velocity remains well below its terminal velocity, we
can ignore air resistance. In this case, we expect the acceleration to be constant and the
velocity to change linearly in time according to the equation
v(t) = v0 gt ,
(7.1)
where g is the local gravitational acceleration. We can fit the data graphically, say by
plotting it as shown below in Fig. 4.6 and then drawing a line through the data. When
we draw a straight line through a data, we try to minimize the distance between the points
and the line, globally averaged over the whole data set.
5
0
velocity (m/s)
5
10
15
20
25
0.0
0.5
1.0
1.5
time (s)
2.0
2.5
3.0
n
X
(7.2)
131
where yi and f (xi ; a, b, c, ...) are the values of the experimental data and the fitting function, respectively, at xi , and S is the square of their difference summed over all n data
points. The quantity S is a sort of global measure of how much the the fit f (xi ; a, b, c, ...)
differs from the experimental data yi .
Notice that for a given set of data points {xi , yi }, S is a function only of the fitting
parameters a, b, ..., that is, S = S(a, b, c, ...). One way of defining a best fit, then, is to
find the values of the fitting parameters a, b, ... that minimize the S.
In principle, finding the values of the fitting parameters a, b, ... that minimize the S is a
simple matter. Just set the partial derivatives of S with respect to the fitting parameter
equal to zero and solve the resulting system of equations:
S
=0,
a
S
= 0 , ...
b
(7.3)
Because there are as many equations as there are fitting paramters, we should be able to
solve the system of equations and find the values of the fitting parameters that minimize S.
Solving those systems of equations is straightforward if the fitting function f (x; a, b, ...)
is linear in the fitting parameters. Some examples of fitting functions linear in the fitting
parameters are:
f (x; a, b) = a + bx
f (x; a, b, c) = a + bx + cx2
x
(7.4)
x2
f (x; a, b, c) = a sin x + be + ce
For fitting functions such as these, taking the partial derivatives with respect to the fitting
parameters, as proposed in (7.3), results in a set of algebraic equations that are linear in
the fitting paramters a, b, ... Because they are linear, these equations can be solved in a
straightforward manner.
For cases in which the fitting function is not linear in the fitting parameters, one can
generally still find the values of the fitting parameters that minimize S but finding them
requires more work, which goes beyond our immediate interests here.
132
Chapter 7. Functions
Finding the best fit in this case corresponds to finding the values of the fitting parameters
a and b for which S(a, b) is a minimum. To find the minimum, we set the derivatives of
S(a, b) equal to zero:
!
X
X
X
S
=
2(yi a bxi ) = 2 na + b
xi
yi = 0
a
i
i
i
!
(7.6)
X
X
X
X
S
2
=
2(yi a bxi ) xi = 2 a
xi + b
xi
xi yi = 0
b
i
i
i
i
Dividing both equations by 2n leads to the equations
a + b
x = y
1X
1X 2
xi =
xi yi
a
x+b
n i
n i
(7.7)
1X
xi
n i
1X
y =
yi .
n i
(7.8)
where
x
=
n
x
i i
a = y b
x.
P
P
Noting that n
y = i y and n
x = i x, the results can be written as
P
(xi x
) yi
b = Pi
) xi
i (xi x
(7.9)
(7.10)
a = y b
x.
While Eqs. (7.9) and (7.10) are equivalent analytically, Eq. (7.10) is preferred for numerical calculations because Eq. (7.10) is less sensitive to roundoff errors. Here is a Python
function implementing this algorithm:
1
2
3
133
xavg = x.mean()
slope = (y*(x-xavg)).sum()/(x*(x-xavg)).sum()
yint = y.mean()-slope*xavg
return slope, yint
4
5
6
7
X yi a bxi 2
i
(7.12)
Finding the minimum for 2 (a, b) follows the same procedure used for finding the minimum of S(a, b) in the previous section. The result is
P
(xi x
) yi /i2
b = Pi
) xi /i2
(7.13)
i (xi x
a = y b
x.
134
Chapter 7. Functions
where
P
xi /i2
x
= Pi
1/i2
Pi
yi /i2
y = Pi
2 .
i 1/i
(7.14)
For a fit to a straight line, the overall quality of the fit can be measured by the reduced
chi-squared parameter
2r =
2
n2
(7.15)
where 2 is given by Eq. (7.11) evaluated at the optimal values of a and b given by
Eq. (7.13). A good fit is characterized by 2r 1. This makes sense because if the
uncertainties i have been properly estimated, then [yi f (xi )]2 should on average be
roughly equal to i2 , so that the sum in Eq. (7.11) should consist of n terms approximately
equal to 1. Of course, if there were only 2 terms (n=2), then 2 would be zero as the best
straight line fit to two points is a perfect fit. That is essentially why 2r is normalized
using n 2 instead of n. If 2r is significantly greater than 1, this indicates a poor fit to
the fitting function (or an underestimation of the uncertainties i ). If 2r is significantly
less than 1, then it indicates that the uncertainties were probably overestimated (the fit and
fitting function may or may not be good).
We can also get estimates of the uncertainties in our determination of the fitting parameters a and b, although deriving the formulas is a bit more involved that we want to get into
here. Therefore, we just give the results:
1
) xi /i2
i (xi x
P 2 2
x /
a2 = b2 Pi i 2i .
i 1/i
b2 = P
(7.16)
The estimates of uncertainties in the fitting parameters depend explicitly on {i } and will
only be meaningful if (i) 2r 1 and (ii) the estimates of the uncertainties i are accurate.
You can find more information, including a derivation of Eq. (7.16), in Data Reduction
and Error Analysis for the Physical Sciences, 3rd ed by P. R. Bevington & D. K. Robinson, McGraw-Hill, New York, 2003.
135
Fit to y = ax + b
velocity (m/s)
5
10
15
20
25
0.0
0.5
1.0
1.5
time (s)
2.0
2.5
3.0
Figure 7.3: Fit using 2 least squares fitting routine with data weighted by error bars.
7.4 Exercises
1. Write a function that can return each of the first three spherical Bessel functions
jn (x):
sin x
x
sin x cos x
j1 (x) = 2
x
x
3
sin x 3 cos x
j2 (x) =
1
2
x
x
x2
j0 (x) =
(7.17)
Your function should take as arguments a NumPy array x and the order n, and
should return an array of the designated order n spherical Bessel function. Take
care to make sure that your functions behave properly at x = 0.
Demonstrate the use of your function by writing a Python routine that plots the
three Bessel functions for 0 x 20. Your plot should look like the one below.
Something to think about: You might note that j1 (x) can be written in terms of
j0 (x), and that j2 (x) can be written in terms of j1 (x) and j0 (x). Can you take
136
Chapter 7. Functions
advantage of this to write a more efficient function for the calculations of j1 (x) and
j2 (x)?
2.
(a) Write a function that simulates the rolling of n dice. Use the NumPy function
random.random_integers(6), which generates a random integer between 1 and 6 with equal probability (like rolling fair dice). The input of your
function should be the number of dice thrown each roll and the output should
be the sum of the n dice.
(b) Roll 2 dice 10,000 times keeping track of all the sums of each set of
rolls in a list. Then use your program to generate a histogram summarizing the rolls of two dice 10,000 times. The result should look like
the histogram plotted below. Use the MatPlotLib function hist (see
http://matplotlib.org/api/pyplot_summary.html) and set the number of bins in
the histogram equal to the number of different possible outcomes of a roll of
your dice. For example, the sum of two dice can be anything between 2 and
12, which corresponds to 11 possible outcomes. You should get a histogram
that looks like the one below.
(c) Repeat part (b) using 3 dice and plot the resulting histogram.
3. Write a function to draw a circular smiley face with eyes, a nose, and a mouth.
One argument should set the overall size of the face (the circle radius). Optional
arguments should allow the user to specify the (x, y) position of the face, whether
7.4. Exercises
137
sum of dice
the face is smiling or frowning, and the color of the lines. The default should
be a smiling blue face centered at (0, 0). Once you write your function, write
a program that calls it several times to produce a plot like the one below (creative improvisation is encouraged!). In producing your plot, you may find the call
plt.axes().set_aspect(1) useful so that circles appear as circles and not
ovals. You should only use MatPlotLib functions introduced in this text. To create
a circle you can create an array of angles that goes from 0 to 2 and then produce
the x and y arrays for your circle by taking the cosine and sine, respectively, of the
array. Hint: You can use the same (x, y) arrays to make the smile and frown as you
used to make the circle by plotting appropriate slices of those arrays. You do not
need to create new arrays.
4. In the section Example: linear least squares fitting, we showed that the best fit of a
line y = a + bx to a set of data {(xi , yi )} is obtained for the values of a and b given
by Eq. (7.10). Those formulas were obtained by finding the values of a and b that
minimized the sum in Eq. (7.5). This approach and these formulas are valid when
the uncertainties in the data are the same for all data points. The Python function
LineFit(x, y) in the section Example: linear least squares fitting implements
Eq. (7.10).
(a) Write a new fitting function LineFitWt(x, y) that implements the formulas given in Eq. (7.14) that minimize the 2 function give by Eq. (7.12).
This more general approach is valid when the individual data points have dif-
138
Chapter 7. Functions
ferent weightings or when they all have the same weighting. You should also
write a function to calculate the reduced chi-squared 2r defined by Eq. (7.12).
(b) Write a Python program that reads in the data below, plots it, and fits it using the two fitting functions LineFit(x, y) and LineFitWt(x, y).
Your program should plot the data with error bars and with both fits with
and without weighting, that is from LineFit(x, y) and LineFitWt(x,
y, dy). It should also report the results for both fits on the plot, similar to
the output of the supplied program above, as well as the values of 2r , the reduce chi-squared value, for both fits. Explain why weighting the data gives a
steeper or less steep slope than the fit without weighting.
Velocity vs time data
for a falling mass
time (s)
velocity (m/s)
2.23
139
4.78
123
7.21
115
9.37
96
11.64
62
14.23
54
16.55
10
18.70
-3
21.05
-13
7.4. Exercises
uncertainty (m/s)
16
16
4
9
17
17
12
15
18
139
23.21
-55
10
140
Chapter 7. Functions
CHAPTER
EIGHT
CURVE FITTING
One of the most important tasks in any experimental science is modeling data and determining how well some theoretical function describes experimental data. In the last
chapter, we illustrated how this can be done when the theoretical function is a simple
straight line in the context of learning about Python functions and methods. Here we
show how this can be done for a arbitrary fitting functions, including linear, exponential,
power law, and other nonlinear fitting functions.
(8.1)
141
This equation is nonlinear in t and in the fitting parameter and thus cannot be fit using
the method of the previous chapter. Fortunately, this is a special case for which the fitting
function can be transformed into a linear form. Doing so will allow us to use the fitting
routine we developed for fitting linear functions.
We begin our analysis by transforming our fitting function to a linear form. To this end
we take the logarithm of Eq. (8.1):
ln N = ln N0
t
.
With this tranformation our fitting function is linear in the independent variable t. To
make our method work, however, our fitting function must be linear in the fitting parameters, and our transformed function is still nonlinear in the fitting parameters and N0 .
Therefore, we define new fitting parameters as follows
a = ln N0
(8.2)
b = 1/
Now if we define a new dependent variable y = ln N , then our fitting function takes the
form of a fitting function that is linear in the fitting parameters a and b
y = a + bx
where the independent variable is x = t and the dependent variable is y = ln N .
We are almost ready to fit our transformed fitting function, with transformed fitting parameters a and b, to our transformed independent and dependent data, x and y. The last
thing we have to do is to transform the estimates of the uncertainties N in N to the uncertainties y in y (= ln N ). So how much does a given uncertainty in N translate into an
uncertainty in y? In most cases, the uncertainty in y is much smaller than y, i.e. y y;
similarly N N . In this limit we can use differentials to figure out the relationship
between these uncertainties. Here is how it works for this example:
y = ln N
y
N
y =
N
y =
(8.3)
N
.
N
Equation (8.3) tells us how a small change N in N produces a small change y in y. Here
we identify the differentials dy and dN with the uncertainties y and N . Therefore, an
uncertainty of N in N corresponds, or translates, to an uncertainty y in y.
142
Lets summarize what we have done so far. We started with the some data points {ti , Ni }
and some addition data {Ni } where each datum Ni corresponds to the uncertainty in
the experimentally measured Ni . We wish to fit these data to the fitting function
N (t) = N0 et/ .
We then take the natural logarithm of both sides and obtain the linear equation
ln N = ln N0
(8.4)
y = a + bx
with the obvious correspondences
x=t
y = ln N
a = ln N0
(8.5)
b = 1/
Now we can use the linear regression routine with 2 weighting that we developed in the
previous section to fit (8.4) to the transformed data xi (= ti ) and yi (= ln Ni ). The inputs
are the tranformed data xi , yi , yi . The outputs are the fitting parameters a and b, as well
as the estimates of their uncertainties a and b along with the value of 2 . You can obtain
the values of the original fitting parameters N0 and by taking the differentials of the last
two equations in Eq. (8.5):
a
N0 = N0
a =
N0
N0
(8.6)
b
b = = 2
The Python routine below shows how to implement all of this for a set of experimental
data that is read in from a data file.
Figure 8.1 shows the output of the fit to simulated beta decay data obtained using the
program below. Note that the error bars are large when the number of counts N are
small. This is consistent with what is known as shot noise (noise that arises from counting
discrete events), which obeys Poisson statistics. You will study sources of noise, including
shot noise, later in your lab courses. The program also prints out the fitting parameters of
the transformed data as well as the fitting parameters for the exponential fitting function.
import numpy as np
import matplotlib.pyplot as plt
143
7
6
ln(N)
5
4
3
2
0
20
40
60
80
100
144
145
(8.7)
We follow the same approach we used for the exponential fitting function and first take
the logarithm of both sides of (8.7)
ln P = ln P0 + ln s .
(8.8)
We recast this in the form of a linear equation y = a + bx with the following identifications:
x = ln s
y = ln P
a = ln P0
(8.9)
b=
Following a procedure similar to that used to fit using an exponential fitting function, you
can use the tranformations given by (8.9) as the basis for a program to fit a power-law
fitting function such as (8.8) to experimental data.
the linear optimization problem). We can gain some insight into this nonlinear optimization problem, namely the fitting of a nonlinear fitting function to a data set, by considering
a fitting function with only two fitting parameters. That is, we are trying to fit some data
set {xi , yi }, with uncertainties in {yi } of {i }, to a fitting function is f (x; a, b) where a
and b are the two fitting parameters. To do so, we look for the minimum in
2 (a, b) =
X yi f (xi ) 2
i
Note that once the data set, uncertainties, and fitting function are specified, 2 (a, b) is
simply a function of a and b. We can picture the function 2 (a, b) as a of landscape with
peaks and valleys: as we vary a and b, 2 (a, b) rises and falls. The basic idea of all
nonlinear fitting routines is to start with some initial guesses for the fitting parameters,
here a and b, and by scanning the 2 (a, b) landscape, find values of a and b that minimize
2 (a, b).
There are a number of different methods for trying to find the minimum in 2 for nonlinear fitting problems. Nevertheless, the method that is most widely used goes by the
name of the Levenberg-Marquardt method. Actually the Levenberg-Marquardt method
is a combination of two other methods, the steepest descent (or gradient) method and
parabolic extrapolation. Roughly speaking, when the values of a and b are not too near
their optimal values, the gradient descent method determines in which direction in (a, b)space the function 2 (a, b) decreases most quicklythe direction of steepest descent
and then changes a and b accordingly to move in that direction. This method is very
efficient unless a and b are very near their optimal values. Near the optimal values of
a and b, parabolic extrapolation is more efficient. Therefore, as a and b approach their
optimal values, the Levenberg-Marquardt method gradually changes to the parabolic extrapolation method, which approximates 2 (a, b) by a Taylor series second order in a and
b and then computes directly the analytical minimum of the Taylor series approximation
of 2 (a, b). This method is only good if the second order Taylor series provides a good
approximation of 2 (a, b). That is why parabolic extrapolation only works well very near
the minimum in 2 (a, b).
Before illustrating the Levenberg-Marquardt method, we make one important cautionary
remark: the Levenberg-Marquardt method can fail if the initial guesses of the fitting parameters are too far away from the desired solution. This problem becomes more serious
the greater the number of fitting parameters. Thus it is important to provide reasonable
initial guesses for the fitting parameters. Usually, this is not a problem, as it is clear from
the physical situation of a particular experiment what reasonable values of the fitting parameters are. But beware!
The scipy.optimize module provides routines that implement the LevenbergMarquardt non-linear fitting method. One is called scipy.optimize.leastsq.
8.2. Nonlinear fitting
147
148
120
100
80
60
40
20
00
15
10
frequency (THz)
20
25
import
import
import
import
numpy as np
matplotlib.pyplot as plt
matplotlib.gridspec as gridspec
scipy.optimize
5
6
7
8
9
10
11
12
13
149
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
150
59
60
61
62
63
64
65
66
67
68
69
70
71
72
plt.show()
The above code also plots the difference between the data and fit, known as the residuals
in the subplot below the plot of the data and fit. Plotting the residuals in this way gives
a graphical representation of the goodness of the fit. To the extent that the residuals vary
randomly about zero and do not show any overall upward or downward curvature, or any
long wavelength oscillations, the fit would seem to be a good fit.
Finally, we note that we have used the MatPlotLib package gridspec to create the two
subplots with different heights. The gridspec are made in lines 3 (where the package
is imported), 36 (where 2 rows and 1 column are specified with relative heights of 6 to 2),
39 (where the first gs[0] height is specified), and 54 (where the second gs[1] height
is specified). More details about the gridspec package can be found at the MatPlotLib
web site.
151
120
100
80
/2fw2
a = 57.1 2.8
b = -2.55 0.80
c = 0.02 0.03
P = 77.5 3.5
fp = 11.1 0.1
fw = 1.8 0.1
r2 = 2.03
60
40
20
00
15
10
frequency (THz)
20
25
15
10
frequency (THz)
20
25
residuals
20
0
200
152
8.3 Exercises
1. When a voltage source is connected across a resistor and inductor in series, the
voltage across the inductor Vi (t) is predicted to obey the equation
V (t) = V0 et
(8.10)
where t is the time and the decay rate = R/L is the ratio of the resistance R to
the inductance L of the circuit. In this problem, you are to write a Python routine
that fits the above equation to the data below for the voltage measured across an
inductor after it is connected in series with a resistor to a voltage source. Following
the example in the text, linearize the (8.10) and use a linear fitting routine, either
the one you wrote from the previous chapter or one from NumPy or SciPy.
(a) Find the best values of and V0 and the uncertainties in their values and
V0 .
(b) Find the value of 2r for your fit. Does it make sense?
(c) Make a semi-log plot of the data using symbols with error bars (no line) and
of the fit (line only). The fit should appear as a straight line that goes through
the data points.
(d) If the resistor has a value of 10.0 k, what is the value of the inductance and
its uncertainty according to your fit, assuming that the error in the resistance
is negligibly small.
Data for decay of voltage across an inductor
in an RL circuit
Date: 24-Oct-2012
Data taken by D. M. Blantogg and T. P. Chaitor
time (ns)
0.0
32.8
65.6
98.4
131.2
164.0
196.8
229.6
262.4
295.2
328.0
360.8
8.3. Exercises
voltage (volts)
5.08e+00
3.29e+00
2.23e+00
1.48e+00
1.11e+00
6.44e-01
4.76e-01
2.73e-01
1.88e-01
1.41e-01
9.42e-02
7.68e-02
uncertainty (volts)
1.12e-01
9.04e-02
7.43e-02
6.05e-02
5.25e-02
4.00e-02
3.43e-02
2.60e-02
2.16e-02
1.87e-02
1.53e-02
1.38e-02
153
393.6
426.4
459.2
492.0
3.22e-02
3.22e-02
1.98e-02
1.98e-02
8.94e-03
8.94e-03
7.01e-03
7.01e-03
2. Small nanoparticles of soot suspended in water start to aggregate when salt is added.
The average radius r of the aggregates is predicted to grow as a power law in time
t according to the equation r = r0 tn . Taking the logarithm of this equation gives
ln r = n ln t + ln r0 . Thus the data should fall on a straight line if ln r is plotted vs
ln t.
(a) Plot the data below on a graph of ln r vs ln t to see if the data fall approximately on a straight line.
Size of growing aggregate
Date: 19-Nov-2013
Data taken by M. D. Gryart and M. L. Waites
time (m)
size (nm)
unc (nm)
0.12
115
10
0.18
130
12
0.42
202
14
0.90
335
18
2.10
510
20
6.00
890
30
18.00
1700
40
42.00
2600
50
(b) Defining y = ln r and x = ln t, use the linear fitting routine you wrote for the
previous problem to fit the data and find the optimal values for the slope and
y intercept, as well as their uncertainties. Use these fitted values to find the
optimal values of the the amplitude r0 and the power n in the fitting function
r = r0 tn . What are the fitted values of r0 and n? What is the value of 2r ?
Does a power law provide an adequate model for the data?
3. In this problem you explore using a non-linear least square fitting routine to fit the
data shown in the figure below. The data, including the uncertainties in the y values,
are provided in the table below. Your task is to fit the function
2
/2 2
+C
(8.11)
154
45
40
35
30
25
20
150
10
15
20
time (ms)
25
30
35
40
8.3. Exercises
155
should be much larger than the one you found for part (c). The program has
found a local minimum in 2 one that is obviously is not the best fit!
(d) Setting the fitting parameters A, B, C, and to the optimal values you found
in part (b), plot 2r as a function of for spanning the range from 0.05 to
3.95. You should observe several local minima for different values of 2r ; the
global minimum in 2r should occur for the optimal value of you found in
part (b).
Data for absorption spectrum
Date: 21-Nov-2012
Data taken by P. Dubson and M. Sparks
time (ms) signal uncertainty
0.2
41.1
0.9
1.4
37.2
0.9
2.7
28.3
0.9
3.9
24.8
1.1
5.1
27.8
0.8
6.4
34.5
0.7
7.6
39.0
0.9
8.8
37.7
0.8
10.1
29.8
0.9
11.3
22.2
0.7
12.5
22.3
0.6
13.8
26.7
1.1
15.0
30.4
0.7
16.2
32.6
0.8
17.5
28.9
0.8
18.7
22.9
1.3
19.9
21.7
0.9
21.1
22.1
1.0
22.4
22.3
1.0
23.6
26.3
1.0
24.8
26.2
0.8
26.1
21.4
0.9
27.3
20.0
1.0
28.5
20.1
1.2
29.8
21.2
0.5
31.0
22.0
0.9
32.2
21.6
0.7
33.5
21.0
0.7
34.7
19.7
0.9
35.9
17.9
0.9
37.2
18.1
0.8
38.4
18.9
1.1
156
CHAPTER
NINE
NUMERICAL ROUTINES:
SCIPY AND NUMPY
SciPy is a Python library of mathematical routines. Many of the SciPy routines are Python
wrappers, that is, Python routines that provide a Python interface, for numerical libraries and routines originally written in Fortran, C, or C++. Thus, SciPy lets you take
advantage of the decades of work that has gone into creating and optimizing numerical
routines for science and engineering. Because the Fortran, C, or C++ code that Python
accesses is compiled, these routines typically run very fast. Therefore, there is no real
downsideno speed penaltyfor using Python in these cases.
We have already encountered one of SciPys routines, scipy.optimize.leastsq,
for fitting nonlinear functions to experimental data, which was introduced in the the chapter on Curve Fitting. Here we will provide a further introduction to a number of other
SciPy packages, in particular those on special functions, linear algebra, finding roots of
scalar functions, discrete Fourier transforms, and numerical integration, including routines for numerically solving ordinary differential equations (ODEs). Our introduction to
these capabilities does not include extensive background on the numerical methods employed; that is a topic for another text. Here we simply introduce the SciPy routines for
performing some of the more frequently required numerical tasks.
One final note: SciPy makes extensive use of NumPy arrays, so NumPy should always be
imported with SciPy
157
1.0
Bessel
0.5
0.0
0.5
1.00
1.0
0.8
0.6
0.4
0.2
0.00.0
15
10
20
Error
0.5
1.0
2.0
2.5
Legendre
0.5
0.0
0.5
1.01.0
0.5
0.0
Gamma
0.5
1.0
10
8
6
4
2
0
2
4
Airy
0.4
0.2
0.0
0.2
0.4
1.5
1.0
100
80
60
40
20
0
20
15
10
Laguerre
158
1
2
3
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# gamma function
x = np.linspace(-3.5, 6., 3601)
g = sp.special.gamma(x)
g = np.ma.masked_outside(g, -100, 400)
ax2 = fig.add_subplot(322)
ax2.plot(x,g)
ax2.set_xlim(-3.5, 6)
ax2.set_ylim(-20, 100)
ax2.text(0.5, 0.95,Gamma, ha=center, va=top,
transform = ax2.transAxes)
30
31
32
33
34
35
36
37
38
# error function
x = np.linspace(0, 2.5, 256)
ef = sp.special.erf(x)
ax3 = fig.add_subplot(323)
ax3.plot(x,ef)
ax3.set_ylim(0,1.1)
ax3.text(0.5, 0.95,Error, ha=center, va=top,
transform = ax3.transAxes)
39
40
41
42
43
44
# Airy function
x = np.linspace(-15, 4, 256)
ai, aip, bi, bip = sp.special.airy(x)
ax4 = fig.add_subplot(324)
ax4.plot(x,ai, x,bi)
159
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
# Legendre polynomials
x = np.linspace(-1, 1, 256)
lp0 = np.polyval(sp.special.legendre(0),x)
lp1 = np.polyval(sp.special.legendre(1),x)
lp2 = np.polyval(sp.special.legendre(2),x)
lp3 = np.polyval(sp.special.legendre(3),x)
ax5 = fig.add_subplot(325)
ax5.plot(x,lp0, x,lp1, x,lp2, x,lp3)
ax5.axhline(color="grey", ls="--", zorder=-1)
ax5.set_ylim(-1,1.1)
ax5.text(0.5, 0.9,Legendre, ha=center, va=top,
transform = ax5.transAxes)
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
# Laguerre polynomials
x = np.linspace(-5, 8, 256)
lg0 = np.polyval(sp.special.laguerre(0),x)
lg1 = np.polyval(sp.special.laguerre(1),x)
lg2 = np.polyval(sp.special.laguerre(2),x)
lg3 = np.polyval(sp.special.laguerre(3),x)
ax6 = fig.add_subplot(326)
ax6.plot(x,lg0, x,lg1, x,lg2, x,lg3)
ax6.axhline(color="grey", ls="--", zorder=-1)
ax6.axvline(color="grey", ls="--", zorder=-1)
ax6.set_xlim(-5,8)
ax6.set_ylim(-5,10)
ax6.text(0.5, 0.9,Laguerre, ha=center, va=top,
transform = ax6.transAxes)
78
79
plt.show()
The arguments of the different functions depend, of course, on the nature of the particular
function. For example, the first argument of the two types of Bessel functions called
in lines 10-13 is the so-called order of the Bessel function, and the second argument is
the independent variable. The Gamma and Error functions take one argument each and
produce one output. The Airy function takes only one input argument, but returns four
outputs, which correspond the two Airy functions, normally designated Ai(x) and Bi(x),
and their derivatives Ai0 (x) and Bi0 (x). The plot shows only Ai(x) and Bi(x).
160
The polynomial functions shown have a special syntax that uses NumPys polyval
function for generating polynomials. If p is a list or array of N numbers and x is an array,
then
polyval(p, x) = p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x +
p[N-1]
import scipy.linalg
a = array([[-2, 3], [4, 5]])
a
array([[-2, 3],
[ 4, 5]])
In [5]: scipy.linalg.det(a)
Out[5]: -22.0
The inverse of a matrix is computed using the scipy.linalg.inv function, while the
product of two matrices is calculated using the NumPy dot function:
9.2. Linear algebra
161
In [6]: b = scipy.linalg.inv(a)
In [6]: b
Out[6]: array([[-0.22727273,
[ 0.18181818,
In [7]: dot(a,b)
Out[7]: array([[ 1.,
[ 0.,
0.13636364],
0.09090909]])
0.],
1.]])
2 4
6
x1
4
A = 1 3 9 , x = x2 , b = 11 .
8 5 7
x3
1
Next we construct the array A and vector b as NumPy arrays:
In [8]: A = array([[2, 4, 6], [1, -3, -9], [8, 5, -7]])
In [9]: b = array([4, -11, 2])
-3.17391304])
162
10.2173913 ,
-3.17391304])
2 4
6
A = 1 3 9
1 2
3
where the third row is a multiple of the first row. Lets try it out and see what happens.
First we change the bottom row of the matrix A and then try to solve the system as we did
before.
In [11]: A[2] = array([1, 2, 3])
In [12]: A
Out[12]: array([[ 2, 4, 6],
[ 1, -3, -9],
[ 1, 2, 3]])
In [13]: scipy.linalg.solve(A,b)
LinAlgError: Singular matrix
In [14]: Ainv = scipy.linalg.inv(A)
LinAlgError: Singular matrix
163
eigenvectors. Lets return to the matrix we were using previously and find its eigenvalues
and eigenvectors.
A = array([[2, 4, 6],[1, -3, -9],[8, 5, -7]])
In [15]: A
Out[15]: array([[ 2, 4, 6],
[ 1, -3, -9],
[ 8, 5, -7]])
In [16]: lam, evec = scipy.linalg.eig(A)
In [17]: lam
Out[17]: array([ 2.40995356+0.j, -8.03416016+0.j,
-2.37579340+0.j])
In [18]: evec
Out[18]: array([[-0.77167559, -0.52633654, 0.57513303],
[ 0.50360249, 0.76565448, -0.80920669],
[-0.38846018, 0.36978786, 0.12002724]])
0.50360249, -0.38846018])
1.21365861, -0.93617101])
In [22]: lam[0]*evec[:,0]
Out[22]: array([-1.85970234+0.j, 1.21365861+0.j,
-0.93617101+0.j])
Thus we see by direct substitution that the left and right sides of Ax = x : are equal. In
general, the eigenvalues can be complex, so their values are reported as complex numbers.
164
Then we can solve the generalized eigenvalue problem by entering B as the optional
second argument to scipy.linalg.eig
In [23]: lam, evec = scipy.linalg.eig(A,B)
The solutions are returned in the same fashion as before, as an array lam whose entries
are the eigenvalues and a matrix evac whose rows are the eigenvectors.
In [24]: lam
Out[24]: array([-1.36087907+0.j, 0.83252442+0.j,
-0.10099858+0.j])
In [25]: evec
Out[25]: array([[-0.0419907 , -1.
, 0.93037493],
[-0.43028153, 0.17751302, -1.
],
[ 1.
, -0.29852465, 0.4226201 ]])
165
are relatively simple and appropriate for the most common types of nonlinear equations.
p
Plots of tan x and (8/x)2 1 vs x are shown in the top plot in the figure Crossing
functions, albeit with x replaced by
p . The solutions to this equation are those x values
where the two curves tan x and (8/x)2 1 cross each other. The first step towards
obtaining a numerical solution is to rewrite the equation to be solved in the form f (x) = 0.
Doing so, the above equation becomes
p
tan x (8/x)2 1 = 0
Obviously the two equations above have the same solutions for x. Parenthetically we
mention that the problem of finding the solutions to equations of the form f (x) = 0 is
often referred to as finding the roots of f (x).
Next, we plot f (x) over the domain of interest, in this case from x = 0 to 8. We are only
interested in positive solutions and for x > 8, the equation has no real solutions as the
argument of the square root becomes negative. The solutions, the points where f (x) = 0
are indicated by green circles; there are three of them. Another notable feature of the
function is that it diverges to at x = {0, /2, 3/2, 5/2}.
Brent method
One of the workhorses for finding solutions to a single variable nonlinear equation is the
method of Brent, discussed in many texts on numerical methods. SciPys implementation of the Brent algorithm is the function scipy.optimize.brentq(f, a, b),
which has three required arguments. The first f is the name of the user-defined function
166
tanx
true roots
false roots
x
Figure 9.2: Roots of a nonlinear function
to be solved. The next two, a and b are the x values that bracket the solution you are
looking for. You should choose a and b so that there is only one solutions in the interval between a and b. Brents method also requires that f(a) and f(b) have opposite
signs; anperror message is returned if they do not. Thus to find the three solutions to
tan x (8/x)2 1 = 0, we need to run scipy.optimize.brentq(f, a, b)
three times using three different values of a and b that bracket each of the three solutions.
The program below illustrates the how to use scipy.optimize.brentq
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
def tdl(x):
y = 8./x
return np.tan(x) - np.sqrt(y*y-1.0)
# Find true roots
rx1 = sp.optimize.brentq(tdl, 0.5, 0.49*np.pi)
rx2 = sp.optimize.brentq(tdl, 0.51*np.pi, 1.49*np.pi)
rx3 = sp.optimize.brentq(tdl, 1.51*np.pi, 2.49*np.pi)
167
168
True roots:
f(1.39547) = -6.39e-14
f(4.16483) = -7.95e-14
f(6.83067) = -1.11e-15
False roots:
f(1.57080) = -1.61e+12
f(4.71239) = -1.56e+12
f(7.85398) = 1.16e+12
The Brent method finds the three true roots of the equation quickly and accurately when
you provide values for the brackets a and b that are valid. However, like many numerical
methods for finding roots, the Brent method can produce spurious roots as it does in the
above example when a and b bracket singularities like those at x = {/2, 3/2, 5/2}.
Here we evaluated the function at the purported roots found by brentq to verify that
the values of x found were indeed roots. For the true roots, the values of the function
were very near zero, to within an acceptable roundoff error of less than 1013 . For the
false roots, exceedingly large numbers on the order of 1012 were obtained, indicating a
possible problem with these roots. These results, together with the plots, allow you to
unambiguously identify the true solutions to this nonlinear function.
The brentq function has a number of optional keyword arguments that you may find
useful. One keyword argument causes brentq to return not only the solution but the
value of the function evaluated at the solution. Other arguments allow you to specify a
tolerance to which the solution is found as well as a few other parameters possibly of
interest. Most of the time, you can leave the keyword arguments at their default values.
See the brentq entry online on the SciPy web site for more information.
Other methods for solving equations of a single variable
SciPy provides a number of other methods for solving nonlinear equations of a
single variable. It has an implementation of the Newton-Raphson method called
scipy.optimize.newton. Its the racecar of such methods; its super fast but less
stable that the Brent method. To fully realize its speed, you need to specify not only the
function to be solved, but also its first derivative, which is often more trouble than its
worth. You can also specify its second derivative, which may further speed up finding the
solution. If you do not specify the first or second derivatives, the method uses the secant
method, which is usually slower than the Brent method.
Other methods, including the Ridder (scipy.optimize.ridder) and bisection
(scipy.optimize.bisect), are also available, although the Brent method is generally superior. SciPy lets you use your favorite.
169
170
= f1 (t, y1 , ..., yn )
= f2 (t, y1 , ..., yn )
=
..
.
= fn (t, y1 , ..., yn ) .
We also need n initial conditions, one for each variable yi . Here we have a second order
ODE so we will have two coupled ODEs and two initial conditions.
We start by transforming our second order ODE into two coupled first order ODEs. The
transformation is easily accomplished by defining a new variable d/dt. With this
definition, we can rewrite our second order ODE as two coupled first order ODEs:
d
=
dt
d
1
= + sin + d cos t .
dt
Q
In this case the functions on the right hand side of the equations are
f1 (t, , ) =
f2 (t, , ) =
1
+ sin + d cos t .
Q
Note that there are no explicit derivatives on the right hand side of the functions fi ; they
are all functions of t and the various yi , in this case and .
The initial conditions specify the values of and at t = 0.
SciPys ODE solver scipy.integrate.odeint has three required arguments and
many optional keyword arguments, of which we only need one, args, for this example.
So in this case, odeint has the form
odeint(func, y0, t, args=())
The first argument func is the name of a Python function that returns a list of values of
the n functions fi (t, y1 , ..., yn ) at a given time t. The second argument y0 is an array (or
list) of the values of the initial conditions of y1 , ..., yn ). The third argument is the array
of times at which you want odeint to return the values of y1 , ..., yn ). The keyword
argument args is a tuple that is used to pass parameters (besides y0 and t) that are
needed to evaluate func. Our example should make all of this clear.
After having written the nth -order ODE as a system of n first-order ODEs, the next task
is to write the function func. The function func should have three arguments: (1) the
list (or array) of current y values, the current time t, and a list of any other parameters
params needed to evaluate func. The function func returns the values of the derivatives dyi /dt = fi (t, y1 , ..., yn ) in a list (or array). Lines 5-11 illustrate how to write
func for our example of a driven damped pendulum. Here we name the function simply
f, which is the name that appears in the call to odeint in line 33 below.
The only other tasks remaining are to define the parameters needed in the function,
bundling them into a list (see line 22 below), and to define the initial conditions, and
9.4. Solving ODEs
171
bundling them into another list (see line 25 below). After defining the time array in lines
28-30, the only remaining task is to call odeint with the appropriate arguments and a
variable, psoln in this case to store output. The output psoln is an n element array
where each element is itself an array corresponding the the values of yi for each time in
the time t array that was an argument of odeint. For this example, the first element
psoln[:,0] is the y0 or theta array, and the second element psoln[:,1] is the
y1 or omega array. The remainder of the code simply plots out the results in different
formats. The resulting plots are shown in the figure Pendulum trajectory after the code.
1
2
3
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
4
5
6
7
8
9
10
11
12
13
14
15
# Parameters
Q = 2.0
d = 1.5
Omega = 0.65
16
17
18
19
# Initial values
theta0 = 0.0
# initial angular displacement
omega0 = 0.0
# initial angular velocity
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# Plot results
fig = plt.figure(1, figsize=(8,8))
172
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
plt.tight_layout()
plt.show()
The plots above reveal that for the particular set of input parameters chosen Q = 2.0,
d = 1.5, and Omega = 0.65, the pendulum trajectories are chaotic. Weaker forcing
(smaller d) leads to what is perhaps the more familiar behavior of sinusoidal oscillations
with a fixed frequency which, at long times, is equal to the driving frequency.
(9.1)
173
theta
15
10
5
0
5
10
15
200
50
100
time
150
200
50
100
time
150
200
omega
2
1
0
1
2
30
3
omega
2
1
0
1
2
30
3
theta
174
where f is the Fourier transform variable; if t is time, then f is frequency. The inverse
transform is given by
Z
G(f ) ei 2f t df
g(t) =
(9.2)
Here we define the Fourier transform in terms of the frequency f rather than the angular
frequency = 2f .
The conventional Fourier transform is defined for continuous functions, or at least for
functions that are dense and thus have an infinite number of data points. When doing
numerical analysis, however, you work with discrete data sets, that is, data sets defined for
a finite number of points. The discrete Fourier transform (DFT) is defined for a function
gn consisting of a set of N discrete data points. Those N data points must be defined at
equally-spaced times tn = nt where t is the time between successive data points and
n runs from 0 to N 1. The discrete Fourier transform (DFT) of gn is defined as
Gl =
N
1
X
gn ei (2/N ) ln
(9.3)
n=0
where l runs from 0 to N 1. The inverse discrete Fourier transform (iDFT) is defined
as
N 1
1 X
gn =
Gl ei (2/N ) ln .
N
(9.4)
l=0
The DFT is usually implemented on computers using the well-known Fast Fourier Transform (FFT) algorithm, generally credited to Cooley and Tukey who developed it at AT&T
Bell Laboratories during the 1960s. But their algorithm is essentially one of many independent rediscoveries of the basic algorithm dating back to Gauss who described it as
early as 1805.
175
width = 2.0
freq = 0.5
t = np.linspace(-10, 10, 101)
# linearly space time array
g = np.exp(-np.abs(t)/width) * np.sin(2.0*np.pi*freq*t)
dt = t[1]-t[0]
G
f
f
G
=
=
=
=
fftpack.fft(g)
# FFT of g
fftpack.fftfreq(g.size, d=dt) # frequenies f[i] of g[i]
fftpack.fftshift(f)
# shift frequencies from min to max
fftpack.fftshift(G)
# shift G order to coorespond to f
The DFT has real and imaginary parts, both of which are plotted in the figure.
The fft function returns the N Fourier components of Gn starting with the zerofrequency component G0 and progressing to the maximum positive frequency component G(N/2)1 (or G(N 1)/2 if N is odd). From there, fft returns the maximum negative component GN/2 (or G(N 1)/2 if N is odd) and continues upward in frequency
until it reaches the minimum negative frequency component GN 1 . This is the standard
way that DFTs are ordered by most numerical DFT packages. The scipy.fftpack
function fftfreq creates the array of frequencies in this non-intuitive order such that
f[n] in the above routine is the correct frequency for the Fourier component G[n]. The
arguments of fftfreq are the size of the the orignal array g and the keyword argument d that is the spacing between the (equally spaced) elements of the time array (d=1
if left unspecified). The package scipy.fftpack provides the convenience function
fftshift that reorders the frequency array so that the zero-frequency occurs at the
middle of the array, that is, so the frequencies proceed monotonically from smallest (most
176
g(t)
0.8
0.6
0.4
0.2
0.0
0.2
0.4
0.6
0.810
10
0
t
real part
imaginary part
5
G(f)
10
0
5
10 3
0
f
177
9.6 Exercises
1. Use NumPys polyval function together with SciPy to plot the following functions:
(a) The first four Chebyshev polynomials of first kind. Plot these over the interval
from -1 to +1.
2
(b) The first four Hermite polynomials multiplied by ex /2 . Plot these on the
interval from -5 to +5. These are the first four wave functions of the quantum
mechanical simple harmonic oscillator.
178
APPENDIX
INSTALLING PYTHON
For scientific programming with Python, you need to install Python and three scientific
Python libraries: NumPy, SciPy, and MatPlotLib. There are many more libraries you can
install, but Python along with NumPy, SciPy, and MatPlotLib are those that are essential
for scientific programming.
Express, is completely free to everybody and contains all the libraries you will need for
this manual. The other, Canopy Basic, contains nearly every library you are ever likely
to need for scientific computing. It is free to academic users; others pay a fee. Go to
https://www.enthought.com/products/canopy/ and press the Get Canopy button, which
will take you to a page where you can either download Canopy Express or request and
academic license, which will allow you to download Canopy Basic.
Spyder provides a completely open source programming environment for
Python.
The entire Spyder distribution is free to all and can be found at
https://code.google.com/p/spyderlib/. It also includes nearly all the scientific libraries
you are likely to need for scientific computing.
In this manual, we assume you are using Canopy, but the Spyder interface is very similar
to Canopy so that most users should have no difficulty using Spyder with this manual. If
you choose to use Spyder, launch Spyder and then go to the Preferences menu and then
under the Console menu, select the Advanced settings tab; tick the box Start
an IPython kernel at startup (it may already be selected, in which case you
need to do nothing). You only need to do this once, which sets up the IPython console
when Spyder is launched. Once that is done, you should be able to follow everything
written in this manual.
10
11
180
12
13
Do not
14
15
16
your_first_name = Dana
your_last_name = Martin
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
plt.plot([0,1], r, [1,0], b)
plt.text( 0.5, 0.8, {0:s} {1:s}.format(your_first_name, your_last_name),
horizontalalignment=center,
size = x-large,
bbox=dict(facecolor=purple, alpha=0.4))
plt.text( 0.5, 0.1,
{1:s}\nscipy {2:s}\nnumpy {3:s}\nmatplotlib {4:s}\non {5:s}\n{6:s}.format(
classname,
term,
scipy.__version__,
numpy.__version__,
matplotlib.__version__,
platform.platform(),
socket.gethostname()
) ,
horizontalalignment=center
)
filename = your_last_name + _ + your_first_name + _ + term + .png
plt.title(*** E-mail the saved version of this plot, ***\n +
"{0:s}" to pine@nyu.edu.format(filename), fontsize=12)
plt.savefig(filename)
plt.show()
181
182
APPENDIX
IPYTHON NOTEBOOKS
IPython notebooks are useful for logging your work. Here we suggest that you use them
for logging and turning in homework assignments. You may also find them useful in
other contexts such as laboratory work. The IPython Notebook interface is similar to
Mathematica or Maple.
183
This will launch the IPython notebook web application and will display the IPython Notebook Dashboard as a page in your default web browser.
184
IPython notebooks in that directory will appear in a list on the IPython Notebook Dashboard. Clicking on one of them will launch that notebook.
185
186
187
allows the user to embed headings, explanatory notes, mathematics, and images.
comment. We can enter any text we wish, including mathematical expressions using the markup language Latex. (If you do not already know Latex, you can get
a brief introduction at these sites: http://en.wikibooks.org/wiki/LaTeX/Mathematics or
ftp://ftp.ams.org/pub/tex/doc/amsmath/short-math-guide.pdf.) Here we enter the following text:
The total distance $x$ traveled during a trip can be
obtained by integrating the velocity $v(t)$ over the
duration $T$ of the trip:
\begin{align}
x = \int_0^T v(t)\, dt
\end{align}
After entering the text, pressing Shift-Enter yields the result shown in Annotation
using a Markdown cell.
189
Prompt (PC), you could write the contents of any text file using the command cat
filename (Mac) or type filename (PC). You can execute the same operation from the
IPython prompt using the Unix (Mac) or DOS (PC) command preceded by an exclamation
point, as described in the section on System shell commands.
190
This will open the IPython Notebook Dashboard in your web browser, where you should
see a list of all the IPython notebooks in that directory (folder). Click on the name of the
notebook you want to open. It will appear in a new tab on your web browser as before.
Note that while all the input and output from the previous saved session is present, none
of it has been run. That means that none of the variables or other objects has been defined
in this new session. To initialize all the objects in the file, you must rerun the file. To
rerun the file, press the Cell menu and select Run All, which will re-execute all the
cells. You will have to re-enter any screen input requested by the notebook scripts. Now
you are ready to pick up where you left off the last time.
191
192
APPENDIX
PYTHON RESOURCES
This text provides an introduction to Python for science and engineering applications but
is hardly exhaustive. There are many other resources that you will want to tap. Here I
point out several that you may find useful.
psy-pi for SciPy, like everyone else. Who says I have to be consistent? (see Emerson)
http://matplotlib.org/api/pyplot_summary.html The Plotting Commands
Summary page for MatPlotLib. It has a search feature and links to all the
MatPlotLib documentation, which I use a lot. You can go the the main
MatPlotLib page, http://matplotlib.org/, but frankly, its less useful. The
site http://www.loria.fr/~rougier/teaching/matplotlib/ is also useful for
learning some MatPlotLib tricks.
http://ipython.org/ I go to this page mostly to learn about IPython Notebook (http://ipython.org/notebook.html) but its also useful if you need
information about the IPython interpreter, especially if you want to find
out more about IPython magic commands.
http://www.enthought.com/ I get my latest version of Python and all the
packages I need for scientific computing here. One stop shopping and
everything is free for academic users. They offer three distributions:
Express, Basic, and Professional. Express is free to all users and contains all the Python libraries, NumPy, MatPlotlib, SciPy, etc, that are
described in this manual. Basic includes all the packages Enthought
supports, which is likely to be everything you will ever need. Professional adds support services. Basic is free to academic users. One nice
feature of Canopy is its package manager, which makes it childs play
to update or add Python packages. This is a very nice feature, especially
for beginners. Canopy displaces the older Enthought EPD packages.
https://code.google.com/p/spyderlib/ Get the latest version of Spyder, the
alternative IDE to Enthoughts Canopy. Sypder is completely open
source and has a number of nice features, like introspection, not currently available in Canopy. Its package manager isnt as nice as
Canopys, but its pretty good. Spyder is very popular and can be installed easily on all platforms.
http://www.scipy.org/Mailing_Lists Go here if you want to sign up for a
mailing list for NumPy or SciPy, or if you want to report a bug. Mailing
lists give you access to a community of developers and users that can
often provide expert help. Just remember to be polite and respectful of
those helping you and also to those posting questions.
https://lists.sourceforge.net/lists/listinfo/matplotlib-users The
list for MatPlotLib. See paragraph immediately above.
194
mailing
C.2 Books
There are a lot of books on Python and there is no way I can provide reviews for all of
them. I have found that the book by Mark Lutz, Learning Python, published by OReilly
Media does the trick for most people. It doesnt have anything special for scientific programming, and thus does not cover the NumPy, SciPy, or MatPlotLib packages, but for
just about everything else, its an excellent resource. It gives a good introduction to object
oriented programming, or OOP, which I say little about in this text. The 3rd edition of
the book covers Python 2 while the 4th and 5th (current) editions cover Python 3. You
are probably better off getting the latest edition as everybody will soon be using Python
3. If you are using Python 2, as we do in this text, you can easily enough figure out the
differences between Python 3 and 2.
C.2. Books
195
196
INDEX
array (NumPy), 29
printing, 61
assignment operator, 17
Canopy
Code Editor, 20
installing Python, 179
IPython pane, 6
tab completion, 20
window, 5
conditionals, 99
applied to arrays, 119
curve fitting, 140
linear, 130
exponential function, 141
power law function, 146
power-law function, 146
with weighting, 133
nonlinear, 146
dictionary, 29, 46
functions
arguments
**kwargs, 122
*args, 122
keyword, 121
positional, 121
variable number, 122
fast array processing, 119
MatPlotLib
module, 14
MatPlotLib functions
ayhline, axhline, 79
figure, 79
legend, 79
plot, 79
savefig, 79
show, 79
tight_layout, 90
xlabel, ylabel, 79
module
importing, 24
MatPlotLib, 14
NumPy, 14
SciPy, 14
NumPy
array, 35
functions, 16
module, 14
output, 55
screen, 57
writing data to a file, 66
reserved words, 19
roots of equations, see solving non-linear
equations
SciPy
module, 14
scripts, 20
solving non-linear equations, 165
Bisection method, 169
Brent method, 166
Newton-Raphson method, 169
Ridder method, 169
systems of nonlinear equations, 169
Spyder
installing Python, 179
string, 29
strings, 30
tab completion
Canopy, 20
IPython, 9
variable, 17
variable names, 19
plots, 73
basic, 75
error bars, 81
interactive, 74
line and symbol specifiers, 80
log-log, 91
logarithmic axes, 89
masked arrays, 85
semi-log, 89
setting axis limits, 83
subplots, 87
programs, 20
Python
installing, 179
module, 14
random numbers, 48
198
Index