Math 21b: Linear Algebra Spring 2007
Mathematica Assignment

Support

Welcome to the Mathematica computer project! In case problems with this assignment, please email math21b@fas.harvard.edu.  Mathematica can demand a lot of resources from computers. Therefore, save the notebook frequently during your work. The actual assignments are in the  magenta colored boxes. This lab is due at the last lecture of the semester. You will have to print out the assignment part with your solutions and hand them in. All problems are listed at the end. There are 7 problems. You have to have 5 correct to get full credit. All problems have sample code in this notebook somewhere.

Cells and Expressions

The notebook is subdivided into cells, bounded by a blue bracket.  This is a "text" cell, the next  is an "input" cell. We compute the dot product of two vectors:

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

To evalute an expression, click anywhere on the expresson, hold down the <Shift> key and press <Enter>.  This is 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 blue for better readability.

Basic  Computations

Mathematica allows to work with vectors and matrices in the same way as with numbers. A dot . is the matrix multiplication or dot product. The following examples should be selfexplanatory.

{{1,4},{0,1}}.{2,3}

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}}]

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

A={{1,2},{3,4}}; Eigenvalues[A]

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

A={{1,2},{3,4}}; QRDecomposition[A]

Eigenvalues

With Eigenvalues[A], we can compute the eigenvalues, with Eigensystem[A], both the eigenvalues and the eigenvectors. What happens for a random matrix? Here, we plot the complex eigenvalues of a random 50x50 matrix.

n=50;M=Table[2 RandomReal[]-1,{n},{n}];A=M+Transpose[M];
vec[z_]:={Re[z],Im[z]};EV=Eigenvalues[M];
Show[Graphics[Table[Point[vec[EV[[i]]]],{i,Length[EV]}]],AspectRatio->1]

It will be one of the homework problem to formulate a conjecture based on your experiments with larger matrices.

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 approaces 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]

Mechanical systems

Here is an example, on how one can build a more elaborate program which shows the motion of a spinner. The following code is more to illstruate how one can make relatively quickly animations of systems in a computer algebra system. It produces a movie. Run it, grab the bracket of the result and chose under the above Cell menu the entry "Animate selected graphics". It will show you how the system moves.

<<"Graphics`Shapes`";ClearAll[A,a,x,y,r,R];
MM=100;a=0.5`;A={{-1,-1},{-2,a}};
S=Transpose[Eigensystem[A][[2]]];lam1=Eigensystem[A][[1,1]];lam2=Eigensystem[A][[1,2]];
y[t_]:={Sin[lam1 t]+Sin[lam2 t],Cos[lam1 t]+Cos[lam2 t]};Scale[{u_,v_}]:={-(10+((20 (u+2))/4)),10 v};x[t_]:=Scale[S.y[t]];
r[h_]:={5,0,h};R[t_]:={{Cos[t],-Sin[t],0},{Sin[t],Cos[t],0},{0,0,1}};
ball[h_]:={AbsolutePointSize[20],RGBColor[1,0,0],Point[{0,0,h}]};rod[h_,theta_,color_]:={AbsolutePointSize[10],color,Cuboid[R[theta].r[h]-{0.8`,0.8`,0.8`},R[theta].r[h]+{0.8`,0.8`,0.8`}],Line[{{0,0,h},R[theta].r[h]}]};spring[h_]:=ParametricPlot3D[{Sin[t],Cos[t],((h t)/60)},{t,1,60},DisplayFunction->Identity];cube[h_]:={RGBColor[1,1,0],Cuboid[{-2,2,-0.01`},{2,-2,-2}]};
spinner[h_,theta_]:={spring[h],Graphics3D[{cube[h],ball[h],rod[h,theta,RGBColor[1,0,0]],rod[h,theta+((2 π)/3),RGBColor[0,1,0]],rod[h,theta+((4 π)/3),RGBColor[0,0,1]]}]};
SpinnerPict[k_]:=Show[spinner[x[((2 π k)/MM)][[1]],x[((2 π k)/MM)][[2]]],AspectRatio->1,AspectRatio->1,DisplayFunction->$DisplayFunction,Boxed->False,Axes->False,PlotRange->{{-10,10},{-10,10},{-30,0}},ViewPoint->{1.019`,-1.924`,2.59`}]Do[SpinnerPict[i],{i,1,MM}];

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.

Differential equations

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

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]

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

Fourier 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,30}];
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,5],f[x]},{x,-Pi,Pi}]

Sound

Here is how you can play a tune with Mathematica with a given wave form. In the example below, we take the Sin function as the wave form and play a tune.

PlaySong[hull_,tune_,name_,ground_]:=Module[{},u=ToCharacterCode[name];
soundfilename=FromCharacterCode[Join[u,{46,119,97,118}]];
scale[n_]:=ground*2^(n/12);beatlength=1/5;
songlength=Length[tune]*beatlength;
frequency[x_]:=tune[[1+Floor[x/beatlength]]];
song[t_]:=hull[scale[frequency[t]]*t];
P=Play[song[t],{t,0,songlength}];
Export[soundfilename,P,"WAV"];]

t1={0,0,0,2,1,2,1,2,1,2,0,0,0,4,10,7,7,5,4,2,2,0,4,7,4,0,0,0,0,0,2,0};
f1[x_]:=Abs[Sin[x]];n1="sin";g1=2000;
PlaySong[f1,t1,n1,g1]

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:

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.

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

Assignment

To get full credit for this  assignment, you have to hand in solutions to the following problems. We kept the assignment short. 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. You need to turn in correct answers to 5 problems of the 7 problems.

1) What is the distribution of the eigenvalues of a random matrix, where each entry is a random number between -1,1. Make experiments with larger matrices (sample code is above) and write down your observation in a sentence. How does the eigenvalue cloud look like in your experiments?
2) Find the determinant of the matrix  A=Table[i^j,{i,50},{j,50}] and compute its prime decomposition                         with the function FactorInteger.
3) Consider the 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}} as in the text. The matrix B=A/65 is a stochastic matrix. What do you observe, if you take powers of B? You get full credit if you write down the observation. If you have time, can you explain why for every magical square , this phenomenon takes place?
4) Find the solution of the ordinary differential equation f''[x]+f'[x]+f[x]==Cos[x]+x^2 with initial conditions  f[0]==1;f'[0]==0.
5)  Find the Fourier series of the function  f[x_]: = x^5 on [-Pi,Pi]. You need at least 20 coefficients of the Fourier series.
6) See wether Mathematica can solve the transport equation  d/dt  u(t,x) = d/dx u(t,x) with DSolve.
7) Write you own little tune with the synthesizer provided. 
        
If you find something cool during your experiments, feel free to include it also. In order to work on your project it is a good idea to save this notebook first as a different document, do the assignment directly in that notebook and print out the relevant pages at the very end on a printer. The assignments have to be printed out and turned in the last class.

Harvard University, Mathematics 21b Spring, 2007, Oliver Knill