Mathematica Project , Math21b  Spring 2016

Support

Welcome to the Mathematica computer project! In case problems with this assignment, please email Oliver at knill@math.harvard.edu.  During work, save the notebook frequently. The actual assignments are in the last blue box at the end of this notebook. This lab is due Friday, April 29. 
Submit electronically by sending the completed notebook. All problems are listed at the end of this document.  All problems have sample code in this notebook somewhere. Please only submit the code with answers. 
The total work should be minimal. It can be done in a 2-3 hours if you are in a hurry. Of course we hope you get interested in one or the other topic and hang around for more.  If you think you need too much time, network around or ask. Don’t hesitate to ask.

Cells and Expressions

The notebook is subdivided into cells. In this notebook, code cells are yellow. The text you right right now is a "text" cell, the next  gren cell is an "input" cell. 
To evalute an expression, click anywhere on the expresson, hold down the <Shift> key and press <Enter>. Here we compute the dot product of two vectors

{3,4,5,6}.{2,3,4,5}

Or the product of a matrix with a vector

{{3,4,5,6},{1,1,1,1}}.{2,1,1,1}

The dot product is the same than  the matrix product of v with transpose w. Matrices are arrays,  lists of lists:

{{3,4,5,6}}.Transpose[{{2,3,4,5}}]

The next example shows how we compute the inverse of a matrix:

Inverse[{{3,2,5},{3,2,1},{1,1,1}}]

In this notebook, input cells are  boxed in yellow for better readability.

Basic Computations

The following examples should be selfexplanatory. We start with multiplying a matrix with a vector.
If you want to see a result, say of a matrix operation, like it would appear in a text book, add // TraditionalForm

{{1,2,3},{3,4,5},{4,5,6}}.Inverse[{{2,3,4},{1,1,1},{1,3,6}}]

A={{2,3},{1,2}}; b={3,4}; Inverse[A].b

Solve[{x+y==1, x-y==3},{x,y}]

Det[{{2,3,5,6},{1,2,3,4},{1,1,1,1},{12,3,2,1}}]

Tr[{{a,b},{c,d}}]

Transpose[{{1,2,3},{3,4,5},{1,2,1}}] //MatrixForm

A={{1,2},{3,4}}; QR=QRDecomposition[A];Q=QR[[1]]; R=QR[[2]]; {TraditionalForm[Q],TraditionalForm[R]}

Here is the multiplication table matrix. We add “MatrixForm” to see it formatted as a matrix

Table[k *l,{k,10},{l,10}] //MatrixForm

Eigenvalues

With Eigenvalues[A], you compute the eigenvalues, with Eigensystem[A], both the eigenvalues and the eigenvectors. What happens for a random matrix? Lets make an experiment and plot the complex eigenvalues of a random 50x50 matrix. You might want to try with larger matrices to see a pattern.

A=Table[Random[]-1/2,{250},{250}]; 
eigenvalues=Eigenvalues[A]; 
point[z_]:=Point[{Re[z],Im[z]}]; Graphics[Map[point,eigenvalues]]

Sergey Sadov of the  Memorial university in Canada asked: 
Suppose A is an NxN matrix with entries independently and uniformly distributed on [0, 1].
What is the distribution of the absolute value of the largest eigenvalues of A for large N?

le:=Module[{},A=Table[Random[],{100},{100}]; Max[Abs[Eigenvalues[A]]]];

Lets do an experiment, take 1000 matrices, compute the largest eigenvalue in each case and look at the distribution. This computation takes some time.

s=Sort[Table[le,{1000}]];

If we plot these data, we see the distribution function:

ListPlot[s]

The density function can be visualized as follows:

ListPlot[BinCounts[s,{Min[s],Max[s],0.1}],Ticks->None]

While we are at it, lets also compute some statistical properties:

{Mean[s],Median[s],StandardDeviation[s],Min[s],Max[s]}

Here is how we compute the eigenvalues of a random graph:

Graph spectra

s=RandomGraph[{20,100}]
N[Eigenvalues[KirchhoffMatrix[s]]]

Mathematica can deal with graphs. Here are classes of graphs

s=WheelGraph[10]
N[Eigenvalues[KirchhoffMatrix[s]]]

Distribution of Ground state of a Graph

A connected graph always has the eigenvalue 0 as the eigenvector [1,1,1,1...,1] shows. The second eigenvalue (usually nozero) is interesting for physisists.The Next procedure computes the 2nd eigenvalue of a random graph with n vertices where an edge is turned on with probability p.

SecondEigenvalueRandomGraph[n_,p_]:=Sort[N[Eigenvalues[1.0*KirchhoffMatrix[RandomGraph[{n,Floor[p n(n-1)/2]}]]]]][[2]];

And here we see the distribution of that eigenvalue for a 1000 random graphs with 20 vertices and edge probability 1/2:

BarChart[BinCounts[Table[SecondEigenvalueRandomGraph[20,0.5],{1000}],0.1]]

LargestEigenvalueRandomGraph[n_,p_]:=Max[N[Eigenvalues[1.0*KirchhoffMatrix[RandomGraph[{n,Floor[p n(n-1)/2]}]]]]];

BarChart[BinCounts[Table[LargestEigenvalueRandomGraph[200,0.5],{10000}],0.1]]

Determinants

What is the distribution of the determinant of a random matrix? To see the distribution of a random 10x10 matrix, lets do 100 experiments and plot the results in a sorted manner.

S=Sort[Table[Det[RandomReal[{0,1},{10,10}]],{100}]];
ListPlot[S]

This distribution approaches the arcsin-distribution.

Plot[ArcSin[x],{x,-1,1}]

Here is how we compute the determinant of the 100 x 100 matrix which contains the entry  i*j at the entry (i,j).

A=Table[i*j,{i,100},{j,100}];
Det[A]

Differential Equations

Differential equations can be solved with the command "DSolve". Here is an example, where one does not specify the initial conditions:

ClearAll[f];DSolve[f'''[x]-2 f''[x] + f'[x]==Exp[3x],f[x],x]

Here is an example, where the initial conditions are specified. This is the solution of the forced driven harmonic oscillator:

S=DSolve[{f''[x]+f'[x]+f[x]==Cos[x],
f[0]==10,f'[0]==0},f[x],x]

You can then plot the solution

g[t_]:=S[[1,1,2]] /. x->t; Plot[g[t],{t,0,20}]

Sometimes, it is good to clear all variables. We have used f,g before and the next solution of a differential equation would make problems if

Remove["Global`*"]

S=DSolve[{f''[x]== -f[x] +g[x], g''[x]==f[x] -g[x],
f[0]==10,f'[0]==0,g[0]==1,g'[0]==0},{f[x],g[x]},x]

Markov Processes

A vector is called a probability vector if each entry is nonnegative and the sum of all its entry is 1. A matrix A for which each column is a probability vector is called a Markov matrix or stochastic matrix.

A={{1/2,1/5,1/4,1/3},{1/6,1/5,1/4,0},{1/6,1/5,1/4,1/3},{1/6,2/5,1/4,1/3}};
MatrixForm[A]

A  Markov matrix has an eigenvalue 1 because its transpose has the eigenvector [1,1,..., 1]. The eigenvector of the eigenvalue 1 of A  is called the steady state.

N[Eigensystem[A]]

When iterating the matrix A, it converges to a matrix which has the eigenvector as column vectors.

N[MatrixPower[A,10]]

Here is a magical square:

A={{17,24,1,8,15},{23,5,7,14,16},{4,6,13,20,22},{10,12,19,21,3},{11,18,25,2,9}};
MatrixForm[A]

If we divide it by the sum 65 over each column, we obtain a stochastic matrix.

Data fitting

The following example produces random data and fits them with a polynomial of degree 3.

h[x_]:=x+2*Random[];M=400;
data=Table[{k/M,h[k/M]},{k,M}];
f=Function[y,Fit[data,Table[z^j,{j,0,3}],z]/.z->y];
S1=ListPlot[data];
S2=Plot[f[x],{x,0,1},PlotStyle->{Red}];
Show[{S1,S2}]

Products of random matrices

If we multiply random 2x2 matrices, how fast does the product grow? The answer depends on the size of the random numbers. Lets take random numbers between -c and c:

F[c_]:=Module[{},randommatrix:=Table[c(2 Random[]-1),{2},{2}]; 
s=IdentityMatrix[2]; lambda=0; Do[s=randommatrix.s; u=Max[Abs[Flatten[s]]]; s=s/u; lambda=lambda+Log[u],{100}]; lambda]

Plot[F[c],{c,0.5,3}]

For which range of c, does the product appear to grow exponentially?

Nonlinear Systems

Here is an example, of how to plot a phase portrait of a nonlinear system. It is a Volterra system which is a model in population dynamics.

VectorPlot[{0.4*x-0.4*y*x,-0.1*y+0.2*x*y},{x,0,2},{y,0,2},VectorScale->0.1,VectorColorFunction->Hue,StreamPoints->20]

Fourier Series

Lets compute some Fourier series. In this example, we compute the series of an odd function so that we only have to compute the sin-series.

ClearAll[f,a,aa,s]
f[x_]:=If[x<0,1,-1];
a[n_]:=Integrate[(1/Pi) f[x] Sin[n x],{x,-Pi,Pi}];
aa=Table[a[k],{k,1,10}]
s[x_,n_]:=Sum[aa[[k]] Sin[k x],{k,1,n}];

For every n, we have a Function s[x,n] which is a Fourier approximation of the function f:

Plot[{s[x,9],f[x]},{x,-Pi,Pi}]

If you plot several plots, you observe that while the Fourier approximation converges pointwise, there is an overshoot near the discontinuity. This is called the Gibbs phenomenon.

Plot[{f[x],s[x,1],s[x,5],s[x,10]},{x,-Pi,Pi}]

Partial Differential Equations

Here is how Mathematica can solve a partial differential equation. Here it is the wave equation. The result will be given in terms of variables C[1],C[2] which can be general functions:

ClearAll[u];
DSolve[{D[u[x,t],{t,2}]==D[u[x,t],{x,2}]},u,{x,t}]

Note that the solution capabilities of Mathematica are rather limited in partial differential equations.
Here is a numerical computation of the one dimensional wave equation with inhomogeneous right hand side on the unit interval with boundary conditions zero at 0 and 1.

ClearAll[u]
s=NDSolve[{(∂_({t,2}))u[x,t]==(∂_({x,2}))u[x,t]-20 Sin[2 Pi x],u[x,0]==Sin[4Pi  x]+ Sin[6 π x],Derivative[0,1][u][x,0]==0,u[0,t]==0,u[1,t]==0},u,{x,0,1},{t,0,2}];
Plot3D[Evaluate[u[x,t]/.s[[1]]],{t,0,1},{x,0,1}]

Sound and Fourier

Instead plotting a function, you can also play it. A song is just a function.

Play[ Sin[1000  Sqrt[x]+800x^2],{x,0,2}]

If you superimpose several functions, you can hear several "melodies" together.

w=Play[Abs[Sin[200 Sin[x^2]]]+Sin[1000 x] + Cos[2000 Sqrt[x]],{x,0,6}]

You can export it as a Wave file and use it as a cellphone ring tone ....

Export["sound.wav",w,"WAV"]

The smoothness of a function depends on how nice it sounds. A discontinuous function sounds horrible. The fourier series decays slowely.

Play[ Sign[Sin[1000*x]],{x,0,2}]

Sound manipulation

Here is an involution on the linear space of all songs.

a=Import["http://www.courses.fas.harvard.edu/~math21b/labs/yesterday.wav"]; 
type=a[[1,2]];n=Length[a[[1,1,1]]];
song={Table[a[[1,1,1,n-k]]/2,{k,0,n-1}],Table[a[[1,1,2,n-k]]/2,{k,0,n-1}]};
c=Sound[SampledSoundList[song,type]]; Export["yesterday2.wav",c,"WAV"]
Show[c]

Image manipulation

Quite a few image manipulation routines are built in. The first example is a filter. Every point is averaged with weights of neighboring points.

AA1=Import["http://www.courses.fas.harvard.edu/~math21b/labs/hall.jpg"]; Show[AA1]
AA2=Dilation[ImageMultiply[LaplacianFilter[AA1,1],1],1] // ImageAdjust

The next procedure splits the picture into small pieces of size 80 pixels and then produces a graph connecting closest pieces.

im=Flatten[ImagePartition[AA1,80]];
cm = Map[Mean[Mean[ImageData[#]]] &,im];
nf=Nearest[cm ->im];
GraphPlot[Flatten[Table[Thread[im[[i]]->nf[cm[[i]],3]],{i,Length[im]}]],
VertexRenderingFunction->(Inset[#2,#,Center,0.5] &),SelfLoopStyle->None]

Here is an example on how to produce a movie morphing one picture to an other: the morph will take some time.

A1=Import["http://www.courses.fas.harvard.edu/~math21b/labs/ladybug.png"];
A2=Import["http://www.courses.fas.harvard.edu/~math21b/labs/rose.png"];Blend[{A2,A1},0.6]

Movie manipulation

The following few lines reads in a movie, identifies and colors morphological components and exports a new movie. This computation will take a while. We simplified procedures of http://blog.wolfram.com/2008/12/01/the-incredible-convenience-of-mathematica-image-processing but  import and export 
directly GIF movies.

A=Import["http://www.courses.fas.harvard.edu/~math21b/labs/geese.gif"];n=Length[A];
c=Prepend[Table[i->N[Chop[Append[Apply[List,ColorConvert[Hue[i/2.3],RGBColor]],0.5]]],{i,100}],0->{0,0,0,1}];
G[n_]:=StringTake["00000"<>ToString[n],-4]<>".gif";
T[B_]:=Image[MorphologicalComponents[ColorNegate[B],0.1]/.c,ColorSpace->"RGB"];
Export["s.gif",Table[T[A[[k]]],{k,n}],"GIF"];

Auto Tune Pitch correction

The following procedures simulates "autotune". A sound file is split into pieces, a Fourier transform performed to cut out unwanted frequencies. Then the new sound is exported. The original "autotune" (which is also called the photoshop of sound), does not merely filter frequencies but move them.

A=Import["http://www.courses.fas.harvard.edu/~math21b/labs/yesterday.wav"];M=5000;S=A[[1,1,1]];n=Length[S];
tunes=Floor[Table[440*2^(k/12),{k,-36,55}]];
filter=Table[If[MemberQ[tunes,k]  ||k>400,1,0],{k,M}];
ClearAll[S1];S1[l_]:=Table[S[[l*M+k]],{k,M}];
ClearAll[S2];S2[l_]:=FourierDCT[filter*FourierDCT[S1[l],2],3];
S3=Flatten[Table[S2[l],{l,0,Floor[n/M]-1}]];
U=Sound[SampledSoundList[{S3,S3},64000]];
 Export["yesterday1.wav",U,"WAV"]
Show[U]

Random music

Lets play some random tune:

T[x_]:=Module[{},y=x+RandomInteger[4]-2;y];
s=NestList[T,5,40];
song1=Sound[SoundNote[#,RandomChoice[{1/4,2/4}],RandomChoice[{"Cello"}]]&/@s,20]

You can export this into a midi file:

Export["song1.midi",song1,"Midi"]

Here is an example with random accords, where time is also chosen at random:

r:=RandomChoice[{4,8,12}]
g[x_]:=RandomChoice[{0.5,0.3,0.2}->{x,{x[[1]]},{x[[1]],x[[2]]}}];
accord:=Module[{},r1=r; {{1+r1,5+r1,9+r1},{1,5,9}}];
scale={1+r,5+r,9+r};
s:=g[RandomChoice[accord]]+RandomChoice[scale]
A=Table[{s,Random[Integer,3]+1},{40}];
sn[{a_,t_}]:=SoundNote[s,t];
song2=Sound[{"Piano",sn/@A},12]

Export["song2.midi",song2,"Midi"]

Motivated by Mathematica documentation we cabn use sound generation using cellular automata:

s1={"Piano",SoundNote[DeleteCases[3 Range[31] Reverse[#],0]-48,.1]&/@Transpose[CellularAutomaton[60,{{1},0},30]]};
s2={"Cello",SoundNote[DeleteCases[3 Range[31] Reverse[#],0]-48,.1]&/@Transpose[CellularAutomaton[90,{{1},0},30]]};
Sound[Join[s1,s2]]

Orchestra music

We are getting a bit off topic here. But here is a way to generate random music:

cminor={1,4,8,13};cmajor={1,5,8,13};e={1,1,1,1}; dim=Length[e];
cisminor=cminor+e;cismajor=cmajor+e;dminor=cminor+2 e;dmajor=cmajor+2 e;
disminor=cminor+3 e;dismajor=cmajor+3 e;eminor=cminor+4 e;emajor=cmajor+4 e;
fminor=cminor+5 e;fmajor=cmajor+5 e;fisminor=cminor+6 e;fismajor=cmajor+6 e;
gminor=cminor+7 e;gmajor=cmajor+7 e;gisminor=cminor+8 e;gismajor=cmajor+8 e;
aminor=cminor+9 e;amajor=cmajor+9 e;aisminor=cminor+10 e;aismajor=cmajor+10 e;
hminor=cminor+11 e; hmajor=cmajor+11 e;minors={cminor,cisminor,dminor,disminor,eminor,fminor,fisminor,gminor,gisminor,aminor,aisminor,hminor};
majors={cmajor,cismajor,dmajor,dismajor,emajor,fmajor,fismajor,gmajor,gismajor,amajor,aismajor,hmajor};
ll=Flatten[{minors,majors}];
accords=Partition[ll,dim]; 
sn[{s_,t_}]:=SoundNote[s,t];

songlength=10;
T[x_]:=x+RandomChoice[{3,5,8}]-5;
melody=NestList[T,15,songlength];
melody={1,5,6,6,5,1,3,3,1};    songlength=Length[melody];(* this is a deterministic song *)

a1=Table[sn[{accords[[melody[[k]]]],{k-1,k}}],{k,songlength}];
v2=Flatten[Table[{RandomSample[accords[[melody[[k]]]]],RandomSample[accords[[melody[[k]]]]]},{k,songlength}]];
u2=Table[r=Random[];If[r<0.4,v2[[k]]+12,If[r<0.2,v2[[k]]+24,v2[[k]]]],{k,Length[v2]}];
a2=Table[sn[{u2[[l]],{(l-1)/8,l/8}}],{l,Length[u2]}];
u3=Flatten[Table[RandomSample[accords[[melody[[k]]]]],{k,songlength}]];
a3=Table[sn[{u3[[l]],{(l-1)/4,l/4}}],{l,Length[u3]}];
u4=Flatten[Table[RandomSample[accords[[melody[[k]]]]],{k,songlength}]];
a4=Table[sn[{u4[[l]],{(l-1)/4,l/4}}],{l,Length[u4]}];

ss1=Prepend[a1,"Piano"];
ss2=Prepend[a2,"Harp"];
ss3=Prepend[a3,"Cello"];
ss4=Prepend[a4,"PanFlute"];
sss=Sound[Join[ss1,ss2,ss3,ss4]]

Random poems

Apropos creativity:  This random poem generator was adapted and modified from a Java program by Nandy Millan.
http://www.cs.bham.ac.uk/~nxm/mscPoetry/Poetry/CGPoetry.html. It generates a poem, which can be considered as a vector of word vectors. A more sophisticated poem generator would not do random choices but use a Markov chain approach to generate the sentences. Change the words or if you want the rules.

ClearAll["Global`*"];
Adjectives={placid,mellow,nice,sad,happy,tight,noisy,cold,beautiful,delicate,shining,
sweet,soft,graceful,pretty,desiring,languishing,fragile,silver,golden,
graceful,deafening,spiritless,heartbroken,sealed,stormy,carnal,frenetic,burning};
Adverbs={soon,late,desperately,passionately,slowly,softly,carefully,tenderly,gently,quickly,kindly};
Articles={the,a}; Ditverbs={gave,told,provided}; Opronouns={me,you,him,her,us,them};Spronouns={i,you,he,she,we};
Itverbs={fell,died,felt,anguished,shivered,surrendered,shouted,whispered,danced,dreamed};
Nouns={heart,girl,mouth,soul,hand,saint,maiden,love,flower,garden,moon,sun,bird,breeze,rose,kiss,misery,despair,passion,desire,warmth,fire,flame,fantasy,wing,kingdom,light,shadow};
Preps={from,for,with,in} ;Pronouns={i,you,he,she,it,we,they,me,him,her,us,them};
Tverbs={loved,hated,touched,wanted,felt,desired,cuddled,fondled,enjoyed};
Verbs={rose,coped,sang,cried,hoped,died,shouted,felt,danced,desired,dreamed};

Adjective:=RandomChoice[Adjectives];Adverb:=RandomChoice[Adverbs];Ditverb:=RandomChoice[Ditverbs];
Itverb:=RandomChoice[Itverbs];Noun:=RandomChoice[Nouns];Opronoun:=RandomChoice[Opronouns];
Prep:=RandomChoice[Preps];Pronoun:=RandomChoice[Pronouns];Spronoun:=RandomChoice[Spronouns];
OPronoun:=RandomChoice[OPronouns];Tverb:=RandomChoice[Tverbs];Verb:=RandomChoice[Verbs];
Article=RandomChoice[Articles]; Subject:=RandomChoice[{{Article,Noun},{Spronoun},{Article,Adjective,Noun}}];
Predicatelist:={{Adverb},{Prep,Article,Adjective,Noun},{Prep,Opronoun},{Article,Adjective,Noun},
{Opronoun},{Article,Adjective,Noun,Adverb},{Opronoun,Adverb},{Opronoun,Article,Adjective,Noun}};
Verblist:={{Itverb},{Itverb},{Itverb},{Tverb},{Tverb},{Tverb},{Tverb},{Ditverb}};
Predicate:=RandomChoice[Table[{Verblist[[i]],Predicatelist[[i]]},{i,Length[Verblist]}]];
Object:=RandomChoice[{{},{Adverb},{Subject},{Opronoun},{Prepsubject}}];
Verbobject:={Verb,Object}; Prepsubject:={Prep,Subject}; Subjectpredicate:={Subject,Predicate};
Continuation:=RandomChoice[{{},{"and",Subjectpredicate}}];
Sentence:=Flatten[RandomChoice[{{Subject,Verbobject,Continuation}}]]; Poem=Table[Sentence,{10}];
Poem //MatrixFormRandomChoice[Spronouns];
OPronoun:=RandomChoice[OPronouns];Tverb:=RandomChoice[Tverbs];Verb:=RandomChoice[Verbs];
Article=RandomChoice[Articles]; Subject:=RandomChoice[{{Article,Noun},{Spronoun},{Article,Adjective,Noun}}];
Predicatelist:={{Adverb},{Prep,Article,Adjective,Noun},{Prep,Opronoun},{Article,Adjective,Noun},
{Opronoun},{Article,Adjective,Noun,Adverb},{Opronoun,Adverb},{Opronoun,Article,Adjective,Noun}};
Verblist:={{Itverb},{Itverb},{Itverb},{Tverb},{Tverb},{Tverb},{Tverb},{Ditverb}};
Predicate:=RandomChoice[Table[{Verblist[[i]],Predicatelist[[i]]},{i,Length[Verblist]}]];
Object:=RandomChoice[{{},{Adverb},{Subject},{Opronoun},{Prepsubject}}];
Verbobject:={Verb,Object}; Prepsubject:={Prep,Subject}; Subjectpredicate:={Subject,Predicate};
Continuations:=RandomChoice[{{},{and,Subjectpredicate}}];
Sentence:=Flatten[RandomChoice[{{Subject,Verbobject,Continuations}}]]; Poem=Table[Sentence,{5}];
Style[Grid[Poem],"Section"]

Assignment

To get full credit for this  assignment, you have to hand in solutions to the following problems. Once you get again a grip at Mathematica, it should take you not long to do it. You have sample code for each of the problems in the text above.  Find the right code sample and modify. 

1) Get the Kirchoffmatrix of a cyclic graph  (CycleGraph) with 7 elements. Use the command “Normal” if you want to see the matrix which will be given as a SparseArray. 

2) Find the eigenvalues of the complete graph (CompleteGraph) with 7 elements.

3) Find the QR decomposition  of the  7x7 matrix which has entries A(i,j) =2^(i*j) everywhere. 

4) Produce a chart of the distribution of the largest eigenvalue of a random graph with 70 vertices where the probability of an edge added is 0.7; 

5) Find the solution of the  differential equation f''[x]-f'[x]+f[x]==x^7 Sin[x] with initial conditions  f[0]=2,;f'[0]==0.

6)  Find the first 7 terms of the Fourier series of the function  f[x_]: = ArcTan[x] on [-Pi,Pi]. 

7) Fit the primes    data=Table[{k,N[Prime[k]/k]},{k,100000}]; with functions {1,x,Log[x]}. The result will be a function which suggests the prime number theorem. 

8) Create a small Orchestra music piece. You can use the above code. Just modify the melody or instrument or more.  

In order to work on your project it is a good idea to copy paste sample code from above to do the assignment to the solution below, then delete the part above this box.
Send the notebook by email until Friday April 29 by email to knill@math.harvard.edu.

Oliver Knill,  Harvard University,  April, 2016