Machine Assistance

Laplacian of graphs Mathematica gives back a compressed matrix (SparseArray) with Kirchhhoffmatrix. With the command Normal, we get a matrix.

L=KirchhoffMatrix[CycleGraph[5]] 
L=Normal[L]
MatrixForm[L] 

                    2    -1   0    0    -1

                    -1   2    -1   0    0

                    0    -1   2    -1   0

                    0    0    -1   2    -1

                    -1   0    0    -1   2

Prime numbers The Chebyshev function:

Theta[x_] := Sum[If[Prime[k] <= x, Log[Prime[k]], 0], {k, x}]


 Music

Example 1: (one instrument)

f[k_]:=Table[Mod[k*j + 1, 24],{j,k}];
s={"Flute",Table[SoundNote[f[k],Mod[k Sqrt[2], 1]],{k,10}]};
Sound[s, {0, 5}]


Example 2: (two instruments after each other) 

f[n_, k_] := Table[Mod[n k, 24], {j, 10}]; a[t_, k_] := Mod[k t, 1];
s1 = {"Flute", Table[SoundNote[f[5, k], a[Sqrt[2], k]], {k, 10}]};
s2 = {"Violin", Table[SoundNote[f[7, k], a[Sqrt[2], k]], {k, 10}]};
Sound[{s1, s2}, {0, 5}]


Variation 3: (playing two things together) 

f[n_, k_] := Table[Mod[n k, 24], {j, 20}]; a[t_, k_] := Mod[k t, 1];
s1 = Sound[{"Flute", Table[SoundNote[f[5, k], a[Sqrt[2], k]], {k, 20}]}, {0, 4}];
s2 = Sound[{"Violin",Table[SoundNote[f[8, k], a[Sqrt[2], k]], {k, 20}]}, {0, 4}];
Sound[{s1, s2}, {0, 4}]


You find more examples
here

Operators Mathematica is an example of a high level programming language. It can also manipulate functions. The operator D, the differentiation operator is defined as D. We could reprogram it as
```
T[f_]:=f'
T[Sin]
```
So, in the last command, the operator T took as input the function Sin and returned the function Cos. Now the operator can be iterated. For example
```
T[f_]:=f';
f=Function[x,Sin[x^2]]; 
Last[NestList[T,f,5]]
```
Applies the differentiation operator five times to a function f. Now we can go a level higher and define an operator, which maps an operator to an operator. Here is an example:
```
F[S_]:=Function[f,S[S[S[S[f]]]]];
T[f_]:=f';
F[T][Sin] 
```
Take an operator S and iterates it four times. We have applied F to the differentiation operator and applied it to Sin. What we got back is Sin, because if we differentiate Sin 4 times, we get back the Sin function. In mathematics we would simply write D⁴ sin = sin and interpret this that sin is an eigenfunction of the operator D4 with eigenvalue 1.

Plotting Phase portraits Here is an example. It plots the Murray system from Unit 23:

VectorPlot[{x(6-2x-y),y(4-x-y)},{x,-1,6},{y,-1,6}, 
 VectorPoints -> 40, StreamPoints -> 120,
 StreamColorFunction->Hue]

Chaos This is for unit 24. Make sure you understand what each line does and try to test things for as many orbits (m) and lenghts of orbits (n) as possible.

(* Standard map entropy conjecture Entropy(g) >= log(g/2)          *)
pi=N[Pi]; R:=N[2*pi*Random[]];
T[{x_,y_},g_]:=Mod[{2*x-y+g*Sin[x],x},2*pi];  (* Chirikov map      *)
A[{x_,y_},g_]:={{2+g*Cos[x],-1},{1,0}};       (* Jacobean matrix   *)
Lya[{x_,y_},g_,n_]:=Module[{B=B={{1,0},{0,1}},p={x,y},t=0,a},
 Do[B=A[p,g].B;p=T[p,g];a=Abs[B[[1,1]]];If[a>1,B=B/a;t+=Log[a]],{n}];
 (t+Log[Sqrt[Tr[Transpose[B].B]]])/n];        (* Lyapunov exponent *)
entropy[g_,n_,m_]:={Sum[Lya[{R,R},g,n],{m}]/m,N[Log[g/2]]};
entropy[4.0,1000,100] (* g=4, got 100 orbits of length 1000        *)

And here is the line listed in the text of unit 24. It allows to see the orbit of a map. Here 100'000 iterations of T with the starting point of 0.3
```
T[x_] := 1.9 Cos[Exp[x]];
ListPlot[NestList[T, 0.3, 100000]]
```

Spectra of molecules

A=Normal[ChemicalData["Water", "AdjacencyMatrix"]];
L=DiagonalMatrix[Total[A]]-A; Eigensystem[L]

Isospectral graphs

s1 = UndirectedGraph[
   Graph[{71 -> 72, 72 -> 73, 73 -> 74, 74 -> 75, 75 -> 76, 76 -> 77, 
     77 -> 70, 70 -> 71, 71 -> 73, 75 -> 77, 77 -> 1, 1 -> 2, 2 -> 73,
      75 -> 3, 3 -> 4, 4 -> 71, 72 -> 76, 74 -> 5, 76 -> 6, 70 -> 7, 
     72 -> 8}]];
s2 = UndirectedGraph[
   Graph[{71 -> 72, 72 -> 73, 73 -> 74, 74 -> 75, 75 -> 76, 76 -> 77, 
     77 -> 70, 70 -> 71, 71 -> 73, 75 -> 77, 74 -> 1, 1 -> 2, 2 -> 70,
      72 -> 3, 3 -> 4, 4 -> 76, 72 -> 76, 71 -> 5, 73 -> 6, 75 -> 7, 
     77 -> 8}]];
L1 = Normal[KirchhoffMatrix[s1]]; L2 = Normal[KirchhoffMatrix[s2]];
K1 = IdentityMatrix[Length[L1]] + L1;
K2 = IdentityMatrix[Length[L2]] + L2;

Roots of xⁿ -1 You get the n'th root of 1 by writing xⁿ = 1 = e^{2 π i k} with integer k=0,1,2,...,n-1. Looks like a silly thing to do, but now we can get x = e^{2π i k/5}. These roots are all on the unit circle and form a regular polygon. Here is a way to visualize the 22'th roots of 1:
```
A = x /. Solve[x^22 == 1, x]; 
point[x_] := {Hue[Random[]], Disk[x, 0.1]};
Graphics[ Map[point, Transpose[{Re[A], Im[A]}]]]
```

N Queen problem Finding n-Queen positions:

PermutationMatrix[p_]:=Module[{n=Length[p],A},A=Table[0,{n},{n}];
   Do[A[[p[[k]],k]]=1,{k,Length[p]}]; A];
QueenConflicts[p_]:=Sum[Sum[If[Abs[p[[i]]-p[[j]]]==Abs[i-j],1,0],
   {i,j+1,Length[p]}],{j,Length[p]}];
F[n_]:=Module[{P=Permutations[Range[n]]},
    U=Flatten[Position[Map[QueenConflicts,P],0]];
    Table[PermutationMatrix[P[[U[[k]]]]],{k,Length[U]}]];
MatrixForm[F[6]]

Finding n-Super queen positions.

PermutationMatrix[p_]:=Module[{n=Length[p],A},A=Table[0,{n},{n}];
   Do[A[[p[[k]],k]]=1,{k,Length[p]}]; A];
SuperQueenConflicts[p_]:=Sum[Sum[If[Abs[p[[i]]-p[[j]]]==Abs[i-j] ||
   (Abs[i-j]==2 && Abs[p[[i]]-p[[j]]]==1) ||
   (Abs[i-j]==1 && Abs[p[[i]]-p[[j]]]==2),1,0],
   {i,j+1,Length[p]}],{j,Length[p]}];
F[n_]:=Module[{P=Permutations[Range[n]]},
    U=Flatten[Position[Map[SuperQueenConflicts,P],0]];  
    Table[PermutationMatrix[P[[U[[k]]]]],{k,Length[U]}]];
MatrixForm[F[10]]  (* PATIENCE, SOAKING !!!! *)

Quaternions Here is the code from lecture 40 in Math 22a. For some reason, companies like Adobe or Microsoft opted one or two decades ago to annoy everybody by introducing non-ASCII symbols. Copy pasting code from a PDF or word document usually is totally not working. Sometimes, the only way to fix it is to retype the characters as the evil characters (especially hyphens, dashes, quotation marks) look like the real one. You can annoy your friends by giving them code which looks like code which runs, but only does if one RETYPES the characters. Anyway, here is the code which verifies that matrix multiplication corresponds:
```
Import["Quaternions`"]
A[{x_, y_, z_, w_}] := {{ x + I*y, z + I*w}, 
                        {-z + I*w, x - I*y}}; 
Q = Quaternion[a, b, c, d] ** Quaternion[p, q, r, s]; 
Simplify[ A[{a,b,c,d}].A[{p,q,r,s}]==A[Table[Q[[k]],{k,4}]]]
```
QR decomposition
```
{Q,R} = QRDecomposition[{{1,3},{0,4}}]
A=Table[Random[Integer,5],{10},{10}]; {Q,R} = QRDecomposition[A];
Transpose[Q].Q ==IdentityMatrix[Length[Q]]
A==Transpose[Q].R
{MatrixForm[Transpose[Q]], MatrixForm[R]}
```
For some reason, Mathematica produces the transpose of Q in the first entry. Mathematica does also some strange things if the matrix is not invertible. We have defined QR decomposition to have Q the same shape than A and R always being square. We could keep that also in the case of non-invertible matrices and just add zero rows but Mathematica prefers not to do that. Try it out:
```
A=Table[1,{5},{5}]; {Q,R}=QRDecomposition[A]
A==Transpose[Q].R
```
It produces two row matrices. As mathematica shows the transpose of Q (for some really strange reasons), we have to read the result in the way that Q is a 5 x 1 matrix (a column vector) and R is a 1 x 5 matrix, (a row vector). The result is still that Q^T R is equal to A. Which makes sense. We could have put additional 0 columns and rows to keep the property that Q has the same shape than A and R is square.
Boolean matrices How can one generate all Boolean matrices of a certain size. Here it is
```
 Tuples[{0,1},{3,3}]
 MatrixForm[%]
```
Be warned that the numbers grow fast. For n=7, your computer goes up to flames. See this exhibit.
Random Matrices: In unit 6, you have to experiment with random matrices. You can plot a matrix with "Matrixplot". Here is how to plot a random matrix or a matrix A(i,j)=i+j;
```
MatrixPlot[Table[Random[],{100},{100}]]
MatrixPlot[Table[i + j, {i, 100}, {j, 100}]]
```
Matrix algebra: Computer algebra systems are fantastically proficient in computations with algebra. Here are examples.
```
   A={{3,4,5},{1,1,1},{1,2,3}}; A.A.A+Transpose[A]
   MatrixExp[{{0,-t},{t,0}}]
   A=MatrixPower[{{1,1},{1,0}},20]
   B=Inverse[A]
   B+3 A 
```
Note that A.A is not the same than A². The later just squares all elements. Similarly, 1/A is not the inverse of A, but just inverts every entry. We have to write Inverse[A] to get the inverse matrix.
```
   A={{a,b},{c,d}}; 
   A.A 
   A^2
```

Math 22b Spring 2019

22b Linear Algebra and Vector Analysis

Machine Assistance