- Here is Mathematica code written over the last few years and optimized in 2006 to fit
in 6 lines which computes the Laplacian L = d d* + d* d for an arbitrary discrete
geometry G. The frame work used here is the language of simplicial complexes.
They are even simpler than graphs.
A simplicial complex is just a set of non-empty finite sets which have the
property that it is closed under the operation of taking non-empty finite subsets. Graphs define
simplicial complexes, the so called Whitney complex: look at all complete subgraphs and
take the set of vertices of each of this subgraphs. This is a simplicial complex. For a
triangular graph for example, the complex is G={{1,2,3},{1,2},{1,3},{2,3},{1},{2},{3}}.
There are 3 sets with one elements which are the vertices.
Then there are three sets with two elements which are the edges. Then there is one
set with three elements which are the triangles. If we write down G as above, we also have
already fixed an orientation of each of its elements.
A discrete differential form is just a function on G.
When restricted to sets of size 1, it is a 0-form, if restricted to sets of size 2, it is a
1-form, if restricted to sets of size 3, it is a 2-form, etc. The Laplacian L leaves these
k-forms invariant. On 0-forms, it is a discrete analogue of the classical Laplacian L0 = div(grad) which appears in Gauss's reformulation L0 V = 4π σ of the laws of gravity which leads to the Newton laws (as we have seen in the divergence lecture).
Then there is the Laplacian L1 which acts on 1-forms.
The equation L1 F = j is completely equivalent to the
Maxwell equations if things are done in space-time R4.
But we can do the same physics in any discrete network like our small world Gaia from Lecture 37.
Here is the code. The procedure R[n,m] generates a random complex. I added the
complex of Gaia, the discrete world [PDF] we have looked at in Lecture 37. You have computed three
eigenvalues in that lecture (they were all 2). There is no accident that 2 appeared in each sector. There is a general McKean-Singer super symmetry:
the non-zero eigenvalues on the even forms are the same
than the non-zero eigenvalues on the odd forms. See
this paper, where it was noticed to hold in the
discrete. The code below computes also the dimension of the kernels of the Li
which are called the Betti numbers. The first one b0 tells how many connectivity
components the world has. The second one b1 is zero if and only if the world is
simply connected. In that case, any vector field F which has zero curl is a gradient field.
The story is not only similar, it is exactly the same than the continuum.
(* "Cohomology in 6 lines", February 16, 2018, Oliver Knill *)
Generate[A_]:=Delete[Union[Sort[Flatten[Map[Subsets,A],1]]],1]
R[n_,m_]:=Module[{A={},X=Range[n],k},Do[k:=1+Random[Integer,n-1];
A=Append[A,Union[RandomChoice[X,k]]],{m}];Generate[A]]; G=R[5,7];
G={{1,2,3},{2,3,4},{1,2},{1,3},{2,3},{2,4},{3,4},{1},{2},{3},{4}}; (*Gaia *)
G=Sort[G]; n=Length[G]; Dim=Map[Length,G]-1;f=Delete[BinCounts[Dim],1];
Orient[a_,b_]:=Module[{z,c,k=Length[a],l=Length[b]}, If[SubsetQ[a,b] &&
(k==l+1),z=Complement[a,b][[1]];c=Prepend[b,z];Signature[a]*Signature[c],0]];
d=Table[0,{n},{n}]; d=Table[Orient[G[[i]],G[[j]]],{i,n},{j,n}];
Dirac=d+Transpose[d]; L=Dirac.Dirac; f=Prepend[f,0]; m=Length[f]-1;
U=Table[v=f[[k+1]];Table[u=Sum[f[[l]],{l,k}];L[[u+i,u+j]],{i,v},{j,v}],{k,m}];
Cohomology=Map[NullSpace,U]; Betti=Map[Length,Cohomology]
MatrixForm[Dirac] (* A matrix containing grad, curl and its adjoints*)
MatrixForm[L] (* The full Laplacian *)
Map[Eigenvalues,U] (* The eigenvalues of L_0, L_1 and L_2 *)
- Here is one possible computation of the volume π2/2
of the 3-ball x2 + y2 + z2 + w2 ≤ 1 in R4:
Integrate[1,{x,-1,1},{y,-Sqrt[1-x^2],Sqrt[1-x^2]},
{z,-Sqrt[1-x^2-y^2],Sqrt[1-x^2-y^2]},
{w,-Sqrt[1-x^2-y^2-z^2],Sqrt[1-x^2-y^2-z^2]}]
This is the sledge hammer method using Euclidean coordinates which Mathematica can do but which
we would have a bit of difficulty with without doing a couple of trig substitutions.
Better is to use spherical or Hopf coordinates. (see section 26). The Jacobean determinant is
r3 sin(2 φ)/2. Now just integrate over the toral variables θ1
and θ1 as well as φ from 0 to π/2. Even more elegant is to use that 2π2 is the surface area
(you have computed that in Homework Problem 26.5), then note that this is A(r) = 2π^2 r3
if the radius parameter r is added. Now integrate ∫0R A(r) dr = 2π2 R4/4 = 2π^2 R4/4
which is for R=1 the expression π^2/2.
- Here is the proof of the ``important formula" in the proof of Stokes
theorem:
Simplify[{Ry-Qz,Pz-Rx,Qx-Py}.Cross[{xu,yu,zu},{xv,yv,zv}]==
{Px*xu+Py*yu+Pz*zu,Qx*xu+Qy*yu+Qz*zu,Rx*xu+Ry*yu+Rz*zu}.{xv,yv,zv}-
{Px*xv+Py*yv+Pz*zv,Qx*xv+Qy*yv+Qz*zv,Rx*xv+Ry*yv+Rz*zv}.{xu,yu,zu}]
- Here is the proof that the planimeter vector field has curl 1:
s = Solve[{(x - a)^2 + (y - b)^2 == 1, a^2 + b^2 == 1}, {a, b}];
{A, B} = First[{a, b} /. s]; F = {-(y - B), x - A};
Simplify[ Curl[F, {x, y}]]
- In unit 15, we want to verify that functions solve a PDE.
Here is an example showing that the gaussian = normal distribution solves the heat equation:
f=PDF[NormalDistribution[m,Sqrt[2 t]],x];
Simplify[D[f,t]==D[f,{x, 2}]]
This actually also works
In Wolfram alpha, but you have to rerun the computation.
- For the heat equation, here is a Wolfram Alpha
verification
The direct code
f=E^(-(x-m)^2/(t))/(Sqrt[Pi*t]);Simplify[D[f,{x,2}]==4 D[f,t]]
does currently not work in Wolfram alpha.
- In unit 10, we have computed some distortion factors.
Here is how you can check it with Mathematica.
First the case of polar coordinates
{f1,f2}={r*Cos[t], r*Sin[t]};
df = {{D[f1,r], D[f1,t]},
{D[f2,r], D[f2,t]}};
Simplify[Det[df]]
And here are the spherical coordinates:
{f1, f2, f3} = {r Sin[s] Cos[t], r Sin[s] Sin[t], r Cos[s]};
df = {{D[f1, r], D[f1, s], D[f1, t]},
{D[f2, r], D[f2, s], D[f2, t]},
{D[f3, r], D[f3, s], D[f3, t]}};
Simplify[Det[df]]
We urge you still to do this calculation once by hand! You only appreciate
the machine assistance if you have gone through the simplifications yourself.
- A special page
About HW 12.
- In Units 8 and 9, it is good to be able to plot a curve. Here
is an example on how to plot a spiral
ParametricPlot[ {t Cos[t],t Sin[t]},{t,0,10Pi}]
And here is the Wolfram Alpha. Wolfram alpha also right away computes the arc length. In Mathematica, you
would do it
Integrate[ Norm[D[{t Cos[t], t Sin[t]}, t]], {t, 0, 10 Pi}]
In this case the integral was known. Numerically you get an answer with
NIntegrate[ Norm[D[{t Cos[t], t Sin[t]}, t]], {t, 0, 10 Pi}]
- In Unit 5, you deal with surfaces having holes. Here is how you draw a bagle
ContourPlot3D[(3+x^2+y^2+z^2)^2-16(x^2+y^2)==0,{x,-3,3},{y,-3,3},{z,-3,3}]
In Wolfram Alpha.
or two bagles
ContourPlot3D[((3+x^2+y^2+z^2)^2-16(x^2+y^2))((3+(x-4)^2+y^2+z^2)^2-16((x-4)^2+y^2))==0,{x,-3,7},{y,-3,3},{z,-3,3}]
In Wolfram Alpha.
- In Unit 2, we row reduced. You can get to the unique row reduced echelon form as follows.
This is an example with a 3x4 matrix:
RowReduce[{{2,2,2,2},{1,1,1,1},{3,4,6,1}}]
And here is the implementation in
Wolfram Alpha
- In Unit 5, you are asked to glue together some manifolds. Here is an
example how you can glue together two spheres. Here is how you plot
in Mathematica. If we have a product of two expressions, then one of them
has to be zero. We see therefore the union of two spheres. Now change the
constant 0 on the right with say 1 to get the two spheres merged.
ContourPlot3D[(x^2+y^2+z^2-1)((x-1)^2+y^2+z^2-1)==0,{x,-3,4},{y,0,3},{z,-2,2}]
You can also use wolfram alpha to check this:
Click here.
- Here is an example where the machine verifies some algebraic identities, like
the one appearing in Unit 4. Here is a proof that v is perpendicular to v x w.
v={v1,v2,v3}; w={w1,w2,w3};
Simplify[Cross[v,w].v==0]
If you don't have installed Mathematica, you can also verify that
with wolfram alpha:
Click here.
