Multivariable calculus Computer algebra Project

Oliver knill, Harvard University, knill@math.harvard.edu, MathS21a, Summer 2007

Welcome to the Mathematica project of Summer 2007! I'm are excited to use Mathematica 6 the first time in a course and we will use of some of its new features in this tutorial. This notebook will be your guide to get started. It also contains the assigment at the very end. If you work on a project, save it frequently and make backups.Have fun!.

Content in a mathematica file is organized as cells. Cells can be evaluated by grabing the bracket to the right, holding down the shift key and hitting return. Try it out in the next cell:

Plot[ Sin[x]/x,{x,-30,30}]

After evaluating the cell, an output cell has been added to the input cell. The output cell can be made interactive as in the following example:

Manipulate[Plot[Sin[x]/x,{x,-c,c}],{c,1,20}]

Lets get started. I suggest you 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 evalute it. All the cells will be evaluated.

A calculator

Anything you can do on a graphics calculator you can do with a computer algebra system, only better and more accurately. Here is an exmple to compute a numerical expression

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

It did not compute it numerically and left terms like square roots of integers intact. To get a numerical varlue or a numerical value jor a value with accuracy 200 digits, access the previous expression (called %) and

N[%,200]

Lets play with this:

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

Or look at the list of prime numbers

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

or binominal expressions

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

Some single variable calculus

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

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,Floor[n]}],{n,1,10}]

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

f[x_]:=Sin[x]*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},{0,2}}],{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}]

A Taylor expansion:

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

But lets move to multivariable calculus!

Graphs

After evaluation of the next cell, you get a graphics object, which you can manipulate with the mouse.

Plot3D[ Sin[x + y^2],{x,-2,2},{y,-2,2}]

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[2 t] Sin[3 s] {Cos[t] Sin[s],Sin[t] Sin[s],Cos[s]},{t,0,2 Pi},{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,2 Pi},{s,0,Pi},Boxed->False,Axes->False,PlotRange->{{-a,a},{-a,a},{-a,a}}],{n,-3,3},{m,-3,3}]

Parametrized Curves

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

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

r[t_]:={Cos[3 t],Sin[4 t],Cos[5 t]}; Clear[S1,S2];

S1:=ParametricPlot3D[Evaluate[r[t]],{t,0,2 Pi},Boxed->False,Axes->False]

S2[t_]:=Show[{S1,Graphics3D[{PointSize[0.05],RGBColor[1,0,0],Point[r[t]]}]}];

Manipulate[S2[t],{t,0,2 Pi}]

Polar Curves

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

Knots

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

Contour maps of a function

Manipulate[ContourPlot[c x^2-y^2-c x y^2, {x,-2,2},{y,-2,2},Frame->False,PlotPoints->40,Contours->10],{c,-1,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.

ContourPlot3D[x^2+y^2-z^2,{x,-3,3},{y,-3,3},{z,-3,3},Contours->{1},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}]

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_]:=2 x^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}]

Vector fields

Plotting vector fields needs an external library. Here is

Needs["VectorFieldPlots`"];
VectorFieldPlot[{2 x,2 y-x^2},{x,-1,1},{y,-1,1}]

Gradient fields can be plotted directly:

GradientFieldPlot[ Sin[x y],{x,-2,2},{y,-2,2}]

In three dimensions, you get an object which you can turn with the mouse:

VectorFieldPlot3D[ {x,y+z,x^2+z},{x,-2,2},{y,-2,2},{z,-2,2}]

And again there is a gradient field version

GradientFieldPlot3D[Sin[x y z],{x,-2,2},{y,-2,2},{z,-2,2}]

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

Show[ContourPlot[x^2-y^2,{x,-2,2},{y,-2,2}],Needs["VectorFieldPlots`"];

GradientFieldPlot[x^2-y^2,{x,-2,2},{y,-2,2}]]

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,2 Pi}]

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,y^2};

s[t_]:={Cos[t],Sin[t]};

LineIntegral[G,s,0,2 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=Curl[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}}]};

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[0,1,0],Disk[{0.2`,0.5`},0.1`]};

A3={RGBColor[0,0,1],Polygon[{{-0.3`,0.2`},{-0.1`,0.3`},{-0.3`,0.5`}}]};

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[1,0,0],Line[Table[f[k],{k,10,30}]]};

rp:=Point[{0.5`+RandomReal[]/2,1/2 (RandomReal[]-0.5`)}];

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]

Sound

Lets play a 2000 Herz sound for 5 seconds:

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

We can try to manipulate this:

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

Polyhedra

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["Dodecahedron","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]

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.

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

SpinnerPict[k_]:=Show[spinner[x[2 Pi k/MM][[1]],x[2 Pi 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}];

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

Image manipulation

A graphics can be saved in many file formats. Here is an example

Export["dodecahedron.jpg",Graphics3D[{PolyhedronData["Dodecahedron","Faces"]},Boxed->False],"JPEG"]

We can import the picture back as a raster image.

A=Import["dodecahedron.jpg"];

Lets darken the picture a bit.

Show[Graphics[Raster[0.8*A[[1,1]],A[[1,2]],A[[1,3]],A[[1,4]]],A[[2]],A[[3]]]]

Assignment

To get full credit for this Mathematica assignment, you have to hand in :

1) A printout of a parametric surface r (u, v) = (x (u, v), y (u, v), z (u, v)) of your choice.

2) A printout of an implicit surface g (x, y, z) = 0 of your choice.

3) A printout of a vector field in two or three dimensions of your choice.

4) Compute the arc length of the rurve r[t_] := {t^2, Sin[t^3],t^3}, where t goes from 0 to 1.

5) A printout of a graphics object of your choice, which involves discs, lines, polygons or points.

Your examples should be different from any example which appears in this notebook.

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 during the last class.

Oliver Knill, July 20, 2007, email: knill@math.harvard.edu