Mathematical Modeling & Mathematica

Lab 6:  Graphing in 3-Space

John E. Cruthirds and L. Richard Hitt
Department of Mathematics and Statistics
University of South Alabama

Objective

This lab explores the capabilities of Mathematica  for graphing curves and surfaces in 3-space.  In addition to graphing functions of 2-variables like f(x,y), parametric curves, parametric surfaces, and surfaces in cylindrical and spherical coordinates will also be explored.

Introduction

In Lab #1 you explored some of the graphing capabilities of Mathematica .  In this lab we will take a closer look at the capabilities of Mathematica  with regard to graphing in 3-space.  The primary tool for graphing functions of 2-variables is the Mathematica  function Plot3D.  The next section studies this function along with various of its options.   After skill with Plot3D is developed, successive sections examine graphing in cylindrical and spherical coordinates, graphing parametric surfaces, coutour graphs, and animations of such graphics.

Using Plot3D

As you probably recall from Lab 1, the Plot3D function is used to graph functions of two variables as follows.

g = Plot3D[Sin[x y], {x, -2, 2}, {y, -2, 2}]

There are many options to the Plot3D function.  We can list the options together with their default values.

Options[Plot3D]

Fortunately, the defaults chosen by Mathematica  gave a reasonable rendition of the above graph.  You will find that you do not have to become familiar with very many of the options of Plot3D to be effective in graphing in 3-space.  We can illustrate the effect of a few of the options by making some adjustments on the previous graph.  Since the options we are adjusting below do not affect the sample points computed by Plot3D in the previous graph, we can just pass the options along to Show.

Show[g, AxesLabel -> {"x", "y", "z"}, 
      Ticks -> {Automatic, Automatic, {-1, 0, 1}}, 
      PlotLabel -> "Sin[x y]", 
      FaceGrids -> {{0, 1, 0}, {-1, 0, 0}}]

The Plot Algorithm

To gain some insight into the method used by the Plot3D function to construct graphs, let's consider what would happen if we graphed the same function from the previous section using a larger region in the xy-plane.

Plot3D[Sin[x y], {x, -5, 5}, {y, -5, 5}, 
      AxesLabel -> {"x", "y", "z"}, 
      Ticks -> {Automatic, Automatic, {-1, 0, 1}}, 
      PlotLabel -> "Sin[x y]"]

What a mess.  In looking at the graph, it is apparent that there is no adaptive algorithm operating with Plot3D which automatically increases the density of sampled points whenever sharp bends are encountered like there is with the Plot function for plane graphs.  Instead, what happens is that there are PlotPoints number of equally spaced points sampled in both the x and y directions.  From the above list of options, we see that PlotPoints has a default value of 15.  Thus, there are 225 sample points used to construct the above graph.  From the computed graph, we see that the number of sample points needs to be increased.  We can do this with the following.

Plot3D[Sin[x y], {x, -5, 5}, {y, -5, 5}, 
      PlotPoints -> 51, AxesLabel -> {"x", "y", "z"}, 
      Ticks -> {Automatic, Automatic, {-1, 0, 1}}, 
      PlotLabel -> "Sin[x y]"]

Now a much clearer picture of the function emerges.

Miscellaneous Options for Plot3D

From the list of options for Plot3D in the previous subsection, you see that there are many which are available.  Nevertheless, we have been able to obtain reasonably good graphs of functions by changing the defaults of only a few of the options.  In this subsection, we will consider some of the options we have not used yet.

Let's begin with the default graph of a particular cubic polynomial and see what we might want to change to improve it.  The surface occuring as the graph of this funciton is often called a "monkey saddle".  Can you see why?

Clear[g]g = Plot3D[(x^3) - 3 (y^2) x, {x, -1, 1}, {y, -1, 1}]

That's not bad as it is.  But suppose we want to eliminate the bounding box and the axes.  Additionally, we can try to uncompress the z-axis relative to the other two axes, and modify the perspective so we can see more of the underside of the surface.  To avoid recomputing all the points, we use the Show command.

Show[g, Boxed -> False, Axes -> None, 
      BoxRatios -> {1, 1, 1.5}, ViewPoint -> {1.2, -2, 1.2}]

In order to understand the ViewPoint value, you have to understand that Mathematica  uses relative (not absolute) coordinates for this option.  In particular, Mathematica  treats the bounding box as being one relative unit long in each dimension.  So the center of the bounding box has relative coordinates (.5,.5,.5).  Now you can follow how the adjustment on the ViewPoint option was made.  When you need to experiment with the ViewPoint, a convenient way to do this is to choose the 3D ViewPoint Selector from the Prepare Input selection of the Action menu.

The BoxRatios option is similar to the AspectRatio option used in two-dimensional plots.  It determines the ratios of the side lengths of the box bounding the graphic.

Contour and Density Plots

A contour graph (sometimes called a level curve graph) of a function  f(x,y)  is a graph of curves in the xy-plane of the form  k = f(x,y)  for various constant values of  k .  For example, weather maps often show isobaric (constant pressure) or isothermal (constant temperature) curves, which are the level curves for the barametric pressure or temperature function, respectively.  In Mathematica, contour graphs are produced with the ContourPlot function.

If we entered the command

ContourPlot[ x^3 - 3 y^2 x, {x,-1,1}, {y,-1,1} ]

a contour plot of the previous graphic would be produced.  However, because the information stored by Mathematica  for Plot3D (which generates what Mathematica  calls a SurfaceGraphic object) and ContourPlot (which generates a ContourGraphic object) are so similar, we can actually convert directly from one form to the other, thereby saving having to recompute all the sample values of the function.  The following command does this, and produces the exact same result (only faster) as if you had entered the ContourPlot command above.

Show[ContourGraphics[g]]

In addition to producing several level curves for the function, Mathematica  also uses shading to indicate the relative z-values.  Compare the shading used in the contour plot with the actual graph to see how the shading works.

Of course, the ContourPlot function has many options which could be viewed with the Options[ContourPlot] command, but usually the default settings are adequate.  If you need to change the number of contours which are graphed, that can be done with the Contours option (the default is 10).  You can also insist that Mathematica  try to smooth the contours by setting ContourSmoothing -> Automatic (the default is None so that no smoothing of curves is attempted).  Be warned, though, that smoothing contours can take both time and memory, so you should have plenty of both before attempting smoothing.

A density graph of a function of two variables is similar to a contour graph.  Shading which indicates the relative z-values is done in rectangular regions in the xy-plane, but no contour curves are drawn.  So a density graph can be computed more quickly than a contour graph.  The following command converts the SurfaceGraphics object stored in the variable g to a DensityGraphics object.  As with the above discussion on contour plotting, the same result could be obtained using the DensityPlot command, except that all the values of the sample points would have to be recomputed.

Show[DensityGraphics[g]]

Parametric Plots

In this section we will examine various curves and surfaces which can be described using equations involving parameters.

Parametric Curves in 3-Space

Mathematica  can graph curves which are described parametrically in the plane or in 3-space .  For 3-space curves, you need to describe each of x, y, and z as a function of a parameter, say t.  For example, if we have

(x, y, z) = (t, t2, t3)

we can graph the resulting parametric curve in 3-space with the Mathematica  function ParametricPlot3D as follows.

ParametricPlot3D[{(t^2), t, (t^3)}, {t, -1.5, 1.5}]

This curve is called a "twisted cubic".

Notice that the bounding box helps provide a sense of depth to the curve.  Without the box, all you see is a curve in the plane of the computer screen.  Even with the box, however, you still can't detect the nuances of the curve from a single projection into the viewing plane.  To get a better idea of what the curve is doing, you could, for example, look at several different perspectives.  Alternatively, you look at the unit tangent and normal vectors moving along the curve.  Both of these techniques are addressed in the Animation section later in the lab.

Parametric Surfaces

Surfaces can often be described using equations for the three rectangular coordinates which involve two parameters.  In the first lab, you saw such an example where a hyperboloid of one sheet was graphed using the cylindrical coordinate parameters z and θ.  We will begin with other examples which use cylindrical and spherical coordinates.  Then we will look at a completely different type of parametric surface.

Cylindrical and Spherical Coordinates

We will start with the problem of graphing a sphere.  For definiteness, we use the unit sphere centered at the origin in 3-space.

x2 + y2 + z2 = 1

Since the sphere is not the graph of any function f(x,y), we cannot use Plot3D directly to produce the desired graph.  We can't even use Plot3D to graph the upper or lower hemisphere since Plot3D requires a rectangular region in the xy-plane over which to work.  What we need to do is describe the sphere parametrically and utilize the Mathematica  function ParametricPlot3D.  When each of the variables x, y, and z are described as functions of two parameters, ParametricPlot3D produces a surface graph, whereas, in the previous subsection, the same Mathematica  function produced a curve.

To describe the sphere parametrically, we can do this with either cylindrical or spherical coordinates, among other methods.

In the case of cylindrical coordinates (r, θ, z), we have a slight complication in that the cylindrical coordinate z must be expressed as two different functions of r and θ, i.e., z is not uiquely determined by r and θ.  This leads to the following use of ParametricPlot3D.

Clear[x, y, z, r, theta, solution, f1, f2]x = r Cos[theta]; y = r Sin[theta]; solution = Simplify[Solve[(x^2) + (y^2) + (z^2) == 1, z]]; f1 = {x, y, 
        solution[[1, 1, 2]]}; f2 = {x, y, 
        solution[[2, 1, 2]]}; ParametricPlot3D[Evaluate[{f1, f2}], {r, 0, 1}, {theta, 0, 2 π}, 
      BoxRatios -> {1, 1, 1}]

This gives an acceptable graph of the sphere.  But the graph looks a little strange in that the latitude lines do not appear evenly spaced around the sphere.  Do you see why?

This problem can be remedied using spherical coordinates, where things become much simpler as the name of the coordinate system might suggest.

Clear[x, y, z, rho, phi, theta]x = rho Sin[phi] Cos[theta]; y = rho Sin[phi] Sin[theta]; z = rho Cos[phi]; rho = 1; ParametricPlot3D[{x, y, z}, {phi, 0, π}, {theta, 0, 2 π}, 
      BoxRatios -> {1, 1, 1}]

This graph of the sphere is more like what we are accustomed to seeing.  The spherical rectangles appear nice and regular.

For our next example, we will graph the function

x y        
f(x,y) = ---------
        x2 + y2

Note that f(x,y) is not defined at the origin in the plane.  In fact, f(x,y) has an essential discontiniuty there (i.e., the limit of f(x,y) fails to exist at the origin). The function also has the interesting property that its value along any straight line through the origin in the plane is constant except at the origin, where it is undefined.  Mathematica  can easily check this for us.

Clear[f, x, y]f[x_, y_] := ((x y)/((x^2) + (y^2)))y = m x; Simplify[f[x, y]]

So the values of f[x,y] along any line y = mx depend only on m, not on x or on y (except at the origin).

Of course, the line given by  x = 0  isn't covered in the above check, but  f[0,y] = 0  for any nonzero value of y, so f[x,y] is constant along this line as well.  We can plot the function now using Plot3D.

Plot3D[f[x, y], {x, -1, 1}, {y, -1, 1}, 
      AxesLabel -> {"x", "y", "z"}, 
      Ticks -> {Automatic, Automatic, {- .5, 0,  .5}}, 
      BoxRatios -> {2, 2, 1}]

The type of wild behaviour you see around the origin is typical of what happens when Mathematica  attempts to graph a surface through this type of essential discontinuity.  Increasing the value of PlotPoints won't significantly help -- the same problems in the graph will persist, only at a smaller scale.  To get a different ruling on the surface, we can switch to another coordinate system.  For example, in cylindrical coordinates, we get the following.

Clear[x, y, r, theta, f]; f[x_, y_] := ((x y)/((x^2) + (y^2)))x = r Cos[theta]; y = r Sin[theta]; z = Simplify[f[x, y]]

We see that the values of the function depend only on theta, not on r.  This again shows that the values are constant along straight lines in the plane containing the origin.  Now we can graph our function using ParametricPlot3D.

ParametricPlot3D[{x, y, z}, {r, 0, 1}, {theta, 0, 2 π}, 
      BoxRatios -> {1.5, 1.5, 1}, PlotPoints -> {5, 33}, 
      ViewPoint -> {-1.5, -4, 2}, 
      AxesLabel -> {"x", "y", "z"}]

A Moebius Band

For an example of a surface described parametrically which does not use parameters from any of the standard coordinate systems, we will start with a circle of radius 4 centered at the origin of the xy-plane.  At each point of the circle, we will attach a line segment of length 1 so that it is orthogonal to the tangent vector of the circle at the point of attachment.  If each of the line segments is vertical (i.e., parallel to the z-axis), we would obtain a cylinder.  However, we want to rotate the line segments as they go around the circle so that one trip around the circle puts a half-twist (Pi radians) in the segments.  The resulting surface is called a Moebius band.  You can easliy build one with a piece of paper by cutting a long rectangluar strip from the paper, giving it a half-twist, and glueing the ends together.

We want to describe the resulting surface parametrically so that we can graph it.  A natural way to do this is to realize that any point on the Moebius band can be described by moving along the circle of radius 4 to a specified location, then traveling out along the line segment through that location to the desired point.  Thus, we can parameterize the Moebius band with two parameters:  one specifying points around the circle of radius 4, the other specifying points on the line segments.  This leads to the following.

rotate[t_] := {{Cos[t], -Sin[t], 0}, {Sin[t], Cos[t], 0}, {0, 0, 1}}
        pt[s_, t_] := {4, 0, 0} + s {Cos[t], 0, Sin[t]}ParametricPlot3D[
      Evaluate[rotate[t] . pt[s, (t/2)]], {s, -1, 1}, {t, 0, 2 π}, 
      PlotPoints -> {3, 21}, Boxed -> False, Axes -> None]

Animations

You have already seen the general structure of an animation in Mathematica  in previous labs.  In this section, several different 3-space animations are discussed.  You will be referred to separate notebooks which contain the already calculated animations for viewing.  This procedure will save you a lot of time since you won't have to wait for the frames of the animation to be calculated.  You can find these animation notebooks in the same folder as this lab.

Remember that animations can require a lot of free RAM, so you may need to shut down the Mathematica  kernel to view the animations smoothly.

Another point to keep in mind when computing an animation is that each frame of the animation is stored in a page description language called PostScript.  If you would like to view the PostScript commands which describe a graphic, you can unformat the graphics cell by toggling the Formatted option in the Attributes option of the Style menu.  For images of surfaces, you can save a lot of memory by converting the description of the graphic from PostScript to Bitmap PICT (use the Graph menu).  A bitmap of a graphic is simply a pixel by pixel description of the image.  As such, it does not resize well, and does not print well since the printer has a much higher resolution than the computer screen.  But for animations, converting to Bitmap PICT usually frees up a lot of badly needed memory to allow the animation to be viewed more smoothly.  Thus the tradeoff is usually a good one.

The Twisted Cubic

When we graphed the twisted cubic curve in an earlier section in this lab, we proposed two animations which would help us to see the details of the curve.  One was to animate the unit tangent and unit normal vectors moving along the curve.  This animation is contained in the notebook named twisted cubic animation.

The other animation discussed earlier for the twisted cubic was simply to rotate the curve around in 3-space.  This type of animation is located in the notebook twisted cubic rotation.

Rotating the Monkey Saddle

The monkey saddle surface was discussed and graphed earlier in this lab.  A notebook named monkey_saddle_rotation.nb contains an animation of the surface rotating in space.

Advanced Topics

Building up Graphics from Primitives

Mathematica  has several graphics primitives which allow you to plot points, draw lines, draw and shade polygons, etc.  If you look at all the graphics thusfar in this lab, they consist primarily of points, lines, and shaded polygons.   By using high-level Mathematica  functions like Plot3D and ParametricPlot3D, we never have to worry about the details of the construction of the graphics at, for example, the shaded polygon level.

If you need to build up a graphic image using primitive graphic objects, you can do so as illustrated in the following example.

Show[
      Graphics3D[{
          Polygon[{{- .5, -1, -1}, {- .5, 1, -1}, {
                -2, 1, 1}}], 
          Polygon[{{ .5, -1, -1}, { .5, 1, -1}, {2, 1, 1}}], 
          Polygon[{{-2, - .75, 0}, {-2, 1, 0}, {2, 1, 0}, {2, 
                - .75, 0}}], Thickness[ .01], RGBColor[ .1,  .4,  .8], 
          Dashing[{ .02,  .02}], Line[{{0, 0, -1}, {0, 0, 1.25}}], 
          Cuboid[{- .25, - .25, 1.25}, { .25,  .25, 1.75}]}], 
      Boxed -> False, ViewPoint -> { .5, -2.5, 1}]

You can specify colors for the objects, or you can let Mathematica  compute the colors based on the lighting model it uses.  This lighting model is discussed briefly in the next subsection.

An example of an animation of two cubes connected by a bar rotating in 3-space is contained in the auxilary notebook dumbbell animation.  There are some minor difficulties which have to be overcome to do this rotation, so the code in that notebook is more subltle than in the other animations.

Lighting and Color

Mathematica  provides a lot of flexability in the way the colors of the polygons are computed when graphing surfaces.  You can specify the color of the background light (AmbientLight), you can specify the color and locations of point sources of light (LightSources), and you can specify how matted or shiny you want the surface to be (SurfaceColor).  Mathematica  combines all of this information using the Phong lighting model to determine the color of each polygon in a 3-dimensional graphic.  Alternatively, you can turn off the simulated lighting model using Lighting -> False.  You can also specify colors for the polygons using the ColorFunction option.  For details on these specifications, you can consult the reference book Mathematica  (pp. 455 -- 459).  A brief description of the Phong lighting model used by Mathematica  is contained there, but for a complete description, you will need to consult a text on computer graphics.

The Structure of SurfaceGraphics and Graphics3D Objects

You may have noticed that the Plot3D function produces a SufaceGraphics object, whereas the ParametricPlot3D function produces a Graphics3D object.  The reason for the difference is that when the one coordinate is a function of the other two, an efficient description of the polygons which model the surface can be empoyed.  To see how Mathematica  does this, we can look at the actual structure of the output from Plot3D, compressing the result to 2 lines of information using the Short command.

Short[
      InputForm[
        Plot3D[x + y, {x, -1, 1}, {y, -2, 2}, 
          DisplayFunction -> Identity, PlotPoints -> 3]], 2]

So the z values are neatly stored in a PlotPoints by PlotPoints array followed by the list of options.  We can compare this to a Graphics3D description of a surface obatined, for example, using ParametricPlot3D.

Short[
      InputForm[
        ParametricPlot3D[{x, y, x + y}, {x, -1, 1}, {y, -2, 2}, 
          DisplayFunction -> Identity, PlotPoints -> {3, 3}]], 9]

So what we get is a polygon by polygon description of the surface.

Exercises

Exercise 1

Consider the vector-valued function

r(t) = (2 Cos[t], 2 Sin[t], t)

This equation describes a helix currve in 3-space.

(a)  Graph a portion of this curve using an interval for t which contains 0 and which is large enough to illustrate the behaviour of the curve.

(b)  Add unit tangent and unit normal vectors to the curve in part (a) at three or more points along the curve.  [Refer to the twisted cubic rotation or twisted cubic animation notebook if you need help with this.]

Exercise 2

Human DNA molecules are modeled using a double helix.  Each helix of this double helix has diameter approximately 20 Angstroms and pitch approximately 34 Angstroms (an Angstrom is 10^-8 cm).  There are about 2.9*10^8 turns in each helix.  How long would such a helix be if it could be unraveled and pulled taught to achieve its full length?

Exercise 3

A Moebius band can be constructed by taking a long rectangular piece of paper and identifying an opposite pair of edges after a half twist is placed in the paper.

(a)  Graph the surface which is constructed by taking such a rectangular piece of paper and identifying an opposite pair of edges after a full twist (i.e., two half twists) is placed in the paper.  

(b)  Do the same problem as in (a), except use 3 half twists.

Exercise 4

(a)  Graph the equation  4x2 + 4y2 - 9z2 = 1  using cylindrical coordinatess.

(b)  The catenoid is the surface obtained by rotating the catenary given by x = Cosh[z] in the xz-plane about the z-axis.  Graph this catenoid.

Exercise 5

Mathematica  defaults to placing axes along three sides of the bounding box which contains a graphic in 3-space.  Override this default, and use the Line graphics primitive to place the x-, y- and z-axes where they belong -- through the origin in 3-space.  In order to do this, choose any graphic of a surface you like, and place the axes accordingly.

Exercise 6

The following self-intersecting surface was obtained by making all the horizontal cross-sections pieces of limacons:  the lower ones from limacons without loops; one of the middle ones from a cardiod; and the top ones from limacons with loops.  Try to recreate such a folded surface using ParametricPlot3D.

Exercise 7

This exercise deals with the function given by

f (x, y) = ((x*(y^2))/((x^2) + (y^2)))

(a)  Use Mathematica  to show that the limit of f(x,y) as (x,y) approaches the origin along any straight line through the origin is 0.  

(b)  Use Mathematica  to show that the limit of f(x,y) as (x,y) approaches the origin along the parabola given by x=y^2 is 1/2.  Since this value disagrees with the values of limits obtained in part (a), conclude that the limit of f(x,y) fails to exist at the origin.

(c)  Graph f(x,y).  Adjust the PlotPoints and ViewPoint options so that an instructive view of the surface is obtained.

(d)  Use the SpinShow function from the Graphics package to produce an animation of the surface rotating in space.

Exercise 8

Earlier in this lab, a cylindrical coordinate plot of the function whose rectangular equation is

f (x, y) = ((x*y)/((x^2) + (y^2)))

was produced.  Study this graph further by rotating it in space.  You can use the SpinShow function from the Graphics package to do this.  However, you will want to take full advantage of the symmetry of the surface to reduce the number of animation frames that need to be computed.  In order to do this, look at the options to SpinShow.  Reduce the range of the spin using the SpinRange option, and reduce the number of frames using the Frames option (but use at least 8 frames).  If you do this right, you should be able to compute all the needed frames to obtain what appears to be a continuous rotation of the surface.  Rotations of this type usually look better without axes and bounding boxes.

Exercise 9

Let  f(x,y) = (x2 - 3xy - y2) exp(-x2 - .1x - y2) .

(a)  Graph the function.

(b) Produce a contour plot of the function.

(c)  Estimate the coordinates of and the values of the relative extrema of the function from the contour plot in (b).  The following technique should be helpful:  click on the contour graphic; while holding down the command key (the one with the apple on it), the approximate coordinates of the cursor appear in the lower left of the window.