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

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.

Addition.

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

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

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.

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

For more information, see Matplotlib's gallery

Matplotlib will also allow us to construct animations. Check out this animation of a double pendulum and run it as a script in Canopy. How does the motion change as you change the mass of the second pendulum weight? Can you construct your own animation of the cycloid like this one?