Machine Assistance

You can get Mathematica from here. There are instructions there to download and install Mathematica. (There is a site licence for Harvard). You need to make a Wolfram account. Start here. Mathematica is a good program to get assisted by. If you have questions about the program, please ask one of us about it.

• Quantum evolution: Here we see a quantum evolution on a discrete circle:

  KirchhoffMatrix[CycleGraph[200]]; Animate[
   MatrixPlot[Abs[MatrixExp[-K*I*t]]], {t, 0, 100}]
  

• The spherical integration factor

  r={rho Sin[phi] Cos[theta],rho Sin[phi] Sin[theta],rho Cos[phi]}; 
  dr={D[r,phi],D[r,theta],D[r,rho]}; Simplify[Det[dr]]
  

• To investigate functions on a sphere, you can use color as a function value. Here is an example

  ContourPlot3D[
   x^2 + y^2 + z^2 == 4, {x, -2, 2}, {y, -2, 2}, {z, -2, 2},
   ColorFunction -> Function[{x, y, z}, Hue[Sin[5 x + 3 y] - Sin[3 z]]],
    Boxed -> False, Axes -> False, PlotPoints -> 100]
  

• To experiment with functions on the plane having two critical points only, try to start with something like

  Plot3D[ Exp[-(x-1)^2-y^2] + Exp[-(x+1)^2-y^2],{x,-2,2},{y,-2,2}]
  

  which has two maxima and a saddle. Ty to modify this to get rid of the saddle.
• For the verification of a PDE, Karen noticed that the following line does not work

  g=(Sqrt[1/t^3]*x*Exp[-x^2/(4t)])/(Sqrt[1/t]*Exp[-x^2/(4t)]+1);
  FullSimplify[D[g,t] + g D[g,x] == D[g,{x,2}]]
  

  But that when replacing Sqrt[t] with t^(1/2), it works

  g=((1/t^(3/2))*x*Exp[-x^2/(4t)])/(1/t^(1/2)*Exp[-x^2/(4t)]+1);
  FullSimplify[D[g,t] + g D[g,x] == D[g,{x,2}]]
  

  This is strange. For all that I know, Sqrt[t] and t^(1/2) should be equivalent.
• Solve a Lagrange problem with 2 variables

  f=2x^2+4 x y;  g=x^2 y;
  Solve[{D[f,x]==L*D[g,x],D[f,y]==L*D[g,y],g==1},{x,y,L}]
  

  With 3 variables:

  f=2x^2+4 x*y+z;     g=x^2 y + z;   c=1;
  Solve[{D[f,x]==L*D[g,x], D[f,y]==L*D[g,y], D[f,z]==L*D[g,z], g==c},{x,y,z,L}]
  

  With 3 variables and two constraints

  f=x y z;  g=x*y+2 y z+2 x z;  h=2x+2y+4z; d=4; c=1
  Solve[{D[f,x]==L*D[g,x] + M*D[h,x],
         D[f,y]==L*D[g,y] + M*D[h,y],
         D[f,z]==L*D[g,z] + M*D[h,z],
         g==c, h==d},{x,y,z,L,M}]
  

• Here is example code on how to compute the gradient and the discriminant

  f=x^3 y^3- x y^2;
  Grad[f,{x,y}]
  D[f,{x,2}]*D[f,{y,2}]-D[f,{x,1},{y,1}]^2
  

• To find the critical points:

  f=4 x y - x^3 y - x y^3;
  ClassifyCriticalPoints[f_,{x_,y_}]:=Module[{X,P,H,g,d,S}, X={x,y};
  P=Sort[Solve[Thread[D[f,#] & /@ X==0],X]]; H=Outer[D[f,#1,#2]&,X,X];g=H[[1,1]];d=Det[H];
  S[d_,g_]:=If[d<0,"saddle",If[g>0,"minimum","maximum"]];
  TableForm[{x,y,d,g,S[d,g],f} /.P,TableHeadings->{None,{x,y,"D","f_xx","Type","f"}}]]
  ClassifyCriticalPoints[f,{x,y}]
  

• Here is the example written down during the February 23 lecture on PDE to check that the Gaussian distribution with time as standard deviation solves the heat equation. Take this as a template for homework in Unit 14. Be careful with brackets when writing expressions. Use round brackets in algebra, the square brackets only are used for functions like Sin[x]. Also note that built-in functions are capitalized in mathematica. The language is case sensitive. The expression sin[x] does not exist. You could for example define sin[x_]:=x^2 and have your own ``sin" function. Note that := is used in definitions of functions, and == is a logical equal sign. The underline scores to the left indicate that t,x are variables. You write these underscore signs only when defining functions.

  f[t_, x_] := Exp[-x^2/(4 t)]/Sqrt[4 Pi*t];
  D[f[t, x], {t, 1}] == D[f[t, x], {x, 2}] // Simplify
  

  A big simpler is the following verification. We did not define f as a function now but as an object.

  f = Exp[-x^2/(4 t)]/Sqrt[4 Pi*t];
  Simplify[D[f,{t, 1}] == D[f,{x, 2}]]
  

  In order to evaluate f in this case, you would have to say something like

  f /. x->3 /. t->4
  

  which reads that you substitute x=3 and t=4 into the expression. To get the first derivative with respect to t, you could also simply say

  D[f,t]
  

• Here is how you can compute the Mandelbrot set. Two lines suffice. We first define a function which is either 0 or 1, then plot it. Of course this can be done much more effectively with more effective code. For example, you do not want to iterate 200 times if you already know you are not in the set.

  F[c_]:=Module[{z=c},Do[z=N[z^2+c];If[Abs[z]>2,z=10],{200}];If[Abs[z]>2,0,1]];
  ContourPlot[F[x+I*y],{x,-2,1},{y,-1.5,1.5},Axes->True,PlotPoints->100]
  

• Here is how you can draw a cylinder (with a mesh drawn) as when you look at the Schwarz Lantern.

  ParametricPlot3D[{Cos[t],Sin[t],s},{t,0,2Pi},{s,0,2}, 
   Mesh -> {7,5}, Boxed -> False, Axes -> False, MeshStyle -> {Red}]
  

• Here is the visualization from PS set 8. We draw a know and color the position with a color given by the 4th coordinate

  r[t_]:={Sin[4 t],Sin[3 t],Cos[5 t]}; color[t_] = (Cos[7 t] + 1)/2;
  Graphics3D[Table[{Hue[color[t]],Sphere[r[t],0.06]},{t,0,2 Pi,0.001}]]
  

• Here is an example how to verify that some curve is on a surface

  x=3t Cos[t]; y=2t Sin[t]; z=5t; f=x^2/9+y^2/4-z^2/25; Simplify[f==0]
  

• In order to evaluate a command in Mathematica hold down the Shift key, then (while keeping holding down shift) hit the Return key.
• Here is an example how to do a contour plot in Mathematica. Note that Mathmatica is case sensitive and the type of brackets matter. Mind the curly brackets, square brackts, round brackets.

       ContourPlot3D[ (3+x^2+y^2+z^2)^2-16(x^2+y^2)==0,{x,-3,3},{y,-3,3},{z,-3,3}]
  

  Try this to plot merge two surfaces

       ContourPlot3D[ (x^2+y^2-z^2-1)(x^2+y^2+z^2-1)==0, {x,-2,2},{y,-2,2},{z,-2,2}]
  

• 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.