Project 4

Problem 1

Write a function that takes a unit vector $n$ and another vector $w$ and returns the reflection of $w$ across the plane normal to $n$. (For help normalizing a vector, try importing norm from numpy.linalg and dividing through by the norm.)

Bonus: Using the lambda notation we saw yesterday, can you write a function which takes a unit vector $n$ and returns a reflection function instead?

Problem 2

Consider a very simple weather model which says that on a sunny day, there's a 90% chance the next day will be sunny and a 10% chance the next day will be rainy. On rainy days, it says there's a 50% chance the next day will be sunny and a 50% chance it will be rainy. Let's encode this as a simple Markov chain transition matrix:

$$T = \begin{bmatrix} .9 & .5 \\ .1 & .5 \end{bmatrix}$$

In order to do long term weather prediction, we'll iteratively apply $T$ to some initial data $x_0$.

Generate a few choices of random initial data and compute the long term probability distribution. Does it depend on the initial state?

Careful to normalize the initial data x0 so the sum of the probabilties is 1! You can do this by dividing through by np.sum(x0)!

Problem 2

In this problem, we'll write a function to construct the "standard" $2n \times 2n$ symplectic matrix of the form:

$$\Omega = \begin{bmatrix} 0 & I_n \\ -I_n & 0 \end{bmatrix}$$

In order to do this, try putting together some indentity and zeros matrices as a block matrix using Numpy's bmat function.

Problem 3

Consider the energy function describing a simplified pendulum given by

$$E(x, v) = \frac{v^2}{2} - \cos(x)$$

As a consequence of conservation of energy, solutions to the system are constrained to level sets of $E(x, v)$.

Compute $E$ on the domain $[-4 \pi, 4 \pi] \times [-2 \pi, 2 \pi]$ then create a contour or filled contour plot using either plt.contour(E) or plt.contourf(E) respectively.

It's you preference if you prefer the filled or unfilled version; see which one you like more! If you want to play around with the availible colormaps to tweak the look, you can find them here!

Can you tell what some of the different contour regions correspond to physically?

Problem 4

Numpy provides a few convinience functions for doing file processing. Recall that in the last project we wrote out the trajectory of a harmonic oscilator in the format:

t0 x0 v0
t1 x1 v1
t2 x2 v2
...

Although we wrote our own method to read the data and convert it the floats, since the data is very uniform, we can easily leverage Numpy's np.loadtxt('some-file-name.txt') function to read in the data as an $n \times 3$ array whose columns are the t, x and v entries.

Give this a try and print the resulting array to verify that it matches what you have in your file. Notice that it will do both the reading and type conversion for you!

Now, let x be the slice of the $x$-column and v be the slice of the $v$-column. Do the $(x, v)$-pairs roughly lie on a circle centered at the origin?

This is useful when applicable, but not all data sets amenable to this; they must have a certain amount of regularity. In cases where you're mostly working with arrays, you can pair this up with Numpy's np.savetxt('output-file-name.txt', A) to save and load data to file in an easy way.

Problem 5: Challenge

The standard $n$-dimension simplex is the simplex bounding the region around the origin $(0, 0, \ldots, 0)$ and each of the "unit" points $(1, 0, \ldots, 0), (0, 1, \ldots, 0), \cdots, (0, 0, \ldots, 1)$.

It's known that the standard $n$-simplex has volume $1/n!$. What we'd like to do is write a function taking a generic simplex $\Delta$ given by $n+1$ points $(p_0, p_1, \ldots, p_n)$ and compute it's volume.

As a starting point, we can relate the two simplices via the affine map $A$ from the standard $n$-simplex defined by

$$Ax = \begin{bmatrix} p_1 - p_0 & p_2 - p_0 & \cdots & p_n - p_0 \end{bmatrix} x + p_0$$

Use the affine map $A$ and the volume of the standard simplex to compute the volume of $\Delta$.

You may find the det function in numpy.linalg useful for this!

If the description of $A$ is tricky, try drawing a picture in 2d and 3d of what the column vectors look like.

Problem 6: Challenge!

We'll try to implement Conway's Game of Life in a Numpy. The game starts out with a board of cells which are either alive or dead. The game proceeds according to the following update rules:

1. Any live cell with fewer than two live neighbours dies, as if caused by under-population.
2. Any live cell with two or three live neighbours lives on to the next generation.
3. Any live cell with more than three live neighbours dies, as if by overcrowding.
4. Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.

See if you can implement this or some variation using Numpy.

The "Challenge!" as opposed to just "Challenge" means I haven't thought too hard about a solution, so this may be more involved than I anticipated. Try it if you find it interesting. :-)