Principia Matematicae

Dr. C. Chryssomalakos, ICN - UNAM

This is a set of notes on some general aspects of Mathematica. Some of the material is scattered in the references cited in the bibliography, mostly Refs. [2], [3], [5] ---the rest are things unearthed while struggling with the subject.

SetDirectory["E:\chryss\math\workshop"];

Lecture 4: Graphics

Basic Plotting

2-D

The basic 2-D plotting command is Plot

a = 2;

Plot[Sin[a x^2], {x, 0, 2  Pi}]

You can add options to Plot and give a name to its output

\(graph1 = 
      Plot[
Sin[a\ x^2], 
{x, 0, 
          4  Pi}, 
AspectRatio -> 
          Automatic, 
PlotPoints -> 
          200, 
AxesLabel -> {"x", ""}, 
PlotLabel -> \ 
          StringJoin["a=", 
            ToString[a]], 
PlotRange -> {{Pi, 
              2  Pi}, {\(-1.1\), 
              1.1}}, 
AxesOrigin -> {1.5  Pi, 
            0}, 
ImageSize -> \ 
          400, 
PlotStyle -> {
Thickness[ .005], 
Hue[ .05]
}
]\)

Mathematica can handle poles in the arguments of Plot

graph2 = 
      Plot[
Zeta[x], 
{x, -2, 
          6}, 
PlotStyle -> {
Thickness[ .003], 
Hue[ .7]
}
]

Plot, and the other plotting commands, evaluate their argument only after giving numerical values to the plot variable(s). Evaluation can be forced with Evaluate (recommendable if simplification occurs). Compare the two examples below

{
Timing[
Plot[
Integrate[Sin[x^2]^2, {x, 0, y}], {y, 0, 
            1}, 
PlotPoints -> 
            200
]
], 
Timing[
Plot[
Evaluate[
Integrate[Sin[x^2]^2, {x, 0, y}]
], {y, 
            0, 1}, 
PlotPoints -> 
            200
]
]
}

The first Plot evaluates the integral for each numerical value of y (200 times...) and then plots the corresponding point. The second Plot only evaluates the integral once for general y

Integrate[Sin[x^2]^2, {x, 0, y}]

and then plots the resulting function of y.

The actual input that stores the information about a plot can be consulted with InputForm (this can be an enormous list)

InputForm[graph2]

The head of the output is Graphics. The first argument is a list containing 
 Thickness and Hue directives and a Line primitive (itself containing a list of the points it passes through). The second argument is a list of options - in graph2 Mathematica used default values for them.

The information about graph2 can be read (and the plot redrawn) with Show - new options can be specified

Show[graph2, AspectRatio -> Automatic]

Plots can be combined with Show

Show[graph1, graph2]

Mathematica adjusts the scale and range automatically, sometimes missing important parts of the plot (graph2 is almost invisible above) -- this can be fixed manually with PlotRange

Show[graph1, graph2, 
      PlotRange -> {{-2, 6}, {-1, 
            4}}
]

Show does not generate any new points, it only reads and displays plots generated with plotting commands. Using Show, options can be changed repeatedly without having to generate the plot from scratch each time (changing ViewPoint for example).

3-D

The basic 3-D plotting command is Plot3D

\(Plot3D[
Sin[5  Sin[x] y] y^2 
         Cos[ .1  y^1.5], 
{x, 0, 
        Pi}, 
{y, 0, 
        2  Pi}, 
PlotPoints -> {100, 
          100}, 
ViewPoint -> {1,  .8, 
          2}, 
AxesLabel -> {"x", "y", "f"}, 
ViewVertical -> {0, 1, 
          0}, 
ImageSize -> \ 
        350
]\)

More versatile, but slower, is ParametricPlot3D. It can produce a curve or a surface depending on the number of parameters used. In the example below phi, theta parametrize the surface of a torus of radii R, rho

R = 1; rho =  .4;

ParametricPlot3D[
{
(R + 
              rho  Sin[phi]) Cos[theta],  (R + rho  Sin[phi]) 
          Sin[theta],  
rho  Cos[
            phi], 
RGBColor[
Cos[
              phi]^2, 
Sin[
              theta + phi]^2, 
Cos[
              theta]^2
]
}, 
{phi, 0, 2  Pi}, 
{theta, 0, 
        2  Pi}, 
Lighting -> 
        False, 
PlotPoints -> {40, 
          120}, 
ImageSize -> 
        400, 
Boxed -> 
        False, 
Axes -> 
        False
]

The fourth entry in the list to be plotted specifies directives to be used. When plotting a curve these could be Thickness, Hue, etc., and their arguments can depend on the parameter, thus changing the thickness or color of the curve along its length . For a surface, as above, the directive is usually a color specification (the arguments of RGBColor have to be between 0 and 1). Colors only appear in the output when Lighting is turned off. With Lighting on, the objects of the plot are illuminated by point sources, with adjustable color, intensity and direction (see the option LightSources in later examples). The reflective characteristics of the objects plotted can also be adjusted (see SurfaceColor).

Here is another example: put theta = phi in the torus above to get the torus "diagonal"

ParametricPlot3D[
{
(R + 
              rho  Sin[phi]) Cos[phi],  (R + rho  Sin[phi]) 
          Sin[phi],  
rho  Cos[
            phi], 
{
RGBColor[
Cos[phi]^2, 
Sin[
                2  phi]^2, 
Cos[
                phi]^2
], 
Dashing[{ .01,  .02}], 
Thickness[ .003 +  .001 
                 Cos[phi]]
}
}, 
{phi, 0, 2  Pi}, 
Lighting -> 
        False, 
PlotPoints -> 
        40, 
ImageSize -> 
        400, 
Boxed -> 
        False, 
Axes -> 
        False
]

The Structure of Graphics

Graphics have head and arguments. The head can be Graphics, DensityGraphics, SurfaceGraphics, Graphics3D, GraphicsArray etc., depending on the type of the plot. The arguments are divided into graphics primitives and options. Graphics primitives appear first, and can be graphics objects like Point, Line and 
Polygon or graphics directives, like Thickness, RGBColor and SurfaceColor, which affect the way primitives are displayed. This can also be achieved with options, but the latter apply globally, while directives can be restricted to individual primitives. 

Graphics can be generated "manually", i.e., without use of any of the built-in plotting commands of Mathematica. Notice the scope of the directives in the following example

\(graph3 = 
      Graphics[
{
PointSize[ .05], 
Hue[ .6], 
{
PointSize[ .04], 
Hue[ .3], 
Point[{ .15,  .2}], 
Hue[ .9], 
Point[{ .4,  .3}]
}, 
Point[{ .7,  .5}], 
{
PointSize[ .08], 
RGBColor[1, 0, 
              0], 
Point[{2.3, 
                1}]
}, 
{
Thickness[ .009], 
Line[{{0, 0}, { .3, 
                  1}, {2.3, 
                  1.2}}]
}
}, 
AspectRatio -> 
          Automatic, 
Axes -> 
          True, 
AxesLabel -> \ {"horizontal", "vertical"}
]\)

Show[graph3]

Each directive affects all primitives (if relevant) between its position and the closest closing curly bracket, with "deeper" and later directives having priority over "higher" and earlier ones.

graph4 = 
      Graphics3D[
{
Hue[ .3, 
            1,  .4], 
Polygon[
{
{0, 0, 0}, 
{1, 0, 
                0}, 
{1, 1, 0}, 
{0, 1, 
                0}
}
], 
{
RGBColor[1, 0, 
              0], 
Polygon[
{
{0, 0, 0}, 
{0, 0, 
                  1}, 
{0, 1, 
                  0}
}
]
}
}
]

Here no options were specified so Mathematica will use default values.

\(Show[
graph4, 
Lighting -> 
        False, 
Axes -> 
        True, 
AxesLabel -> {"x", "y", "z"}
]\)

Manipulating Graphics

Needs["Graphics`Shapes`"]

Inserting Directives

Can one change PointSize in ListPlot? Check the options available

Options[ListPlot]

PlotStyle is available, so PlotStyle->{PointSize[]} will work (but only globally). We show other ways of doing this (plus more) below. One can always choose to do things manually, starting with InputForm

graph = ListPlot[Table[N[n^Sin[n/50 Pi]], {n, 1, 50}]]

InputForm[graph]

A cheap way to change PointSize is to sacrifice the first point of the graph

graph2 = graph /. Point[{1, 1. }] -> PointSize[0.05]

Now the first entry in the list of primitives in graph2 is a PointSize directive

InputForm[graph2]

Show[graph2]

A better way is to write a function that adds a global directive dir in the beginning of the list of primitives contained in graph, keeping the options

prependToGraph[graph_, dir_] := 
      Graphics[Prepend[graph[[1]], dir], graph[[2]]]

We now add PointSize[0.02] at the beginning of the list of primitives in graph

graph3 = prependToGraph[graph, PointSize[0.02]] // Show

Add a global color directive

graph4 = 
      prependToGraph[
graph3, 
          RGBColor[1, 0,  .6]
] // Show

Try to change again PointSize

graph5 = prependToGraph[graph4, PointSize[ .05]] // Show

Since directives are added at the beginning, PointSize can no longer be changed with prependToGraph (there is another PointSize folowing it). It can still be done though with a substitution

graph4 /. PointSize[_] -> PointSize[0.04] // Show

This is the most powerful trick to get directives inserted "locally" into a graph. One can change for example, the color and size of the points in graph individually, depending on their coordinates, or flip the graph, etc., etc.,...

graph /. 
        Point[{x_, 
              y_}] :>  
{
Hue[Cos[ .03  x]^2], 
            PointSize[ .001  y +  .005], 
Point[{y, 
                x}]
} // Show

The package Shapes includes commands for generating various 3-D objects and moving them around

Graphics3D[
TranslateShape[
Torus[
            3, 1, 10, 10], 
{0, 0, 
            0}
]
] // Show

Here we map a pure function that translates a torus by its argument to a list of twenty random vectors (a new plot is produced every time the command is run)

Graphics3D[
TranslateShape[
Torus[ .3,  .1, 10, 
                10], 
#
] & /@ 
(10 {Random[], Random[], Random[]} &) /@ 
            Range[1, 20]
] // Show

We use the above construction to replace points by toruses (each with its own Hue directive) in graph -- the plot is embedded in the x-z plane.

Show[
Graphics3D[
{
Hue[#[[1, 2]]/
                  60], 
TranslateShape[
Torus[
 .5  #[[1, 
                        1]], 
 .2  #[[1, 2]] + 
                    3, 
10, 
10
], 
{20  #[[1, 1]], 0, 
                  20  #[[1, 
                        2]]}]
} &
 /@ 
          graph[[1]]], 
ImageSize -> 
        400, 
Lighting -> 
        False, 
Axes -> 
        True
]

Examples

From Xfig to Riemman

Needs["Graphics`Shapes`"]

The command Import reads data from a file. nontriv.eps contains a two-dimensional image generated with xfig

psplot = Import["nontriv.eps", "EPS"]

The tricky point is that Import converts the data to a Mathematica expression that produces the same visual output (!). In the example above, the extension .eps instructed Import that the image in the file was in Postscript. psplot now contains Mathematica commands that produce equivalent output

Show[psplot, AspectRatio -> Automatic]

InputForm[psplot]

Replace the 1 with a Line command and embed the graph in 3-D

xT = 3.5;  zT =  .1;  R = 250;

psplot1 = 
psplot /. 
{
AbsoluteThickness[x_] -> AbsoluteThickness[2  x], 
            Text[__] -> 
              Line[{
{450.4 + xT, 65.0 + zT}, {456.4 + xT, 
                    65.0 + zT}, {453.4 + xT, 65.0 + zT}, {453.4 + xT, 
                    80.5 + zT}, {450.4 + xT, 
                    77.0 + zT}
}
]
} /. 
{
{U_Real, 
              V_Real} -> 
{U - 230, V - 70, 0.0}, 
          Line[{{0, 0}, {0, 
                  0}}] -> 
RGBColor[ .1,  .1,  .1]
}

Change the head to Graphics3D (so that 3-D points are recognized)

\(Show[
Apply[Graphics3D, psplot1], 
      AxesLabel -> {"x", "y", "z"}, 
Axes -> True, 
ViewPoint -> {\(-100\), 200, 
          300}, 
ImageSize -> \ 
        350
]\)

This might look the same as before, but the object now lives in 3-D (try changing the ViewPoint).

Project the points in psplot1 stereographically on the Riemann sphere

psplot2 = 
      psplot1 /. 
{
{U_Real, V_Real, 
              0. } -> {U (R + 
                      R (R^2 - (U^2 + V^2))/(R^2 + (U^2 + 
                                  V^2)))/R, 
              V (R + 
                      R (R^2 - (U^2 + V^2))/(R^2 + (U^2 + 
                                  V^2)))/R, 
              R (R^2 - (U^2 + V^2))/(R^2 + (U^2 + 
                          V^2))}, 
RGBColor[__] ->  RGBColor[1. , 0, 0]
}

Show[Graphics3D @@ psplot2]

Draw a cut sphere

sfaira2 = 
ParametricPlot3D[
{
 .99  R Sin[theta] Cos[phi],  .99  R Sin[theta] 
            Sin[phi],  .99  
            R Cos[theta], 
{Hue[ .5 (1 + 
                    Tanh[ .2  theta]), 1,  .95], 
            EdgeForm[]}
}, 
{phi, -Pi/
            4, Pi + Pi/4}, 
{theta, 0, 
          Pi}, 
PlotPoints -> {80, 
            50}, 
Lighting -> 
          False
]

Put everything together

\(\(nonTrivProof = 
Show[
{sfaira2, 
Apply[Graphics3D, psplot2], 
            Apply[Graphics3D, 
              psplot1]
}, 
Boxed -> 
            False, 
Axes -> 
            False, 
ViewPoint -> {200, \(-900\), 
              1000}, AspectRatio -> 
            Automatic, 
ImageSize -> 
            400, 
AxesLabel -> {"x", "y", "z"}, 
Lighting -> 
            False
];\)\)

Export["E:\chryss\math\workshop\stereoProof.eps", 
      nonTrivProof, "EPS"]

Animations

A list of plots can be animated by selecting it and then choosing Animate Selected Graphics from the Cell tab in the menu (or just double-clicking any of the selected frames - click again once to stop the animation).

Off[General::spell, General::spell1]

Needs["Graphics`Colors`"]

Needs["Graphics`Shapes`"]

The following example produces a list of plots of a solution to Schroedinger's equation for a charged particle on a cylinder, in the presence of a radial magnetic field. 

First define the wavefunction, its modulus and its phase

b = 1.8;

psi[n_][phi_, z_, t_] := 
      E^(-I n^2 b^2 Sin[t] - I n b z Sin[t] - 
            1/2 (z - n b Cos[t] + n b)^2  + I n phi)

wave[phi_, z_, 
        t_] := psi[1][phi, z, t] + psi[-1][phi, z, t]

abswave[phi_, z_, t_] := Abs[wave[phi, z, t]]

argwave[phi_, z_, t_] := Arg[wave[phi, z, t]]

Plot the wavefunction, codifying its phase in the color of the surface (see the examples in Hue in the HelpBrowser) - the modulus is proportional to the radial distance of the surface from that of the cylinder, it also affects the brightness (third entry in Hue)

sgra = 
Table[
ParametricPlot3D[
{
x, 
(1 +  .7 
                     abswave[phi, x, t]) 
              Cos[phi], (1 +  .7  abswave[phi, x, t]) 
              Sin[phi], {Hue[
argwave[phi, x, t]/2/
                  Pi, 
1, 
Min[
                  1,  .45 + 1.3  abswave[phi, x, t]]
], 
              EdgeForm[]
}
}, 
{phi, 0, 2  Pi}, 
{x, -4.5, 
            4.5}, 
PlotRange -> 
{{-4.5, 4.5}, {-3, 3}, {-3, 
                3}}, 
PlotPoints -> {50, 
              50}, 
Lighting -> 
            False, 
Boxed -> 
            False, 
Axes -> 
            False
], 
{t, Pi/12, 
          2  Pi - Pi/12, Pi/12}
]

More animations like this can be seen in http://www.nuclecu.unam.mx/~chryss.

Holes

Surfaces with holes can be produced by substituting each Polygon by four smaller ones, leaving a hole in the middle

eps =  .25;

sgra2 = 
sgra[[1]] /. 
Polygon[{p1_, p2_, p3_, 
              p4_}] :> {
EdgeForm[
              Thickness[ .001]], 
FaceForm[
Hue[ .27, 1, 
                1], 
Hue[ .8,  .7,  .4]
], 
Polygon[{p1, p2, 
                p2 + eps (p4 - p2), p1 + eps (p3 - p1)}], 
            Polygon[{p2, p3, p3 + eps (p1 - p3), 
                p2 + eps (p4 - p2)}], 
            Polygon[{p3, p4, p4 + eps (p2 - p4), 
                p3 + eps (p1 - p3)}], 
            Polygon[{p4, p1, p1 + eps (p3 - p1), 
                p4 + eps (p2 - p4)}]
}

pyramidPlot

Polygons can be turned into little pyramids to produce pineapple-looking spheres.

pyramidPlot[graph_, eps1_] := 
      graph /. Polygon[{p1_, p2_, p3_, 
              p4_}] -> {
EdgeForm[Thickness[ .001]], 
            Polygon[
{p1, p2, 
                eps1 (p1 + p2 + p3 + 
                      p4)}
], 
Polygon[
{p2, p3, 
                eps1 (p1 + p2 + p3 + 
                      p4)}
], 
Polygon[
{p3, p4, 
                eps1 (p1 + p2 + p3 + 
                      p4)}
], 
Polygon[
{p4, p1, 
                eps1 (p1 + p2 + p3 + 
                      p4)}
]
}

pineapple = 
Show[
Graphics3D[
pyramidPlot[
            Sphere[2, 18, 
              10],  .28]
], 
DisplayFunction -> $DisplayFunction, 
PlotRange -> 
          All, 
AspectRatio -> 
          Automatic, 
ViewPoint -> {5, 7, 
            2.5}, 
ImageSize -> 
          350, 
LightSources -> {
{{5, 8, 3}, 
              RGBColor[ .1, 0, 1]}, 
{{4, 4, 4}, 
              RGBColor[ .8, 0,  .1]}, 
{{-4, 3, 4}, 
              RGBColor[0,  .7, 
                0]}
}
]

Exercise

What goes wrong at the poles?

Stereoscopia

Plotting an object from two slightly different viewpoints and then superimposing the two images by relaxing the eye muscles produces 3-D perception (see [3])

{g1, 
          g2} = 
Map[
Show[
sgra2, 
ViewPoint -> #, 
ImageSize -> 110, 
              ViewVertical -> {1, 0, 
                  0}, 
DisplayFunction ->  
                Identity, 
Lighting -> True, 
              LightSources -> {
{{5, 8, 3}, 
                    RGBColor[ .1, 0, 1]}, 
{{4, 4, 4}, 
                    RGBColor[ .8, 0,  .1]}, 
{{-4, 3, 
                      4}, RGBColor[0,  .7, 
                      0]}
}
] &, 
{{5, 8, 3}, {5, 7.2, 
              3.5}}
];

DisplayFunction->Identity avoids rendering (i.e., the plotting command only computes the points of the object but does not display it, taking into account lighting, viewpoint etc.).

Show[
GraphicsArray[
{g1, 
          g2}, 
DisplayFunction -> $DisplayFunction, 
        GraphicsSpacing -> - .51
]
]

Curvature Color Coding

The four points defining a typical polygon on a surface do not lie on exactly the same plane. A measure of this non-planarity is given by the norm  of the  cross product of the normals to the two triangles the quadrilateral can be divided into.

Define a trefoil knot as a variation on the torus

R = 1 +  .4  Cos[1.5  theta]; 
    rho =  .25 +  .05  Cos[12  theta];

bubble = 
      ParametricPlot3D[
{
(R + 
                rho  Sin[phi]) Cos[theta],  (R + rho  Sin[phi]) 
            Sin[theta],  
rho  Cos[phi] +  .4 
               Sin[1.5 theta]
}, 
{phi, 
          0, 2  Pi}, 
{theta, 0, 
          4  Pi}, 
Lighting -> 
          True, 
PlotPoints -> {30, 
            400}, 
ImageSize -> 
          400, 
Boxed -> 
          False, 
Axes -> 
          False, 
ViewPoint -> {6, 2, 
            5}
]

Calculate the non-planarity of a quadrilateral. Divide by the tenth power of the length of its side to get an approximate invariant under PlotPoints rescaling (why tenth?). Raise to a small power for range compression.

curv[p1_, p2_, p3_, 
          p4_] := 
1.2 (
(Plus @@ (Cross[
Cross[p2 - p1, p3 - p1], 
                          Cross[p3 - p1, 
                            p4 - p1]
]^2))/
(Plus @@ ((p2 - 
                              p1)^2))^5
)^ .3;

Substitute each polygon in bubble by a colored one

bubble1 = 
      bubble /. 
Polygon[{p1_, p2_, p3_, 
              p4_}] :> {
EdgeForm[
              Thickness[ .001]], 
Hue[
              curv[p1, p2, p3, p4]], 
Polygon[{p1, p2, p3, 
                p4}]
}

cPlot = 
      Show[
bubble1, 
DisplayFunction ->  $DisplayFunction, 
ViewPoint -> {6, 2, 
            5}, 
ImageSize ->  
          460, 
Background ->  
          Automatic, 
Lighting -> 
          False, 
LightSources -> {
{{-3, 4, 3}, 
              RGBColor[1, 0, 1]}, 
{{3, 4, 2.5}, 
              RGBColor[0, 0, 0]}, 
{{-4, 3, 4}, 
              RGBColor[0, 0, 
                0]}
}
]

PlotPoints->{40,600} gives better looking results (but it takes time...).

Bibliography

[1] S. Wolfram, "The Mathematica Book", Cambridge University Press (1999), fourth edition.

[2] J. W. Gray, "Mastering Mathematica", Academic Press (1997).

[3] H. Ruskeepaa, "Mathematica Navigator", Academic Press (1998).

[4] R. Maeder, "Programming in Mathematica", Addison Wesley (1997), third edition.

[5] T. W. Gray, J. Glynn, "Exploring Mathematics with Mathematica", Addison Wesley (1991)

[6] A. Gray, "Modern Differential Geometry of Curves and Surfaces with Mathematica", CRC Press (1998).