Math s21a: Multivariable Calculus, Summer 2006
Mathematica Assignment

Support

Welcome to the Mathematica computer project! In case of technical problems with this assignment, please consult first the  page at http://www.courses.fas.harvard.edu/~math21a/lab.html . If the problem is not listed in the frequently asked questions there and you are still stuck, send e-mail to math21a@fas.harvard.edu. Please be aware that Mathematica is quite a complex program which demands of resources from computers, whether it is a Unix or Linux workstation, a Mac or a PC. Therefore, it is a good idea to save the notebook frequently during your work!   The actual assignment is at the very end of this notebook in a magenta colored box. If you should already be familiar with Mathematica, you can skip ahead to the end. This lab is due during the last lecture.

Cells

Mathematica notebooks are subdivided into cells. Each cell is bounded by a blue bracket at the right of the window.  Cells come in different types. The text you read now is contained in a "text" cell, while the expression below is contained in an "input" cell.

76 - 23

The different types of cells can be distinguished by the font they use: text cells useTimes font, while input cells use by default a boldface Courier font.  Text cells usually contain instructions, descriptions, or explanations.  If you want Mathematica to evaluate a mathematical expression, you must enter it in an input cell.  When creating a new cell, it will be an input cell by default.To create a new cell, position the cursor between two existing cells or after the last cell.  The cursor will then change from a vertical I-beam to a horizontal I-beam.  Click there and Mathematica will draw a horizontal line across the window. Then start typing.

Evaluating Expressions

All commands are entered into Mathematica as text.  To evalute an expression, click anywhere in the expression, hold down the <Shift> key and press <Enter>.  Mathematica will write the answer below the expression.  Your input and Mathematica's output are labeled "In[#]:" and "Out[#]:" respectively, and Mathematica draws the input in boldface to distinguish it from the plain text output.  For example, to add the following two numbers, click the sum below and press <Enter>.

3 + 8

Throughout this Mathematica assignment, there are some expressions to evaluate, along with descriptive text.  Just go through it one cell at a time, reading the text and evaluating each of the expressions one by one.  In this notebook input cells are  boxed in blue for better readability.

Basic Operators

The basic operations are + (addition), - (subtraction), * (multiplication), / (division), and ^ (exponentiation).  Mathematica obeys the usual operator precedence, first performing exponentiation, then multiplication and division, and finally addition and subtraction.  If you want to group your operations differently, use parentheses ().  Try evaluating the following expressions:

12 - 3*5

(12 - 3)*5

6/2*3

6/(2*3)

-5^2

(-5)^2

If you put two expressions next to each other separated only by a space, Mathematica multiplies them, just like in normal mathematical notation.  This can be confusing if you just put two numbers next to each other, but it comes in handy with more complicated expressions.  Here are some examples:

5 2

3 (8 - 4)

Functions and Constants

Mathematica knows many functions and constants.  Some of the functions we will be using are Sin[x], Cos[x], Exp[x] (the exponential (e^x)), Log[x] (the natural logarithm which you have seen in high school named  ln x), Sqrt[x] (square root (Sqrt[x])), and Abs[x] (the absolute value |x|).  Note that the first letter of all built-in functions are capitalized, and you use square brackets [] to apply them, as opposed to regular parentheses ().  Some useful constants that Mathematica knows are Pi (or π) and E (or e ≈ 2.71828, the base for natural logarithms).  Again, be sure to capitalize the first letter.  Alternatively, you can use the palette of templates and symbols to type in some of these operators and constants.  (If you do not see the palette, select from the menu File - Palettes - BasicInput.)  Here are some examples.

Sin[0.123]

Just like in standard mathematical  notation, if you write two things next to each other with only a space between them, Mathematica multiplies them.

Cos[2 Pi]

When you evaluated this next expression, note that Mathematica is clever enough to give you the exact answer, rather than a decimal approximation.

Sin[Pi/4]

If you prefer a numeric approximation, add //N to the end of your expression.

Sin[Pi/4] // N

or use the numerical evaluation as a function:

N[Sqrt[2]]

You can evaluate expressions with arbitrary precision. To get the first 1000 digits of π for example, try

N[Pi, 1000]

By the way, it is believed, that the digits of π occur with equal distribution. The next expression illustrates this
as well as how you can use built in functions to build more complex expressions.

ListPlot[Sort[IntegerDigits[Floor[10^1000 N[Pi, 1000]]]]]

Defining Variables and Functions

Mathematica lets you do both symbolic and numerical calculations using variables and functions.  You can use any string of letters for a variable name.  Exceptions are variable names which are already defined in Mathematica like  N. In addition to the usual variable names x, y, t, you can also use more descriptive names like theta or radius. You can assign variables a value and then use that value in further computations.  Here are some examples:

a = 5

a^2  -  3  a  +  4

b = -3

a b

Notice that if you write two variable names next to each other separated only by a space, Mathematica multiplies them, just like in normal mathematical notation.  However, you must take care to put in the space between your variable names.  Otherwise, Mathematica will think that you are talking about a single variable with a long name!  For instance, in the following cell, instead of multiplying a by b, Mathematica assumes you are using a new variable ab:

ab -  a b

To define a function for Mathematica, we write f[x_] := value, where f is the name of the function, x is the variable, and value is the expression for the function.  Note that when defining a function, you again use square brackets [] instead of parentheses. In Mathematica, this distinction is very important; parentheses () are used to group expressions, while square brackets [] are used when applying a function.  When defining a function, you use colon-equal := instead of a normal equal sign, and you put an underscore _ immediately after the variable name.  This tells Mathematica to evaluate the right side only after plugging in the value of the variable.  Note that Mathematica does not give you an output when you define a function, but once a function is defined, you can apply it to any number you want, or even to a variable or expression.  When you are finished typing out your function definition in the input cell, you must type 'enter' on the keypad or hold down 'shift' and press 'enter' on the regular key board to tell Mathematica that your definition is ready. This tells Matheamtica to read the input cell. Try out the following functions:

f[x_] := x^2 - 3  x  +  4

f[5]

f[t + 1]

Here we define a function of two variables.  Note that your variable and function names can be as long or as descriptive as you like.

DistanceFromOrigin[x_, y_]  :=  Sqrt[x^2 + y^2]

DistanceFromOrigin[3, 5]

An other way to define the distance is to use vectors and the dot product to do the same thing

DistanceFromOrigin[v_]  :=  Sqrt[v . v]

DistanceFromOrigin[{3, 4, 5}]

Graphing Functions

To graph a function, you need to tell Mathematica the function itself, the name of the dependent variable, and the domain of the function. For instance, to plot the function f(x) = (e^(-((x^2)/2))) (called a Gaussian distribution function) in the domain –3 ≤ x ≤ 3 you would type:

f[x_]  :=  Exp[-x^2/2]

Plot[f[x], {x, -3,  3}]

In the Plot command, be sure to write f[x] including the variable name. If you just write f, Mathematica will think that f is a variable name, not a function. As a shorthand, you can skip the step where you define the function separately, and include that in the Plot command itself:

Plot[Exp[-x^2/2], {x, -3,  3}]

Curves in the plane

Parameterizing a Circle

As we have seen, Mathematica's Plot[] command lets us plot the graph of a function.  However, not all curves in the plane can be described as the graph of a function.  In particular, a function can have only one value of y for each value of x.  Suppose we want to draw the unit circle, given by the equation (x^2) +  (y^2) = 1.  To draw it as a graph of a function, we first solve this equation for y to obtain y = ±(Sqrt[1 - (x^2)]).  But since there are two possible values for the square root, we cannot draw it as a graph of a single function.

Instead, we can plot the circle by treating the two coordinates x = cos(t) and y = sin(t) as functions of a third variable t called a parameter.  However, before getting to the circle, it is instructive to first plot the graphs of x and y as functions of t separately:

x[t_] := Cos[t]

Plot[x[t], {t, 0, 2 Pi} ]

y[t_] := Sin[t]

Plot[y[t], {t, 0, 2 Pi}]

Now, to plot the circle, combine them both into a single  vector-valued   function r(t) =( x(t) ,  y(t) ).  We  can represent the vector-valued  function in Mathematica as follows:

r[t_] := {x[t], y[t]}

Use Mathematica to compute the vector r(t) corresponding to various values of t:

r[0]

r[Pi/4]

r[t]

If you base these vectors at the origin, their endpoints describe a circle as t  varies from 0 to 2π.  In Mathematica, you plot this curve using the ParametricPlot[] command, which requires you to specify a vector-valued  function, the dependent variable, and the domain in the following format:

ParametricPlot[r[t],  {t, 0,  2  Pi}]

If we plot a predefined vector-valued function such as r[t], Mathematica generates a warning but proceeds with the parametric plot.  To prevent the warning message from being displayed, use the Evaluate[] version of the function r[t] to tell Mathematica to substitute the definition of r[t] into the function to be plotted, as shown in the example below.

ParametricPlot[Evaluate[r[t]],  {t, 0,  2  Pi}]

Things can be done all at once as in the following example, which shows an ellipse:

ParametricPlot[{2 Cos[t], 3 Sin[t]},  {t, 0,  2  Pi}]

To plot a polar graph  r(t) =t, we would do

r[theta_]  :=  theta

ParametricPlot[
      Evaluate[r[theta] {Cos[theta], Sin[theta]}],  {theta, 0, 2  Pi}, 
      AspectRatio -> Automatic]

Now play with the cells above to graph a spiral that winds around the origin three times. Here is an other example:

ParametricPlot[Cos[5 t] {Cos[t], Sin[t]},  {t, 0, 2  Pi}, 
      AspectRatio -> 1, Axes -> False]

Curves in space

Plotting curves in three dimensions goes very similar:

ParametricPlot3D[{Cos[t], Sin[t], t}, {t, 0, 6 Pi}]

Here is an other example, the  Slinky™. The  option PlotPoints->400 is added to the command to make Mathematica compute the value of the function at more points than usual, giving us better resolution.

r[t_] := {Cos[t] ( 2 + Cos[40  t]), Sin[40  t], 
        Sin[t] (2 + Cos[40  t])}

ParametricPlot3D[r[t],  {t, 0,  Pi}, PlotPoints -> 400]

If there is additional time, here is a task designed especially for fans of late-night infomercials.  Rotato Peeler is an amazing gadget that can peel an orange, apple, potato, or Kiwi in one piece.  It has a blade that sits on the skin of the fruit that slowly moves from the top to the bottom while the fruit turns (courtesy of a hand crank operated by various actors who probably have better things to do than peel fruit).  Assuming the fruit is a perfect sphere, we will draw the spiral path of the blade on the fruit's surface, wrapping around 10 times, using a parametric curve.

r[t_] := {Sin[t] Cos[20  t], Sin[t] Sin[20  t], Cos[t]}

ParametricPlot3D[Evaluate[r[t]], {t, 0, Pi},  
      PlotPoints -> 400]

Note: This curve involves two separate circular motions. The first is the horizontal motion in the xy-plane. The second is the vertical motion involving the radius in the xy-plane and the vertical height along the z axis. This is just half a circle going from the north pole to the south pole.  Combining the two motions, we obtain a potato peeler.

An animation

The next procedure will take as arguments a parametrized curve r[t] as well as the boundary points [a,b] on which the curve is defined. The result is a sequence of graphics.

AnimateCurve[r_, a_, b_] := 
        Module[{}, 
          S1 = ParametricPlot3D[r[t], {t, a, b}, 
              DisplayFunction -> Identity, 
              PlotPoints -> 1000]; 
S2[s_] := 
            Show[Graphics3D[{PointSize[0.02], RGBColor[1, 0, 0], 
                  Point[r[s]]}], 
              DisplayFunction -> Identity]; 
S = 
            Table[Show[{S1, S2[a + (b - a) k/30]}, 
                DisplayFunction -> $DisplayFunction], {k, 
                30}]];

Lets try this out. We define a curve R[t] in space parametrized by [0,2Pi] and then invoke the just defined procedure "AnimateCurve". After evaluating the next cell, you can animate the graphics, by grabing the bracket to the right which encloses all the pictures, then find in the menu above "Animate selected graphics". You will see a movie.

R[t_] := {Cos[3  t], Sin[4  t], Cos[5  t]};
AnimateCurve[R, 0, 2  Pi]

3D Surface Plots

Mathematica uses the Plot3D function to graph a surface defined by a function of two variables.  To plot such a graph, you need to specify the function f(x, y) and the names and range of the two independent variables.

For example, to graph the saddle f(x, y) = (x^2) - (y^2) over the domain -2 ≤ x ≤ 2, -2 ≤ y ≤ 2,  have Mathematica evaluate the following cell:

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

One of the Plot option in Mathematica is PlotPoints-> n.  This option tells Mathematica to use n points to plot the function.  The default is n = 15.  Thus, if we include PlotPoints->30 in the Plot3D function, Mathematica would plot the surface with higher resolution. Here is an example to try:

Plot3D[x^2 - y^2, {x, -2, 2}, {y, -2, 2},  
      PlotPoints -> 30]

You can look up other Plot3D options by typing

?? Plot3D

One questionmark alone gives a rough description of the function, two question marks provides you with 
more detailed information.

?Plot3D

The function f(x, y) = (Sqrt[1 + (x^2) +  (y^2)]) gives an upward-facing  hyperboloid. This is half of the two sided hyperboloid, you have seen in class. To verify this, have Mathematica evaluate the following cell:

Plot3D[Sqrt[1 + x^2 + y^2], {x, -3, 3}, {y, -3, 3}]

For practice, change the function in the Plot3D input above so as the bottom of the hyperboloid is (0, 0, 5) without the changing the shape of the hyperbolic. Tell Mathematica to plot the result. Before using Mathematica to graph the next surface, imagine what the surfaces should look like.

Plot3D[Sin[x + y], {x, -2 Pi,  2  Pi}, {y, -2 Pi, 2  Pi}, 
      PlotPoints -> 30]

Parametric Plots

To plot a surface defined by a function in polar coordinates, we will use the ParametricPlot3D function. For example, graph the surface which is given in cylindrical coordinates as z = (r^2) for 0 ≤ r ≤ 2.

ParametricPlot3D[{r Cos[t],  r Sin[t],  r^2},  {r,  0,  
        2},  {t,  0,  2  Pi}]

Earlier, we used polar coordinates to describe the surface of a paraboloid f(r, θ) = (r^2) and used the ParametricPlot3D function in Mathematica to plot the surface. You can review the picture by evaluating:

ParametricPlot3D[{r Cos[t], r Sin[t], r^2},  {r, 0, 2},  {t, 0, 
        2  Pi}]

For practice with another function, plot the top half of a sphere in polar coordinates  f(x, y) = (Sqrt[1 - (x^2) - (y^2)]).  You will need to adjust the limits of the domain to avoid taking the square root of a negative number. To make this plot, fill in the missing data in the next cell and have Mathematica evaluate the result.
Polar coordinates are useful even when the function is not radially symmetric.  Consider the saddle function f(x, y) = (x^2) - (y^2) which is in polarcoordinates x=r cos(t), y=r sin(t) given as
	\(TraditionalForm\`g(r, 
          t) = \(\((r\ cos\ t)\)\^2 - \((r\ sin  t)\)\^2 = \(\(r\^2\)(\(cos\^2\) t\  - \(sin\^2\) t) = r\^2\ \(cos(2  t)\)\)\)\). 
Here is a graph of the saddle in polar coordinates :

ParametricPlot3D[{r Cos[t], r Sin[t], r^2 Cos[2  t]},  {r, 0, 
        2},  {t, 0, 2  Pi}]

For practice, change the function in the ParametricPlot3D input above so that the resulting plot is a saddle with three instead of two pairs of upward and downward pointing directions.  Have Mathematica plot the result. This graph is also known as a Monkey Saddle. The last example illustrates waves in pond:

ParametricPlot3D[{r Cos[t], r Sin[t], 2  Cos[r]},  {r, 0, 
        6  Pi},  {t, 0, 2  Pi}]

Visualizing knots

The following procedure allows to draw a tube around a curve in space. It uses the unit tangent vector T, the normal vector N and the binormal vector B defined by the parametrized curve r(t).

\(\(\(TubePlot[c_, {t_, a_, b_}, r_, opts___] := 
        Module[{v, w, TT, NN, BB}, 
          v = D[c, t]; \ \ TT = v/Sqrt[v . v]; 
w = 
            D[TT, t]; NN = w/Sqrt[w . w]; 
          BB = Cross[TT, NN]; 
ParametricPlot3D[
            Evaluate[c + r\ Cos[s]\ NN + r\ Sin[s]\ BB], {t, a, b}, {s, 0, 
              2  Pi}, AspectRatio -> 1, Boxed -> False, 
            Axes -> False, opts]];\)\(
\)
    \)\)

Here is an example of a knot. It is a three dimensional Lissajou knot

TubePlot[{Sin[3 t], Cos[5 t], Sin[4 t]}, {t, 0, 2  Pi}, 0.1, 
      PlotPoints -> {100, 40}]

Torus knots are a class of knots which are located on the surface of a torus.

TorusKnotPlot[p_, q_, opts___] := 
        TubePlot[{Cos[t] (1 + Cos[(q/p) t]/2), 
            Sin[t] (1 + Cos[(q/p) t]/2), Sin[(q/p) t]/2}, {t, 
            0, 2 Pi p}, 0.1, opts];

Show[{TubePlot[{Cos[t], Sin[t], 0}, {t, 0, 2  Pi},  .4], 
        TorusKnotPlot[3, 5, PlotPoints -> {128, 16}]}]

Contour Lines and Level Curves

A contour line is the set of points on a surface where the height is constant.  If the surface is described by z = f(x, y), then a contour line is the set of points satisfying the equation z = f(x, y) = c.  Visually, you can think of this as the intersection of the surface with the horizontal plane corresponding to the constant function z = c.  We will use the Show function to combine the graph of the surface and the plane in order to examine the points of intersection. For example, the z = 1 contour line of the saddle f(x, y) = (x^2) - (y^2) is the intersection of the saddle surface with a plane:

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

A level curve is the projection of a contour line down onto the xy-plane.  A contour plot is a 2D graph showing multiple level curves of a function.  This gives you an alternative method of looking at a function of two variables.  If you have gone hiking, you may have seen a topographical map.  A topographical map is just a contour plot of the height function, which lets you add elevation information to an otherwise 2D map.  Another example of a contour plot is a weather map showing isobars, lines of equal barometric pressure.  In this, case the function value is not height, so it is not as natural to graph it in three dimensions.

Mathematica can draw level curves using the ContourPlot command.  To  draw the level curves of the saddle function f(x, y) = (x^2) - (y^2) for -2 ≤ x ≤ 2 and -2 ≤ y ≤ 2 have Mathematica evaluate the next cell:

ContourPlot[x^2 - y^2, {x, -2 , 2 }, {y, -2, 2}, 
      ColorFunction -> Hue]

In the next box, we form a new procedure which allows us to plot the level curves of a function together with 
the gradients. The details are not  relevant for doing this assignment, but if you are interested: we first load an additional library which allows plotting fields. The producure then forms first the contourplot, calculates
the gradient G(x,y) at each point, then forms a second graphics containing the vectors. These two graphic pictures
are combined to one plot and then displayed. The DisplayFunction is changed in order that we do not see the building blocks of the final result.

Get["Graphics`PlotField`"];
ContourWithGradient[F_, a_, b_] := 
        Module[{}, 
          S1 = ContourPlot[F[x, y], a, b, ColorOutput -> RGBColor, 
              ColorFunction -> Hue, 
              DisplayFunction -> Identity]; 
G[x_, 
              y_] := {D[F[u, y], u] /. u -> x, 
              D[F[x, v], v] /. v -> y}; 
S2 = 
            PlotVectorField[G[x, y], a, b, DisplayFunction -> Identity]; 
          Show[{S1, S2}, 
            DisplayFunction -> $DisplayFunction]];

F0[x_, y_] := x^2 - y^2 + x^3 y;
ContourWithGradient[
        F0, {x, -2, 2}, {y, -2, 2}];

Vector Fields

Once you have a vector field such as the gradient, you can have Mathematica plot it for you using the function  PlotVectorField.  Unlike the Mathematica functions we have seen so far, this function  is not in the standard Mathematica package.  To use this add-on function, we need to first load the Graphics`PlotField` package.  (Note the use of backquotes.)

Once  the package is loaded, we can use the PlotVectorField function to plot our vector fields.  Like what we have seen in ParametricPlot3D, we need to tell Mathematica to evaluate the Gradient before plotting the vector field. Tell Mathematica to evaluate the next cell in order to obtain a picture of the gradient of (x^2)-(y^2)

Get["Graphics`PlotField`"]
PlotVectorField[{2  x, 2  y}, {x, -1, 1}, {y, -1, 
        1}]

We can also draw both the level curves and gradient vector field on the same graph.  Use the Show function to combine the level curves and gradient vector field of the saddle f(x, y) = (x^2) - (y^2). To do this first have Mathematica evaluate the next cell to define f(x,y) as (x^2) - (y^2).

f[x_, y_] := x^2 - y^2

Then have Mathematica evaluate

Show[ContourPlot[x^2 - y^2, {x, -2 , 2 }, {y, -2, 2}], 
      PlotVectorField[{2  x, -2 y}, {x, -2, 2}, {y, -2, 
          2}]]

Implicit Surfaces

To draw surfaces of the form g(x,y,z)=0, we load the package ContourPlot3D first. For example, to draw an ellipsoid

Get["Graphics`ContourPlot3D`"]
f[x_, y_, z_] := x^2 + y^2 + z^2 - 1; 
    ContourPlot3D[
      f[x, y, z], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]

You can produce amazing surfaces like this. Just try to go beyond quadrics

f[x_, y_, z_] := x^4 - 2  y^2 + z^2 - 1; 
    ContourPlot3D[
      f[x, y, z], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]

Symbolic calculations

Mathematica can be used as a calculator to do symbolic or numerical calculations. For example, to 
look up integrals

Integrate[Sin[x]^5, x]

or to calculate integrals numerically:

NIntegrate[Sin[1 + Sin[x]], {x, 0, 2  Pi}]

Also higher dimensional integrals can be computed. Here is a numerical integration of a two dimensional
integral:

NIntegrate[Sin[1 + x*y], {x, -1, 1}, {y, -1, 1}]

Here is a numerical integration of a three dimensional integral:

NIntegrate[
      Sin[x] + Cos[y] + z, {x, -1, 1}, {y, -1, 1}, {z, -1, 
        1}]

Differentiation can be done in different ways:

f[x_] := Sin[ Cos[x]]

f'[x]

which is a short hand notation for

D[f[x], x]

To calculate a Taylor series of a function a a point, try

Series[Sin[x], {x, 0, 8}]

Line Integrals

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

Integrate[F[r[t]] . r'[t], {t, a, b}]

For example:

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

It might not always be possible that the integral can be evaluated symbolically. Use NIntegrate for a numerical
computation. Mathematica might give feedback about the accuracy. You can ignore these messages.

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

If we want to use line integrals many times, it might make sense to produce our own procedure LineIntegral
which takes as arguments a vector field, a curve and the integration bounds. The computation is so simple that in practice, it is actually easier to type directly the right hand side instead of remembering the new command.

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

For example:

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

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

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 2D 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 3D 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}]

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

Its your call to see whether A=B.

Graphics for illustration

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 or Postscript:

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

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

A1 = {RGBColor[1, 1, 0], 
          Disk[{0, 0}, 1, {Pi/4, 2  Pi - Pi/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  Pi t], Sin[2  Pi t]}; 
    s[t_] := 0.8*{Cos[2  Pi t], Sin[2  Pi 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 + Random[]/2, (Random[] - 0.5)/
              2}];
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]

Assignment

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

        1) A printout of a graph of a function f(x,y) of two variables of your choice.
        2) A printout of a parametric surface r(u,v) =(x(u,v),y(u,v),z(u,v)) of your choice. 
        3) A printout of an implicit surface g(x,y,z)=0 of your choice
        4) The surface area of the ellipsoid r[u_,v_]:={Cos[u] Sin[v],3 Sin[u] Sin[v], 5 Cos[v]}. 
        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 the last class.

Harvard University, Mathematics s21a Summer, 2006, Oliver Knill