# Calculating solutions with python

Friday, 11:15-12:05, 2016-01-22, Osmond 216

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 Canopy, open the editor, 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
• Multiplaction.

>>> 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
• Substraction is treated like addition of negative numbers.

>>> 1 - 2 - 3
• Exponential powers are represented with double-asterix's.

>>> 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 is a little un-conventional in python version 2, so always be careful with it.

>>> 1 / 0.2
>>> 10 / 3
• Python also 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 also can do arbitrary precision integer arithmetic!

>>> 2**200 - 1
>>> 4**4**4
• Pyton also 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

## 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 interpretter's global namespace and shell scripts with the "import" keyword.

>>> import math

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

>>> help(math)

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

>>> math.exp(0.)
>>> math.pi

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

>>> 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 - interpretter 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
• urllib2 - for http scripting
• turtle - turtle graphics, ala 1970's logo

## 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 "Run" menu and select "Restart Kernel".

Now, there is standard stack of modules we use for scientific computing called the "scipy stack". Parts of this stack of modules are automatically loaded into Canopy's global namespace. These are numpy (which creates fast arrays), scipy (which supplies all sorts of functions and algorithms) for calculations, 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 the standard python interprettor, you can load these modules into the global namespace with the following three lines of 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 we taught you in Math 220, 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 generalized 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)
>>> map(f, range(0,10))

In the code above, range(0,10) returns a list of integers [0,1,2,...9]. 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(). You can use the function g = lambda x : x**3 - 10*x**2 + 31*x - 30 to TEST if your solutions are correct.

## Plotting

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

>>> 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.

>>> 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. `