GRAPHICAL PROGRAMMING
19 Feb 2003

A Mathematica programming example

We can plot a real function with ease:

Plot[Sin[x],  {x,  0,  Pi}]

Suppose we like to represent a complex function Sin[z] graphically, where z is a complex number.  How can we make a plot of a complex function?

One way to do this is to consider the function as a mapping from a two-dimensional point z=(x,y) to a new point f(z) = (Re f, Im f). Such maps are called conformal maps in mathematics.  In fact, we consider the effect of how lines are mapped.  Specifically, we'll write a function which plot the image of the regular grid

under the mapping of f(z).  For the sine function, we have the image:

Question: which line is associated which line in the original regular grids?  The grids get distorted, but which feature(s) of the original grid is (are) still preserved?

I'll show you the code first, then explain slowly some of the features in the program that you have not seen.  I also demonstrate in class how such a program can be built-up and refined step-by-step.  The code is adapted from the standard Mathematica package "Graphics`ComplexMap`".

Complex map package

BeginPackage["Graphics`ComplexMap`"]

Unprotect[ CartesianMap, Lines ]

CartesianMap::usage = "CartesianMap[f, {x0, x1, (dx)}, {y0, y1, (dy)}] plots the image of the cartesian coordinate lines under the function f. The default values of dx and dy are chosen so that the number of lines is equal to the value of the option Lines."
Lines::usage = "Lines -> {lx, ly} is an option of CartesianMap that gives the number of lines to draw."

Begin["`Private`"]

Needs["Utilities`FilterOptions`"]

$Lines = 15; (* global default *)
Options[CartesianMap] = { Lines :> $Lines }

CartesianMap[ func_, {x0_, x1_, dx_:Automatic}, {y0_, y1_, dy_:Automatic}, opts___?OptionQ ] :=
    Module[ {x, y}, 
    	Picture[ CartesianMap, func[x + I y], {x, x0, x1, dx}, {y, y0, y1, dy}, opts ]
    ] /; NumericQ[x0] && NumericQ[x1] && NumericQ[y0] && NumericQ[y1] && 
     (NumericQ[dx] || dx === Automatic) && (NumericQ[dy] || dy === Automatic)

Picture[ cmd_, e_, {s_, s0_, s1_, ds_}, {t_, t0_, t1_, dt_}, opts___ ] :=
    Module[ {hg, vg, lines, nds = ds, ndt = dt},
        lines = Lines /. {opts} /. Options[cmd];
        If[ Head[lines] =!= List, lines = {lines, lines} ];
        If[ ds === Automatic, nds = N[(s1-s0)/(lines[[1]]-1)] ];
        If[ dt === Automatic, ndt = N[(t1-t0)/(lines[[2]]-1)] ];
        hg = Curves[ e, {s, s0, s1, nds}, {t, t0, t1}, opts ];
        vg = Curves[ e, {t, t0, t1, ndt}, {s, s0, s1}, opts ];
        Show[ Graphics[ Join[hg, vg] ], FilterOptions[Graphics, opts],
              AspectRatio -> Automatic, Axes -> True ] ]

Curves[ xy_, spread_, bounds_, opts___ ] :=
    With[{curves = Table[{Re[xy], Im[xy]}, spread]},
        ParametricPlot[ curves, bounds, DisplayFunction -> Identity, Evaluate[FilterOptions[ParametricPlot, opts]] ][[1]] ]

End[ ]
Protect[ CartesianMap, Lines ]
EndPackage[ ]

Here is an example of how to use the code:

CartesianMap[Exp,  {0,  1},  {0,  1},  
      Lines  ->  5]

I'll show the development process interactively in class.  I may not able to explain all the features of the above program, but wen should be able to code a workable program similar to the above.

The process of writing a program in Mathematica

First try, using straight lines

For simplicity, we consider a version that use straight lines only.  Later on we change it to smooth curves.   We consider a set of horizontal lines and a set of vertical lines, and how these lines are get mapped by a complex function f.  Actually, if we use straight lines in the plot for the image lines, we consider the image points of intersection points in the grids only. The code below gives us only one horizontal line.  Note that a line passing through points {x1,y1}, {x2, y2}, {x3,y3} is represented in Mathematica as Line[{ {x1,y1}, {x2, y2}, ... }].

drawlines[{x0_,  x1_,  dx_},  {y0_,  y1_,  dy_}]  :=  
      Module[ {}, 
 
        y =  y0; 
oneline  =  
          Table[ {x,  y}, {x,  x0,  x1,  
              dx}]; 
oneline
]

I did not declare the local variables like oneline, x, y.  I am just lazy, but not a good practice, for in some cases, one can get unexpected results.

drawlines[{0,  1,  0.2},  {0,  1,  0.2}]

Show[Graphics[Line[%]]]

We modify the drawlines so that we can draw n horizontal lines:

drawlines[{x0_,  x1_,  dx_},  {y0_,  y1_,  dy_}]  :=  
      Module[ {}, 
 
        manylines  =  {}; 
Do[
oneline  =  Table[ {x,  y}, {x,  x0,  x1,  
                dx}]; 
manylines  =  
            Append[manylines,  Line[oneline]], 
{y,  
            y0,  y1,  dy}
]; 
 
        manylines
]

ls  =  drawlines[{0, 1, 0.2},  {0,  1,  0.2}]

Show[Graphics[ls]]

Using functional programming style to avoid the explicit loop, a better version might be.

drawlines[{x0_,  x1_,  dx_},  {y0_,  y1_,  dy_}]  :=  
      Module[ {}, 
 
        manylines  =  
          Table[{x, y},  {y, y0, y1, dy}, {x, x0,  x1, 
              dx}]; 
Line  /@  
          manylines
]

Show[Graphics[drawlines[{0, 1, 0.3}, {0, 1, 0.2}]]]

Note that Graphics[ ] is in the form of Graphics[ { g1, g2, ...}] where g1, g2, ... are graphical elements such as Line[], or Point[], or Polygon[] etc.   Adding vertical line is easy:

drawlines[{x0_,  x1_,  dx_},  {y0_,  y1_,  dy_}]  :=  
      Module[ {}, 
  
hlines  =  
          Line  /@ Table[{x, y},  {y, y0, y1, dy}, {x, x0,  x1, 
                dx}]; 
  
vlines  = 
          Line  /@ Table[{x, y},  {x, x0,  x1, dx}, 
{y, 
                y0, y1, dy}]; 
 
        hlines ~Join~ vlines
]

Show[Graphics[drawlines[{0, 1, 0.2},  {0,  1,  0.2}]],  
      AspectRatio ->  Automatic]

Now we are able to draw grids.  The next thing is to do the mapping. Rather than modify the above code, we do the mapping by a separate function, using the data generated by drawlines[ ].   The problem now is to relate each pair of point (x,y) by its image (Re f(x,y), Im f(x,y)).  The Map[ ] can be used.  Some test

f   /@ { {x1,  y1},  {x2,  y2}}

?Apply

Apply[f, {{x1, y1},  {x2,  y2}}, 2]

This also does not work, we need to convert {x,y} into x + I y.  Also, note that the return value from drawlines is { Line[{ {x1,y2}, ..}], ...}

Map[ #[[1]]  +  I #[[2]]  &,  { 
        Line[{{x1, y1},  {x2,  y2}}]},  {3}]

?Map

Map[{ Re[f[#]],  Im[f[#]]} &,  %41]

We put the above together as a module,

cmap[f_,  {x0_,  x1_,  dx_},  {y0_,  y1_,  dy_}]  :=  
      Module[{},  
ls  =  
          drawlines[{x0, x1,  dx},  {y0,  y1,  
              dy}]; 
z  =  
          Map[ (#[[1]]  +  I #[[2]]) &,  
            ls,  {3}]; 
 
        newlines  = 
          Map[ {Re[f[#]],  Im[f[#]]} &,  
            z,  {3}] 
]

cmap[Identity,  {0,  1,  0.5},  {0, 1,  0.5}]

Show[Graphics[cmap[Sin, {0, 1, 0.2},  {0,  1,  0.2}]],  
      AspectRatio ->  Automatic]

We now have a working program.  We can now write the final function ComplexMap[].  However, you see that our own code is fairly complex.  If we use the built-in function ParametricPlot[], we can do it much simpler.

MATHEMATICA
PROGRAMMING

Plotting the complex function, continued

Previously, we have programmed the ComplexMap with straight line segments.  This picture is faithful only at the intersection points of two perpendicular lines.  The second, and better solution is to use smooth curves, generated by ParametricPlot[].   The mapping then is faithfully represented for all the lines in the pre-image grids.  First, let us see how to draw a single line by ParametricPlot:

y0  =  0.0;
ParametricPlot[{x, y0},  {x,  0, 1}]

How about three lines? try

ParametricPlot[ { {x, 0},  {x, 0.2},  {x,  0.4} },  {x,  0,  
        1}]

What about the image of Sin[z] of these lines?

ls  =  Table[ {x,  y},  {y,  0, 0.4,  0.2}]

This gives three horizontal lines in parametric form (suitable for ParametricPlot to use.

ls  =  
      Table[ {Re[Sin[x + I y]],  Im[Sin[x +  I y]]},  {y,  0,  0.4,  
          0.2}]

ParametricPlot[Evaluate[ls],  {x,  0,  1}]

?With

With[ {curves  =  
          Table[ { Re[Sin[x  +  I y]],  
              Im[ Sin[x  +  I y]]},  {y,  0,  0.4,  
              0.2}]}, 
ParametricPlot[
        curves,  {x,  0,  1}]]

The trick here is that if we use With[ ], we do not need Evaluate[ ].  Now we are ready to define a function to draw the image of a set of parallel lines:

Curves[expr_,  spread_,  bounds_]  :=  
      With[{cvs  =  
            Table[ {Re[expr],  Im[expr]},  
              spread]}, 
ParametricPlot[cvs,  
          bounds]
]

Note that expr must be an expression in terms of x and y, e.g. Sin[x + I y]  (not a pure function).  The spread specifies how many lines to draws, and bounds specify the limits of the lines.  This function works both for horizontal lines as well as vertical lines.  E.g.:

Curves[
      Sin[x  +  I y],  {y,  0.0,  1.0,  0.2},  {x,  0.0,  
        1.0}]

Curves[
      Sin[x  +  I y],  {x,  0.0,  1.0,  0.2},  {y,  0.0,  
        1.0}]

Now we can draw the images of set of horizontal lines and vertical lines separately.  But how can we draw them in the same plot?  One simple way is to say:

hg  =  
      Curves[Sin[x  +  I y],  {y,  0.0,  1.0,  0.2},  {x,  0.0,  
          1.0}]

hv  =  
      Curves[Sin[x  +  I y],  {x,  0.0,  1.0,  0.2},  {y,  0.0,  
          1.0}]

Show[hg, hv,  
      AspectRatio ->  Automatic]

However, you get three plots, one for horizontal lines, one for vertical lines, and the final combined lines.   We need to suppress the display of the first two plots.   Before doing this, let us see what is the FullForm of the graphics data structure hg or vg.

FullForm[hg]

We see hg is   Graphics[{ {{Line[{{x,y}, ...}]}}, {{Line[{{..}}]}} }, {opt->r}].
To turn off the display of the graphics, but retaining the data structure, we use an option in plot function DisplayFunction -> Identity:

Curves[expr_,  spread_,  bounds_]  :=  
      With[{cvs  =  
            Table[ {Re[expr],  Im[expr]},  
              spread]}, 
ParametricPlot[cvs,  bounds,  
          DisplayFunction  ->  
            Identity]
]

Now, we don't see the plot, but generate data only.

hg  =  
      Curves[Sin[x  +  I y],  {y,  0.0,  1.0,  0.2},  {x,  0.0,  
          1.0}]

Show[hg]

Show[Graphics[hg[[1]]]]

Now we write a function call Picture to combine the two set of curves:

Picture[
        expr_,  {s_,  s0_,  s1_,  ds_},  {t_,  t0_,  t1_,  
          dt_}]  :=  
      Module[{ hg,  vg}, 
vg  =  
          Curves[expr,  {s, s0, s1, ds},  {t, t0, 
              t1}]; 
hg  =  
          Curves[expr,  {t,  t0,  t1,  ds},  {s,  s0,  
              s1}]; 
Show[ 
          Graphics[ vg[[1]] ~Join~ hg[[1]]],  
          AspectRatio  ->  Automatic,  
          Axes  ->  True]
]

Picture[ 
      Sin[x  +  I y],  {x,  0,  1,  0.2},  {y,  0,  1,  
        0.2}]

Now we define a user interface so that pure function in complex variable z = x + I y can be specified:

CartesianMap[
        func_,  {x0_,  x1_,  dx_},  {y0_,  y1_,  dy_}]  :=  
      Module[{x, y}, 
Picture[
          func[x  +  I y],  {x,  x0,  x1,  dx},  {y,  y0,  y1,  
            dy}]
]

Now we can say:

CartesianMap[ Sin,  {0,  1,  0.2},  {0,  1,  0.2}]

Comparing to the version in the Graphics`ComplexMap` package, we still have a number of features missing.  For example, the increment dx, and dy can be omitted, if so, a globally define $Lines will be used to decided how many lines to draw.  In addition, the package version can add options in the form

CartesianMap[..., opt  ->  value]

where opt and value are these allowed options for Plot or general graphics.   We'll not go into that details.  But show you the last feature: the package feature.   The main advantage of wrapping your code into a package is that the user can only access the public symbols and not the private symbols or definitions.   Also the package creates separate name space to void name conflict (between two set of programs).  Here is the stripped down outline of a package:

\( (*\ Summary : \ 
        This\ package\ is\ such\ and\ such\ a\ package\  ... \ *) 
\(\(BeginPackage["Graphics`ComplexMap`"]\)\(
\)
    \)\)
Unprotect[CartesianMap]

    
\(\(\(CartesianMap::usage\  = "CartesianMap[...] do such and such..."\)\(
\)
    \)\)
Begin["`Private`"]

    
\(	CartesianMap[ ... ]\  := 
		Module[{}, \ \(\(...\)\(.\)\)]\)
\(	Picture[ .. ]\  := \  ... \)
\(	\(\(Curves[ .. ]\  := \  ... \)\(
\)
    \)\)
End[]
Protect[CartesianMap]
EndPackage[]

The above text (code) will be saved as a .m file, you load in the package with <<Graphics`ComplexMap`.   Then only the function Cartesian map is available to the user.

SUMMARY OF MATHEMATICA

(1): Computation can be performed in user-specified precision.
(2): Ability to so symbolic computation with math formulas.
(3): Rich programming environment, support various programming styles.

Hope that you'll find Mathematica a very useful tool in your work.

Reading: Sec. 2.6.

If you have time, try to cover Part 1. A practical introduction, and Part 2. Principles, of the Mathematica book.