Multivariable Calculus Computer Algebra Project (more)

Oliver knill, Harvard Summer School, knill@math.harvard.edu, Maths21a, Summer 2013

This file uses Mathematica 9. Please use the other project called lab.nb for the assignment. This is just an additional bonus file with more examples compiled from previous years. 
Content in a mathematica file is organized as cells. To evaluate, grab the bracket to the right, holding down the shift key and hit return.

Plot[ Sin[x^2/(1-x)],{x,-1,3}]

After evaluating the cell, an output cell has been added to the input cell. You can make it interactive as in the following example:

Manipulate[Plot[Sin[x^2/(1-x)],{x,1-c,1+c}],{c,0.1,1}]

Lets get started. Just read through the examples and evaluate the ones you are interested in which are hopefully all. If you are unpatient, grab the most outer bracket to the right which contains the entire notebook and evaluate it. All the cells will be evaluated.

The function Plot is quite flexible. Here is an example, which shows how to plot two functions at once and fill the part between them:

Plot[{Sin[x],Cos[x]},{x,0,10},Filling->True,FillingStyle->{RGBColor[1,1,0]}]

Lets now look at weather data, and plot the average temperatures of the last 66 years in Boston. Your computer will have to look up some data first:

T=WeatherData["KBOS","MeanTemperature",{{1945},{2011},"Year"}]; 
d=Table[{T[[k,1,1]],T[[k,2]]},{k,Length[T]}];f=Fit[d,{1,x},x]; Show[{Plot[f /. x->y,{y,1945,2010},PlotStyle->{RGBColor[1,0,0],Thickness[0.01]}],Graphics[{PointSize[0.02],RGBColor[0,0,1],Map[Point,d]}]},PlotRange->All]

A calculator

Anything you can do on a graphics calculator you can do with a computer algebra system:

234^10+Sin[Pi/3]^3 +Sqrt[5]

To get a numerical value or a numerical value or a value with accuracy 200 digits, access the previous expression (called %) and

N[%,200]

Lets play with this:

Manipulate[N[Sqrt[5],digits],{{digits,5,"digits"},1,1000}]

Or look at the list of prime numbers

Manipulate[Prime[n],{n,1,300000000,1}]

or binominal expressions

Manipulate[Expand[(x+y)^n],{n,1,60,1}]

Some single variable calculus

One of the most common tasks one does with a computer algebra system  is to differentiate and integrate. Here are several ways to write a derivative:

Clear[f];f[x_]:=Sqrt[x]; f'[x]

D[Sqrt[x],x]

Sqrt'[x]

Or integration

Integrate[ Sqrt[1+x^2],x]

Here we can manipulate algebraic expressions with the parameter n, the number of derivatives.

Manipulate[D[Sin[x]/x,{x,n}],{n,1,20,1}]

The derivative of a function gives the rate of change. Here is an animation which illustrates that

ClearAll["Global`*"];
f[x_]:=Sin[3x]*x; SS1=Plot[f[x],{x,-3,3},DisplayFunction->Identity];
SS2[u_]:=Graphics[{Blue,PointSize[0.03],Point[{u,f[u]}]}];
SS3[u_]:=Graphics[{Red,Line[{{u-5,f[u]-5 f'[u]},{u+5,f[u]+5 f'[u]}}]}]; 
Manipulate[Show[{SS1,SS2[u],SS3[u]},PlotRange->{{-3,3},{-3,4}}],{u,-3,3}]

Here is an experiment which explores the convergence of p-series:

Manipulate[ListPlot[Table[Sum[k^(-p),{k,n}],{n,100}],PlotRange->{0,10}],{p,0,2,0.1}]

A Taylor expansion:

Manipulate[Series[Log[1+x],{x,0,n}],{n,1,50,1}]

But lets move to multivariable calculus!

Graphs of functions of two variables

We have already seen how to graph functions of one variables. Here is an example, where the graph is filled.

Plot[x Sin[1/x],{x,-1,1},Filling->Axis,FillingStyle->Red,Background->Black,PlotStyle->{Thickness[0.01],Yellow}]

But lets move to space. With Plot3D, you can plot a function of three variables. After evaluation of the next cell, you get a graphics object, which you can manipulate with the mouse. The following surface appeared in a homework. It is an example  where the partial derivatives exist everywhere but where the derivatives are not continous:

Plot3D[ Sign[x y] Abs[x y]^(1/3),{x,-2,2},{y,-2,2}]

Also here, we can fill the region between the graph and the xy-plane, if we like. Lets fill it with blue water, with opacity 0.9 so that we can look through the water.

Plot3D[ Sin[x + y^2],{x,-2,2},{y,-2,2},Filling->Axis,PlotStyle->Green,FillingStyle->{Opacity[0.9],RGBColor[0.5,0.5,1]}]

Several graphs at once colored differently and transparent:

Plot3D[Evaluate[Table[Sin[x*y+2n],{n,3}],{x,0,3},{y,0,3},PlotStyle->Table[Opacity[0.8,Hue[i/6]],{i,4}]],Boxed->False,Axes->False,AspectRatio->1]

Here is an example on how to show a surface in a shiny golden metal form which is slightly transpartent:

Plot3D[Sin[u v/3],{u,-4,4},{v,-4,4},Mesh->None,Axes->False,BoxRatios->{1,1,0.5},PlotPoints->35,Boxed->False,PlotStyle->Directive[Opacity[0.8],Orange,Specularity[White,20]]]

Here is a way to fill the region between two graphs.

Plot3D[ {Sin[5 x* y]-4,4+Cos[3 x^2+y^2]},{x,-2,2},{y,-2,2},AspectRatio->1,Filling->Axis,FillingStyle->{RGBColor[1,1,0],RGBColor[0,1,0]},Boxed->False,Axes->False]

Contour plots

We can draw ContourPlots of functions of two variables as follows.

ContourPlot[ x^2-y^3,{x,-2,2},{y,-2,2},Contours->50,ColorFunction->Hue,ContourStyle->{Blue,Thickness[0.002]}]

If you want to see a single contour, the fastest way is the following

ContourPlot[ x^2-y^3==3,{x,-2,2},{y,-2,2},ContourStyle->{Blue,Thickness[0.01]}]

If you want to see more Contours, like 2

ContourPlot[ Sin[x^2-y^3]+Cos[x],{x,-2,2},{y,-2,2},Contours->{0,1},ContourStyle->{Red,Green}]

Here we see a bifurcation. By deforming the function the critical point changes

Manipulate[ContourPlot[c x^2-y^2-c x y^2, {x,-2,2},{y,-2,2},Frame->False,PlotPoints->40,Contours->40],{c,-1,1}]

You can even make a Contour plot interactive. Here is a cool example, which shows how to do that. You drag each of the three chargse and see the Contour plot of the charge tripol field changing:

Manipulate[ContourPlot[(1/Norm[{x,y}-C[[1]]]+1/Norm[{x,y}-C[[2]]]+1/Norm[{x,y}-C[[3]]]),{x,-3,3},{y,-3,3},Contours->20,PlotPoints->50,PlotRange->{0,10},ColorFunction->Hue],{{C,{{-1,-1},{1,1},{0,1}}},{-2,-2},{2,2},Locator},Deployed->True]

Quadrics

Manipulate[ContourPlot3D[ x^2+y^2- z^2==d,{x,-3,3},{y,-3,3},{z,-3,3},Boxed->False,Axes->False,PlotPoints->10],{d,-1,2}]

Parametrized Surfaces

ParametricPlot3D[ Cos[2t] Sin[3 s] {Cos[t] Sin[s],Sin[t]Sin[s],Cos[s]},{t,0,2Pi},{s,0,Pi},Axes->False,Boxed->False,PlotRange->All]

Here is an version which can be manipulated:

a=1;Manipulate[ParametricPlot3D[ Cos[n t] Sin[m s] {Cos[t] Sin[s],Sin[t]Sin[s],Cos[s]},{t,0,2Pi},{s,0,Pi},Boxed->False,Axes->False,PlotRange->{{-a,a},{-a,a},{-a,a}}],{n,-3,3},{m,-3,3}]

Want to be closer to the surface. Just change the View angle. Lets make it also a bit more transparent and remove the mesh

ParametricPlot3D[ Cos[2t] Sin[3 s] {Cos[t] Sin[s],Sin[t]Sin[s],Cos[s]},{t,0,2Pi},{s,0,Pi},Axes->False,Boxed->False,PlotRange->All,ViewAngle->0.2,PlotStyle->{Yellow,Opacity[0.5]},Mesh->False,Background->Black]

Here is the famous Klein bottle:

K[u_,v_]:=UnitStep[Sin[u]] {6 Cos[u](1+Sin[u])+4(1-Cos[u]/2) Cos[u]Cos[v],16 Sin[u]+4(1-Cos[u]/2) Sin[u]Cos[v],4(1-Cos[u]/2) Sin[v]}+(1-UnitStep[Sin[u]]) {6 Cos[u](1+Sin[u])-4(1-Cos[u]/2) Cos[v],16 Sin[u],4(1-Cos[u]/2) Sin[v]};
ParametricPlot3D[Evaluate[K[u,v]],{u,0,2Pi},{v,0,2Pi},Boxed->False,Axes->False,PlotStyle->{Opacity[0.5]}]

Wave fronts

Here is an example of a parametrization of a wave front as you have done in a homewok.

Manipulate[Show[{ParametricPlot3D[{x,y,x^2-y^2-1},{x,-2,2},{y,-2,2}],ParametricPlot3D[{x,y,x^2-y^2-1}+s {2 x,-2 y,-1}/Sqrt[1+4 x^2+4 y^2],{x,-2,2},{y,-2,2}]},PlotRange->{{-5,5},{-5,5},{-5,5}},Boxed->False,Axes->False],{s,-5,5}]

Parametrized Curves

Parametric curves are plotted with the same command than parametrized surfaces. There is only one paremater.

ParametricPlot3D[{(2+Cos[40t])Cos[t],Sin[40t],(2+Sin[40 t]) Sin[t]},{t,0,2Pi},AspectRatio->1]

r[t_]:={Cos[3 t],Sin[4t],Cos[5t]}; 
H1=ParametricPlot3D[Evaluate[r[t]],{t,0,2Pi},Boxed->False,Axes->False];
H2[t_]:=Show[{H1,Graphics3D[{PointSize[0.05],RGBColor[1,0,0],Point[r[t]]}]}];
Manipulate[H2[t],{t,0,2Pi}]

Polar Curves

Simple formulas can produce realistic leaf shapes:

Manipulate[PolarPlot[  Abs[m+Sin[Floor[n ]t]+Cos[Floor[k] 5 t]],{t,0,2Pi},AspectRatio->1],{n,1,10},{m,1,10},{k,1,10}]

Splines

Here is an example on how points of a curve can be manipulated:

Manipulate[Graphics[{Opacity[1,Yellow],Polygon[p[[Last[FindShortestTour[p]]]]]},PlotRange->1],{{p,RandomReal[{-1,1},{10,2}]},{-1,-1},{1,1},Locator}]

or

k=12;rr=RandomReal[{0,1},{k,2}];io=3;
Manipulate[LocatorPane[Dynamic[rr],Dynamic@(ListLinePlot[rr,InterpolationOrder->io,AspectRatio->1,Axes->False]/.Line->Polygon),LocatorAutoCreate->True],{{rr,RandomReal[{0,1},{12,2}]},ControlType->None}]

Knots

Some knots are already built in:

Manipulate[KnotData[{n,2}],{n,5,8,1}]

Implicit surfaces

Implicit surfaces are now plotted with CountourPlot. By default, a contour map is drawn.

ContourPlot3D[x^2+y^2-z^2,{x,-3,3},{y,-3,3},{z,-3,3},Boxed->False,Axes->False]

We can specify which contour we want to see by specifying the set of values. For a single value, this is the simplest

ContourPlot3D[x^2+y^2-z^2==1,{x,-3,3},{y,-3,3},{z,-3,3},Boxed->False,Axes->False]

Again this can be made into a manipulation object.

Manipulate[ContourPlot3D[x^2+y^2-z^2,{x,-3,3},{y,-3,3},{z,-3,3},Contours->{c},Boxed->False,Axes->False],{c,-4,4}]

We can change the shading

F[c_]:=ContourPlot3D[x^2-y^2-z^2+x^5c,{x,-2,2},{y,-2,2},{z,-2,2},Contours->{0.2},Boxed->False,Axes->False,Mesh->True,ColorFunction->ColorData["BrightBands"]];Manipulate[F[c],{c,0,1}]

Surfaces in spherical coordinates

Go Crimson!

SphericalPlot3D[1/(2+Sin[t+s]),{t,0,Pi},{s,0,2Pi},Boxed->False,Axes->False]

Surfaces in cylindrical coordinates

RevolutionPlot3D[{2+Cos[3t],t},{t,0,2 Pi},Boxed->False,Axes->False]

RevolutionPlot3D[{t,6Sqrt[z]-Sqrt[t]},{t,0,2Pi},{z,0,4Pi},Mesh->False,PlotPoints->50,Boxed->False,Axes->False]

Textures

Here is a bagle:

ParametricPlot3D[{Cos[u] (1.5+Cos[v]),Sin[u] (1.5+Cos[v]),Sin[v]},{u,0,2 π},{v,0,2 π},TextureCoordinateFunction->({#4,#5}&),PlotStyle->Directive[Specularity[White,300],Texture[ExampleData[{"ColorTexture","YellowMarble"}]]],Axes->False,Lighting->"Neutral",Mesh->None,Boxed->False]

You want to see what textures are built in already? Here you see:

ExampleData["ColorTexture"]

Second derivative test

The following procedure and example speaks for itself

ClassifyCriticalPoints[f_,{x_,y_}]:=Module[{X,P,H,g,d,S},X={x,y};P=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}/.Sort[P],TableHeadings->{None,{x,y,"D","f_xx","Type","f"}}]]
ClassifyCriticalPoints[4 x y-x^3 y-x y^3,{x,y}]

Lagrange multipliers

mathematica has in general no problems to solve Lagrange problems in with one constraint

F[x_,y_]:=2x^2+4 x y
G[x_,y_]:=x^2 y
Solve[{D[F[x,y],x]==L*D[G[x,y],x],D[F[x,y],y]==L*D[G[x,y],y],G[x,y]==1},{x,y,L}]

or two constraints

F[x_,y_,z_]:=x y z;
G[x_,y_,z_]:=x+y+z;
H[x_,y_,z_]:=x y+x z+y z;
Solve[{D[F[x,y,z],x]==L*D[G[x,y,z],x]+M*D[H[x,y,z],x],D[F[x,y,z],y]==L*D[G[x,y,z],y]+M*D[H[x,y,z],y],D[F[x,y,z],z]==L*D[G[x,y,z],z]+M*D[H[x,y,z],z],H[x,y,z]==750,G[x,y,z]==50},{x,y,z,L,M}]

Surface area

Given a parametrized surface r[u, v], we can calculate its area.The following example finds the surface area of the sphere.

r[u_,v_]:={Cos[u] Sin[v],Sin[u] Sin[v],Cos[v]};
l[v_]:=Sqrt[v.v];
Integrate[l[Cross[D[r[u,v],u],D[r[u,v],v]]],{u,0,2 Pi},{v,0,Pi}]

In the next example, we have a surface, where the integral can not be evaluated symbolically and where we
switch to a numerical computation of the integral with NIntegrate.

r[u_,v_]:={Cos[u/2] Sin[v],Sin[u/3] Sin[v],Cos[3 v]};
l[v_]:=Sqrt[v.v];
NIntegrate[l[Cross[D[r[u,v],u],D[r[u,v],v]]],{u,0,2 Pi},{v,0,Pi}]

Double and triple integrals

You write integrals in the same order than they appear on paper. Mathematica knows how to change the order of integration.

Integrate[ Exp[-y^2],{x,0,1},{y,x,1}]

Some integrals can not be computed in closed form like:

Integrate[ Exp[-y^Sin[y]],{x,0,1},{y,Sin[Sin[x]],1}]

you can still compute the integral numerically:

NIntegrate[ Exp[-y^Sin[y]],{x,0,1},{y,Sin[Sin[x]],1}]

Triple integrals can be done in the same way. Lets compute the volume of a sphere:

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

or the moment of inertia. Of course this would be much simpler in spherical coordinates, but Mathematica can do it the hard way too and even symbolically for a sphere of radius L:

Integrate[ x^2+y^2, {x,-L,L},{y,-Sqrt[L^2-x^2],Sqrt[L^2-x^2]},{z,-Sqrt[L^2-x^2-y^2],Sqrt[L^2-x^2-y^2]}]

Vector fields

Plotting vector fields can be plotted as follows:

VectorPlot[{-y,x+y},{x,-1,1},{y,-1,1}]

Here is a bit fancier version with streamlines:

VectorPlot[{-y,x+y},{x,-2,2},{y,-2,2},StreamPoints->{{Table[{k/10,k/10},{k,30}]}},StreamStyle->Red,StreamScale->None,VectorScale->{Large,Scaled[0.7]}]

Gradient fields can be plotted  as follows:

VectorPlot[Evaluate[D[Cos[x+y],{{x,y}}]],{x,0,1},{y,0,1}]

In three dimensions, you can look at the vector field which you can turn with the mouse:

VectorPlot3D[Evaluate[D[x+y^2-z,{{x,y,z}}]],{x,0,2},{y,0,2},{z,0,2},VectorScale->{Large,Scaled[0.4]},VectorColorFunction->"Rainbow",Boxed->False,Axes->False,Background->Black]

We all know that gradients are perpendicular to the level curves:

Show[ContourPlot[x^2-y^2,{x,-2,2},{y,-2,2},ColorFunction->"BlueGreenYellow",Frame->False,Axes->False],
VectorPlot[Evaluate[D[x^2-y^2,{{x,y}}]],{x,-2,2},{y,-2,2},VectorColorFunction->"Rainbow"],Background->Black]

New in Mathematica 7 is to see the stream lines of a vector fields. Lets show it on a nice blue background:

StreamPlot[{x^2+y^2-x,1+x-y^2},{x,-2,2},{y,-2,2},StreamColorFunction->Hue,Background->RGBColor[0,0,0.4]]

Here is an example on how to show the current wind direction data: the next command will call extrenal functions

Show[Graphics[{RGBColor[0,0.3,0],CountryData["USA","Polygon"]},Background->RGBColor[0,0,0.4]],ListStreamPlot[Table[{{x,y},Through[{Cos,Sin}[WeatherData[{y,x},"WindDirection"] ]]},{x,-125,-67,4},{y,24,50,4}],StreamColorFunction->Hue],ImageSize->800]

Line Integrals

Given a vector field F and a curve r (t), line integrals can be computed with a command like

F[{x_,y_}]:={-y/(x^2+y^2),x/(x^2+y^2)};
r[t_]:={Cos[t],Sin[t]};
Integrate[F[r[t]].r'[t],{t,0,2Pi}]

If you wanted to make this into a procedure, you could define

LineIntegral[F_,r_,a_,b_]:=Integrate[F[r[t]].r'[t],{t,a,b}]

and then give a vector field and a path as an argument:

G[{x_,y_}]:={x+2,y^2};
s[t_]:={Cos[t],Sin[ t]};
LineIntegral[G,s,0, Pi]

Flux integrals

Given a parametrized surface r[u,v],and a vector field F in space,we can compute the flux of the vector field
through the surface as follows:

r[u_,v_]:={Cos[u] Sin[v],Sin[u] Sin[v],Cos[v]};
F[{x_,y_,z_}]:={x,y^2,z-x};
Integrate[F[r[u,v]].Cross[D[r[u,v],u],D[r[u,v],v]],{u,0,2 Pi},{v,0,Pi}]

In the next example, we have a surface, where the integral can not be evaluated symbolically and where we
switch to a numerical computation of the integral with NIntegrate.

r[u_,v_]:={u^2,v^3,Sin[u v]};
F[{x_,y_,z_}]:={x,y^2,z-x};
NIntegrate[F[r[u,v]].Cross[D[r[u,v],u],D[r[u,v],v]],{u,0,1},{v,0,1}]

Greens theorem

The curl of a 2 D vector field is a scalar function.We define a Mathematica procedure which takes a vector field as an argument and returns a function :

Curl2D[G_]:=Module[{c,f},c[F_]:=Function[{x,y},D[F[{x,y}][[2]],x]-D[F[{x,y}][[1]],y]];
f[{x_,y_}]:=Evaluate[c[G][x,y]];Return[f];]

For example:

F0[{x_,y_}]:={x^7 Sin[x y],Cos[y^2]};f0=Curl2D[F0];f0[{x,y}]

Greens theorem assures that the double integral of curl (F) over a region G is the line integral of F along the boundary.Lets check that : first compute the line integral :

F[{x_,y_}]:={x^7 Sin[x y],Cos[y^2]};r[t_]:={Cos[t],Sin[t]};
A=NIntegrate[F[r[t]].r'[t],{t,0,2 Pi}]

Now compute the double integral of the curl of F over the region :

f=Curl2D[F];
B=NIntegrate[-x^8 Cos[x y],{x,-1,1},{y,-Sqrt[1-x^2],Sqrt[1-x^2]}]

These two numbers should agree.To see the error, we ask to give the result with 20 digits :

N[A-B,20]

Stokes theorem

The curl of a 3 D vector field is a vector field.We define a Mathematica procedure which takes a vector field as an argument and returns vector field which is the curl.

Curl3D[G_]:=Module[{c,f},c[F_]:=Function[{x,y,z},{D[F[{x,y,z}][[3]],y]-D[F[{x,y,z}][[2]],z],D[F[{x,y,z}][[1]],z]-D[F[{x,y,z}][[3]],x],D[F[{x,y,z}][[2]],x]-D[F[{x,y,z}][[1]],y]}];
f[{x_,y_,z_}]:=Evaluate[c[G][x,y,z]];Return[f];]

For example

F0[{x_,y_,z_}]:={x z,x y,3 x z};
G0=Curl3D[F0];G0[{x,y,z}]

Stokes theorem assures that the flux of curl (F) through a surface S is the line integral of F along the boundary.Lets check that in a special case where S is the south hemisphere.The boundary of this surfaced is a circle in the xy - plane but note the orientation!

F[{x_,y_,z_}]:={-y,Cos[y^2 z]+z,Cos[x-1]-x};
r[t_]:={Cos[t],-Sin[t],0};
A=Integrate[F[r[t]].r'[t],{t,0,2 Pi}]v

Now compute the double integral of the curl of F over the surface.Note that the parametrization is chosen such that the normal vector points outside.We compute the integral numerically :

G=Curl3D[F];
r[u,v]:={Cos[u] Sin[v],-Sin[u] Sin[v],Cos[v]};
dS[u_,v_]:=Evaluate[Simplify[Cross[D[r[u,v],u],D[r[u,v],v]]]];
g[u_,v_]:=Evaluate[G[r[u,v]].dS[u,v]];
B=NIntegrate[g[u,v],{u,0,2 Pi},{v,Pi/2,Pi}]

It is your call to see whether the two numbers agree.

Graphics

Software like Mathematica is often also used for doing illustrations in books, articles, reports or the web.Here is an example :

A1={RGBColor[1,0,0],Polygon[{{0,0},{5,0},{5,5},{0,5}}]};
A2={RGBColor[1,1,1],Polygon[{{1,2},{4,2},{4,3},{1,3}}]};
A3={RGBColor[1,1,1],Polygon[{{2,1},{2,4},{3,4},{3,1}}]};
S1=SwissFlag=Show[Graphics[{A1,A2,A3}],AspectRatio->1]

This graphics object can be exported in different formats like Tif, Gif, Jpg, PDF or Postscript :

Export["swissflag.tif",SwissFlag,"TIFF"];
Export["swissflag.gif",SwissFlag,"GIF"];
Export["swissflag.jpg",SwissFlag,"JPEG"];
Export["swissflag.pdf",SwissFlag,"PDF"]
Export["swissflag.ps",SwissFlag,"EPS"];

Besides polygons, you can use lines, discs or points.

A1={RGBColor[1,1,0],Disk[{0,0},1,{π/4,2 π-π/4}]};
A2={RGBColor[1,0,0],Disk[{0.2,0.5},0.1]};
A3={RGBColor[0,0,1],Polygon[{{-0.3,0.5},{-0.1,0.4},{-0.3,0.1}}]};
r[t_]:=1.2 {Cos[2 π t],Sin[2 π t]};
s[t_]:=0.8 {Cos[2 π t],Sin[2 π t]};
f[k_]:=If[Mod[k,2]==0,r[(2 k)/100],s[(2 k)/100]];
A4={RGBColor[0,1,0],Thickness[0.03],Line[Table[f[k],{k,10,30}]]};
rp:=Point[{0.8,0}+0.4Random[]{ Cos[2Pi t],Sin[2Pi t]} /. t->Random[]];
A5={RGBColor[0,1,1],Table[rp,{100}]};
S=Show[Graphics[{A1,A2,A3,A4,A5}],PlotRange->{{-1.3`,1.3`},{-1.3`,1.3`}},AspectRatio->1]

Or include pictures from the web:

p=CountryData["Switzerland","Polygon"];
S1=Graphics[Inset[Graphics[{Opacity[.5],RGBColor[0,1,0],p}],{0,0},Center,16]];
A1={RGBColor[1,0,0],Polygon[{{0,0},{5,0},{5,5},{0,5}}]};
A2={RGBColor[1,1,1],Polygon[{{1,2},{4,2},{4,3},{1,3}}]};
A3={RGBColor[1,1,1],Polygon[{{2,1},{2,4},{3,4},{3,1}}]};
S2=Show[Graphics[{A1,A2,A3}],AspectRatio->1];
P=Import["http://www.math.harvard.edu/~knill/oliver/oliver19.jpg","JPG"];
S3=Graphics[Inset[P,{-2,0},Center,3.5]];
Show[{S1,S2,S3},PlotRange->{{-8,8},{-8,8}}]

P1=Import["http://www.math.harvard.edu/archive/21a_fall_08/labs/dog.jpg","JPG"];
R1[a_]:=Graphics[Inset[P1,{-5,0},Center,a]];
P2=Import["http://www.math.harvard.edu/archive/21a_fall_08/labs/cat.jpg","JPG"];
R2[b_]:=Graphics[Inset[P2,{3,0},Center,b]];
Show[{R1[10],R2[10]},PlotRange->{{-8,8},{-8,8}}]

We can manipulate the pictures:

Manipulate[Show[{R1[a],R2[b]},PlotRange->{{-8,8},{-8,8}}],{a,5,15},{b,5,15}]

Sound

Lets play a 2000 Herz sound for 5 seconds:

Play[Cos[2Pi 2000 x],{x,0,5}]

We can try to manipulate this:

Manipulate[Play[Sin[c*1000 x],{x,0,1}],{c,1,10}]

You can also play Sound notes directl. Lets play some random accords with a piano

Sound[Table[SoundNote[ Random[Integer,20]+{0,4,12-Random[Integer,2]},0.2,"Piano"],{k,15}]]

The Mathematica documentation gives the following cool example of a sound generation using cellular automaton 30 :

Sound[SoundNote[DeleteCases[3 Range[31] Reverse[#],0]-48,.1]&/@Transpose[CellularAutomaton[90,{{1},0},30]]]

Objects

Lets draw a snow cube polyhedron, a semiregular polyedron. With the opacity, we can adjust, how transparent the surface is

Graphics3D[{Opacity[0.8],PolyhedronData["PentagonalHexecontahedron","Faces"]},Boxed->False]

Here are all 5 regular polyhedra matched into each other

Graphics3D[{Opacity[0.5],{PolyhedronData["Tetrahedron","Faces"],PolyhedronData["Cube","Faces"],PolyhedronData["Dodecahedron","Faces"],PolyhedronData["Dodecahedron","Faces"],PolyhedronData["Octahedron","Faces"]}},Boxed->False]

This is an example illustrating visualizing data. Here is a horse:

ListSurfacePlot3D[ExampleData[{"Geometry3D","Horse"},"VertexData"],PlotRange->All,MaxPlotPoints->40,Boxed->False,Axes->False]

Prefer a Zebra?

ListSurfacePlot3D[ExampleData[{"Geometry3D","Horse"},"VertexData"],PlotRange->All,MaxPlotPoints->40,Boxed->False,Axes->False,MeshFunctions->(5#1+2#2+#3&),MeshShading->{Blue,Yellow}]

ExampleData["Geometry3D"]

Chemical Data

There is a large amount of Chemical data available. The first time you use, it it will take some time to process. But once the picture appears you can turn the molecule.

ChemicalData["Caffeine","MoleculePlot"]

If you want to know what data are available, you can try:

ChemicalData["Properties"]

To see all molecules available evaluate:

ChemicalData[]

Try for example:

ChemicalData["Suramin","ColorStructureDiagram"]

Financial Data

Also financial data are available. Here is a plot of the Dow Jones data:

DateListPlot[{#1,#2}&@@@FinancialData["^DJI",All],Joined->True,Filling->Bottom,PlotRange->All]

What happend with the DJI since January 1, 1998?

DateListPlot[{#1,#2}&@@@FinancialData["^DJI",{1998,1,1}],Joined->True,Filling->Bottom,PlotRange->All]

What happend with the stock of the company Google between January 1, 2006 and July 1, 2011?

DateListPlot[{#1,#2}&@@@FinancialData["GOOG",{{2006,1,1},{2011,7,1}}],Joined->True,Filling->Bottom,PlotRange->All]

Objects

ListSurfacePlot3D[ExampleData[{"Geometry3D","Horse"},"VertexData"],Lighting->"Neutral",MaxPlotPoints->50,Mesh->False,PlotStyle->{Specularity[White,20],Texture[ExampleData[{"ColorTexture","GiraffeFur"}]]},
TextureCoordinateFunction->(Normalize[{#1,#2,#3}]&),
Boxed->False,Axes->False]

Want to see what is built in already?

ExampleData["Geometry3D"]

S1=ExampleData[{"Geometry3D","Horse"}];
S2=ExampleData[{"Geometry3D","Cow"}]; 
S3=ExampleData[{"Geometry3D","HammerheadShark"}]; 
S4=ExampleData[{"Geometry3D","Triceratops"}]; 
S5=ExampleData[{"Geometry3D","StanfordBunny"}]; GraphicsGrid[{{S1,S2,S3,S4,S5}}]

Linguistic

Look up the meaning of a wrod

WordData["sing"]

For solving crossword puzzles, the following is handy: lookup all words which start with "o" contain "er" and end with "r"

DictionaryLookup["o"~~__~~"er"~~__~~"r"]

Lets define a function, which gives us words close to a given word :

nf=Nearest[DictionaryLookup[]];

Which words are in distance 4 from "Harvard"?

nf["Harvard",4]

Lets see what is close to  the word "Calculus" and illustrate this

w=Nest[Flatten[Thread[#->nf[#,4]]&/@(Last/@#)]&,{""->"calculus"},3];
Show[GraphPlot[w,VertexLabeling->True,AspectRatio->1]]

The Gettyisburg Address :

ExampleData[{"Text","GettysburgAddress"}]

See the words of the Gettysburg address:

w=ExampleData[{"Text","GettysburgAddress"},"Words"]

Show their frequency :

w1=Tally[w]

Lets see only the numbers, sorted  and reversed

A=Reverse[Sort[Last/@w1]];

ListPlot[A, Joined -> True, Filling -> Axis]

What texts are available?

ExampleData["Text"]

We also can get texts from the web : lets get Mark Twains "Adventures of Huckleberry Finn"

w=Import["http://www.gutenberg.org/files/76/76.txt","Words"];

The 1001'th word is

w[[1001]]

The word distribution. A bit more than 1/2 of the words appear once.

ListPlot[Sort[Last /@ Tally[w]]]

Websites

You can grab websites and read it back to you:

Clear[S];S=Import["http://www.cnn.com"];Speak[S]

Animations

The following animation is a mechanical system. Run the cell and you will see an animated picture. You can turn around the picture while it moves.

ClearAll[S,U];
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}];
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+2Pi/3,RGBColor[0,1,0]],rod[h,theta+4Pi/3,RGBColor[0,0,1]]}]};

SpinnerPict[k_]:=Show[spinner[x[2Pi k/MM][[1]],x[2Pi k/MM][[2]]],AspectRatio->1,Boxed->False,Axes->False,PlotRange->{{-10,10},{-10,10},{-30,0}},ViewPoint->{1.019,-1.924,2.590}];

S=Table[SpinnerPict[i],{i,1,MM}];

U=Animate[S[[Floor[k]]],{k,1,MM}]

Image manipulation

Lets get a picture from the web:

P=Import["http://www.math.harvard.edu/archive/21a_fall_08/labs/oliver.jpg","JPG"]

The following matrix of color vectors contains the image data:

A=P[[1]];

We get the (red,green,blue) color values for the pixel at position (3,4) as follows

A[[3,4]]

Lets do some woodoo on Oliver and find the boundaries:

AA2=Dilation[ImageMultiply[LaplacianFilter[P,1],1],1]; Show[AA2]

ColorNegate[P]

Image[MorphologicalComponents[ColorNegate[P],0.5]]

Demonstrations

By opening a new demonstration (under File -> Demonstrations) you are able to produce a program which can be run with Mathreader without Mathematica.You can convert your project on this link : http : // www.wolfram.com/solutions/interactivedeployment/publish
Here are some additional inspirations :

Manipulate[Graphics[{v=p2-p1;w=p3-p1;e={0.2,-0.2};
d=(v[[2]] w[[1]]-v[[1]] w[[2]])/Sqrt[v.v];
p4=p3+d*{-v[[2]],v[[1]]}/Sqrt[v.v];
{RGBColor[0,0,0],FontSize->40,Text["Vector Projection",{0,1.8}]},{RGBColor[0,1,0],Thickness[0.001],Dynamic[Line[{p1-100*v,p2+100*v}]]},{RGBColor[0,0,1],Disk[p1,0.1]},{RGBColor[0,1,0],Disk[p2,0.1]},{RGBColor[1,0,0],Disk[p3,0.1]},{RGBColor[1,1,0],Disk[p4,0.1]},{RGBColor[0,0.4,0],Thickness[0.01],Arrow[{p1,p2}]},{RGBColor[0.5,0,0],Thickness[0.002],Arrow[{p1,p3}]},{RGBColor[1,0.8,0],Thickness[0.012],Arrow[{p1,p4}]},Locator[Dynamic[p1],ImageSize->40],Locator[Dynamic[p2],ImageSize->40],Locator[Dynamic[p3],ImageSize->40]},PlotRange->{{-2,2},{-2,2}}],{{p1,{-1,-0.3}},{-1,-1},{1,1},ControlType->None},{{p2,{1,-0.5}},{-1,-1},{1,1},ControlType->None},{{p3,{-0.3,1.2}},{-1,-1},{1,1},ControlType->None}]

Example 1: The following example is given in the Mathematica Guidbook by Michael Trott : It shows how one can build a Penrose 
Tribar :

Manipulate[With[{b=N[1/7]+A,c=1.12,d=1.12,e=1.12,f=0.9333},p1=Polygon[{{b,0,b},{b,-b,b},{1+b,-b,b},{1+b,-b,b},{1+b,0,b}}];
p2=Polygon[{{b,-b,b},{b,-b,0},{1+b,-b,0},{1+b,-b,b},{1+b,-b,b}}];
p3=Polygon[{{0,0,b},{0,-1,b},{b,-1,b},{b,0,b}}];
p4=Polygon[{{0,0,0},{0,0,b},{0,-1,b},{0,-1,0}}];
p5=Polygon[{{0,-1,b},{b,-1,b},{b,-1,c},{0,-1,d}}];
p6=Polygon[{{0,-1,0},{b,-1,0},{b,-1,b},{0,-1,b},{0,-1,b}}];
p7=Polygon[{{0,-1+b,b},{0,-1,b},{0,-1,e},{0,-1+b,f}}]];
Show[Graphics3D[{RGBColor[1,1,0],{p1,p2,p3,p4,p5,p6,p7}}],PlotRange->All,ImageSize->{600,600},ViewPoint->{-12.625,-10.9375,14.21},Boxed->False],{A,-0.2,0.2}]

Example 2

h[x_]:=If[Abs[x]<1,1,-1];
f[{{a_,b_},{v_,w_}}]:={{a+h[a] v,b+h[b] w},{h[a] v,(h[b]-0.0001) w-0.001}};
DynamicModule[{c={}},r:=0.03*(Random[]-1/2);
EventHandler[Graphics[{RGBColor[1,0,0],PointSize[0.05],Point[Dynamic[c=Map[f,c];Map[First,c]]]},PlotRange->{{-1,1},{-1,1}},ImageSize->{600,600}],"MouseDown":>(AppendTo[c,{MousePosition["Graphics"],{r,r}}])]]

Example 2

Manipulate[If[r==1,S=ContourPlot3D[x^2*P[[1]]+y^2*P[[2]]-z^2==1,{x,-1,1},{y,-1,1},{z,-1,1},Boxed->False,Axes->False]];
If[r==2,S=ContourPlot3D[x^2*P[[1]]+y^2*P[[2]]+z^2==1,{x,-1,1},{y,-1,1},{z,-1,1},Boxed->False,Axes->False]];
If[r==3,S=ContourPlot3D[x^2*P[[1]]+y^2*P[[2]]+z==1,{x,-1,1},{y,-1,1},{z,-1,1},Boxed->False,Axes->False]];
If[r==4,S=ContourPlot3D[x^2*P[[1]]+y^2*P[[2]]==1,{x,-1,1},{y,-1,1},{z,-1,1},Boxed->False,Axes->False]];S,{{r,1,"surface:"},{1->"hyperboloid",2->"ellipsoid",3->"paraboloid",4->"cylinder"}},Control[{{P,{1,1}},{0,0},{4,4},ImageSize->{350,200}}]]

Example 3

Manipulate[If[r==1,S=StreamPlot[{P[[1]] y,P[[2]] x},{x,-2,2},{y,-2,2}]];
If[r==2,S=VectorPlot[{P[[1]] y,P[[2]] x},{x,-2,2},{y,-2,2}]];S,{{r,1,"type:"},{1->"stream lines",2->"vectors alone"}},Control[{{P,{1,0}},{-1,-1},{1,1},ImageSize->{350,200}}]]

Example 4

Manipulate[Graphics[{{x1,y1}=p1;{x2,y2}=p2;{x3,y3}=p3;
R=2*(x3*(y1-y2)+x1*(y2-y3)+x2*(-y1+y3));
m1=(x3^2*(y1-y2)+(x1^2+(y1-y2)*(y1-y3))*(y2-y3)+x2^2*(-y1+y3))/R;
m2=(-(x2^2*x3)+x1^2*(-x2+x3)+x3*(y1^2-y2^2)+x1*(x2^2-x3^2+y2^2-y3^2)+x2*(x3^2-y1^2+y3^2))/R;
center={m1,m2};radius=Sqrt[(center-p1).(center-p1)];
{RGBColor[1,0,0],Dynamic[Disk[center,0.07]]},{RGBColor[0,0,1],Dynamic[{Disk[p1,0.1],Disk[p2,0.1],Disk[p3,0.1]}]},{RGBColor[1,0,0],Thickness[0.007],Dynamic[Circle[center,radius]]},{RGBColor[0,1,0],Thickness[0.004],Dynamic[Line[{p1,p2,p3,p1}]]},Locator[Dynamic[p1],ImageSize->40],Locator[Dynamic[p2],ImageSize->40],Locator[Dynamic[p3],ImageSize->40]},PlotRange->{{-2,2},{-2,2}},ImageSize->{600,600}],{{p1,{1.1,0.6}},{-1,-1},{1,1},ControlType->None},{{p2,{-0.9,0.5}},{-1,-1},{1,1},ControlType->None},{{p3,{-0.3,1.2}},{-1,-1},{1,1},ControlType->None}]

Example 5

Manipulate[Graphics[{Q1={0,0};Q3=Q2+Q4;
c=Floor[100*(Q2[[1]]*Q4[[2]]-Q4[[1]]*Q2[[2]])];
{RGBColor[0,0,1],Disk[Q2,0.1]},{RGBColor[1,0,0],Disk[Q4,0.1]},{RGBColor[0,1,0],PointSize[0.04],Point[Q1]},{RGBColor[1,0,0],Thickness[0.01],Dynamic[Arrow[{Q1,Q2}]]},{RGBColor[0,0,1],Thickness[0.01],Dynamic[Arrow[{Q1,Q4}]]},{RGBColor[1,1,0],Dynamic[Polygon[{Q1,Q2,Q3,Q4,Q1}]]},Locator[Dynamic[Q2],ImageSize->40],Locator[Dynamic[Q4],ImageSize->40],{FontSize->40,Text["Cross Product",{0,1.8}]}},PlotRange->{{-2,2},{-2,2}}],{{Q2,{1,0}},{{-1,-1},{1,1}},ControlType->None},{{Q4,{0,1}},{{-1,-1},{1,1}},ControlType->None}]

Assignment

See the seperate notebook lab.nb.

Oliver Knill, July 15, 2013, email: knill@math.harvard.edu