# Calculating with python

( previous, home, next )

Goals of this laboratory…

• Learn how to do basic math in python
• Introduce you to the concept of a python module
• Show you the numpy, scipy, and matplotlib modules
• Show you how to perform simple calculations and plotting

Start up Spyder and work through examples below to see what you can learn.

## Basic math

In Laboratory 1, you saw a little about how to do math calculations in python. But let’s go through the basics to make sure we’re all on the same page.

  >>> 1 + 2
• Multiplication.

  >>> 2 * 3
• In expressions involving both addition and multiplication, multiplication takes precedence.

  >>> 1 + 2 * 3

This could return 9 or 7, depending on whether you do the addition or multiplication first, but python also always does the multiplication. Python has a precedence order for all operations. Use parentheses to override the default precedence.

  >>> (1 + 2) * 3
• Subtraction is treated like addition of negative numbers.

  >>> 1 - 2 - 3
• Exponential powers are represented with double-asterixes.

  >>> 2**3
>>> 3**3
>>> 9**0.5 + 1
• The percent sign is used to find remainders of a division

  >>> 17 % 5
>>> 7.5 % 2
• Division will default to floating point numbers when needed.

  >>> 1 / 0.2
>>> 10 / 3
• Python has complex numbers built in. Imaginary numbers are post-fixed with a “j” to indicate imaginariness. So, the square root of -1 is represented by 1j or -1j.

  >>> 1j**2
>>> (1-1j) * (1+1j)
>>> 1j**4 - 1
>>> z = 1 + 2j
>>> z.real
>>> z.imag
• Python can do arbitrary precision integer arithmetic!

  >>> 2**200 - 1
>>> 4**4**4
• Python has several relational operators that we can use to test relationships between numbers and variables.

  >>>  2 < 3
>>>  2 > 3
>>>  2 == 3

Note that you use double-equals to test equality in python. Single equals is used only for assignment, and will raise an error if used incorrectly.

  >>> 2 = 3

## Floating point arithmetic

For most math problems, it’s easier to work with decimals rather than integers. Computers do this using something we call floating point arithmetic. It works pretty well, but you should just be aware of some of it’s quarks because you’ll probably encounter them.

Inside a computer CPU, floating point numbers can be stored in 32 bits for single precision or 64 bits for double precision. Mostly, we (and python) use double precision these days, which gives about 16 digits of accuracy for each number. That’s usually good enough. For each line below, try to guess what the output will be before you run it.

[Show code]
>>> u = 0.0
>>> v = 0.000000000000000005781
>>> u == v
>>> u - v

We expect u and v to be different from each other, and python agrees. But not always. Numbers that differ only by a very small relative amounts are considered equal.

[Show code]
>>> x = 0.1
>>> y = 0.100000000000000005781
>>> x == y
>>> x - y

The differences between x and y are so small relative to their sizes that python considers them equal. Another special thing about computer math is that our computers use base-2 (binary) representations instead of base-10 decimal representations. This difference in representation combined with the finite precision of calculations leads to small errors in even simple cases. This is called “representation error”.

[Show code]
>>> 0.1 + 0.7
>>> 0.1 + 0.7 - 0.8
>>> 1.1 ** 2 - 1.21

Notice the result is PRACTICALLY zero, in scientific notation. (see the ‘e’ at the end. That’s used to represent the $$10^{-16}$$ part of the scientific notation, since primitive computers could only use ASCII characters).

Finite-precision arithmetic is one of the important topics covered in Math 451 and our other numerics courses, so we won’t spend much time on it unless essential. But it’s something you’ll probably want to study at some point. In particular, formulas that are the same formally are often not the same in practice.

[Show code]
>>> e = 1e-15
>>> x = 1 + e
>>> y = 1 - e
>>> u, v = (x/7 - y/7), (x-y)/7
>>> u
>>> v

Observe that u and v SHOULD be equal to each other, but in fact only agree to a few decimal places. Which one is closer to the true answer?

## Modules

The python language itself (like all good modern languages) is simple, using only about 30 words. Python’s power comes from the rich collection of well-documented “modules” that are available to perform complex tasks. A module is python script that contains functions (and classes) to perform related tasks, like a library in C. Perhaps the simplest example of this is the “math” module.

You can load a module into your python interpreter’s global namespace and shell scripts with the “import” keyword.

[Show code]
>>> import math

You can use help to see the contents of the module you have imported.

[Show code]
>>> help(math)

To use one of the functions from the module, you put the modules name as a prefix, like

[Show code]
>>> math.exp(0.)
>>> math.pi

Function names can also be imported directly into the working namespace, although this can be dangerous and is discouraged.

[Show code]
>>> from math import sin, pi, exp, log
>>> sin(0), sin(pi/2)
>>> exp(0), exp(1), log(exp(1))

All other modules work the same way, and provide a variety of capabilities. Here are some common modules.

• sys - interpreter and shell functions
• os - operating system functions like path, file testing,
• itertools - tools for fancy for-loops
• string - basic algorithms for manipulating strings
• pickle - tool for saving variables to a file or reloading them
• re - for regular expressions
• urllib - for http scripting

## Scientific computing modules

For this next part of the lab, it will be useful to start with a clean environment, so go up to the menus and select “Restart Kernel”.

While python’s math module is useful for common highschool calculations, it won’t work well in our course because of a variety of limitations. Instead, there is a standard stack of modules we use for scientific computing called the “scipy stack”. Parts of this stack of modules are automatically loaded into Spyder’s global namespace. These are numpy (which creates arrays you can use to represent vectors, matrices, and tensors), scipy (which supplies all sorts of functions and algorithms that work with arrays), and matplotlib (for making plots and pictures). If you know matlab already, there are many similar things to python, but also some differences. If you write a script or use some other python interpreter, you can load these modules into the global namespace with the following three lines of code.

[Show code]
>>> from numpy import *
>>> from scipy import *
>>> from matplotlib.pyplot import *

Now, let’s see how we can use the scipy stack to solve some problems you have seen before.

• Suppose we have the linear system $$x - 2 y = -3, \; y - 2 z = 0, \; -y + z = -1$$ which we want to solve for $$x,y,z$$. In matrix form $$Ax=b$$, $\begin{gather*} \begin{bmatrix} 1 & -2 & 0 \\ 0 & 1 & 2 \\ 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} -3 \\ 0 \\ -3 \end{bmatrix} \end{gather*}$ To solve this, we first make arrays containing the matrix and vector we know.

>>> A = array([[1,-2,0],[0,1,2],[0,-1,1]])
>>> b = array([[-3],[0],[-3]])
>>> print(A)
>>> print(b)

Now, the best methods for solving linear systems turn out to be much fancier than what you learned in your first matrix algebra class, but they are programmed into the linalg submodule of scipy. As you type the command below, note that when you type “(”, a help-dialog automatically pops up reminding you how linalg.solve works so you can enter things correctly – in this case, the matrix is the first argument and the vector is the second argument.

  >>> x = linalg.solve(A, b)

To check our answer, we can calculate $$A x - b$$, and make sure it’s entries are all zero or close to zero. But to do matrix multiplication of arrays like they are matrices, we have to use the .dot() operator insteady of standard multiplication.

  >>> A.dot(x) - b

If instead, you did A*x, this would return the Hadamard product of the matrix and vector.

  >>> A*x

linalg has several useful linear algebra algorithms, including eigvals() for calculating the eigenvalues of a matrix.

  >>> linalg.eigvals(A)
• In python, variables can point to generic objects, be they integers, floats, strings, list, or even functions. This last case is particular useful as it allows us to create and pass functions around to other functions. “def” is the usual way to create functions (see above), but for short functions, there’s a sweet trick called lambda-forms, in honor of Alonzo Church’s Lambda calculus. (this is a badly-choosen keyword, unfortunately, since it clashes with standard math notation frequently, and one of the few flaws in the python language)

>>> f = lambda x : x**2
>>> f(1)
>>> f(2)
>>> list(map(f, range(0,10)))

In the code above, range(0,10) returns an iterator over the integers [0,1,2,...9], map applies f to each element value returned by the iterator, and list executes the iterator and makes the complete list. In python, intervals are usually semi-closed, containing their lower-bound but not their upper-bound. The map(f, ...) function is a functional programming trick meaning “apply the function f to each of the things in the subsequent sequence”.

• Using these lambda functions, we can do things like integration using standard algorithms. In calculus, we learned $\int_{0}^{1} x^2 dx = \left. \frac{x^3}{3} \right|_0^1 = \frac{1}{3}$ Scipy has a module called integration for approximating integrals.

>>> from scipy import integrate
>>> integrate.quad(lambda x : x**2, 0, 1)

The answer returns a tuple with two values. The first is the estimate of the integral, which you’ll see is almost correct. The second is the estimated error of the quadrature (hence the name, quad). Feel free to experiment. There are actually several algorithms listed, though some important ones are also missing.

• There’s lots more in here. For example, see if you can find all three roots of the cubic polynomial equation $$x^3 - 10 x^2 + 31 x - 30 = 0$$ using the function roots(). This is a good chance to practice using the help function mentioned above. You can use the function g = lambda x : x**3 - 10*x**2 + 31*x - 30 to TEST if your solutions are correct.

## Plotting

Before running the following animations in spyder, follow the menus and go to “Tools -> Preferences -> IPython console -> Graphics”, and change “Backend” to “Automatic”. This should set up spyder so the figures below come out correctly.

The scipy stack has a useful set of plotting tools. For a simple example, here is a use of the rational parametric form of a circle.

[Show code]
>>> t = linspace(-10,10,64)
>>> x = (t*t-1)/(t*t+1)
>>> y = 2*t/(t*t+1)
>>> plot(x, y, 'ro-', cos(t), sin(t), 'k:')

The linspace(-10,10,64) function creates an array of 64 evenly spaced points from -10 to 10, including the endpoints. The extra strings in the plot command specify how to draw a line. Can you guess what each character means?

To clear the figure, you can use the command clf(). Of course, plots without labels are often more confusing than helpful. We can add titles and labels with the functions title(), ylabel(), xlabel(), and text.

>>> xlabel('x-values')
>>> title('A circle', fontsize=25)

To save your figure to a .png image file with the “savefig” command.

>>> cd Desktop
>>> savefig('mycircle.png')

Of course, vector formats for images are better, so we usually use portable document format (.pdf), encapsulated postscript (.eps), or scalable vector graphic format (.svg).

For a second example, let’s see if we can draw a cycloid.

[Show code]
>>> figure(2)
>>> clf()
>>> r = 1.
>>> theta =  linspace(0, 8 * pi, 257)
>>> x = r * (t - sin(t))
>>> y = r * (1 - cos(t))
>>> plot(x, y, 'r-')

The figure function is used to create a new figure or select an old figure for the plot. `