10 Interactive Programming
Programming in the interpreter is easy. So is the use of Axiom's graphics facility. Both are rather flexible and allow you to use them for many interesting applications. However, both require learning some basic ideas and skills.
All graphics examples in the gallery section are either produced directly by interactive commands or by interpreter programs. Four of these programs are introduced here. By the end of this chapter you will know enough about graphics and programming in the interpreter to not only understand all these examples, but to tackle interesting and difficult problems on your own. The appendix on graphics lists all the remaining commands and programs used to create these images.
10.1 Drawing Ribbons Interactively
We begin our discussion of interactive graphics with the creation of a useful facility: plotting ribbons of two-graphs in three-space. Suppose you want to draw the two-dimensional graphs of functions all over some fixed range of . One approach is to create a two-dimensional graph for each one, then superpose one on top of the other. What you will more than likely get is a jumbled mess. Even if you make each function a different color, the result is likely to be confusing.
A better approach is to display each of the in three ribbon dimensions as a ``ribbon'' of some appropriate width along the -direction, laying down each ribbon next to the previous one. A ribbon is simply a function of and depending only on .
We illustrate this for defined as simple powers of for ranging between and .
Draw the ribbon for .
Now that was easy! What you get is a ``wire-mesh'' rendition of the ribbon. That's fine for now. Notice that the mesh-size is small in both the and the directions. Axiom normally computes points in both these directions. This is unnecessary. One step is all we need in the -direction. To have Axiom economize on -points, we re-draw the ribbon with option .
Re-draw the ribbon, but with option so that only step is computed in the direction.
The operation has created a viewport, that is, a graphics window on your screen. We assigned the viewport to and now we manipulate its contents.
Graphs are objects, like numbers and algebraic expressions. You may want to do some experimenting with graphs. For example, say
showRegion(vp, "on")
to put a bounding box around the ribbon. Try it! Issue to rotate the figure longitudinal degrees and latitudinal degrees.
Here is a different rotation. This turns the graph so you can view it along the -axis.
There are many other things you can do. In fact, most everything you can do interactively using the three-dimensional control panel (such as translating, zooming, resizing, coloring, perspective and lighting selections) can also be done directly by operations (see Chapter ugGraph for more details).
When you are done experimenting, say to restore the picture to its original position and settings.
Let's add another ribbon to our picture---one for . Since ranges from to for the first ribbon, now let range from to . This puts the second ribbon next to the first one.
How do you add a second ribbon to the viewport? One method is to extract the ``space'' component from the viewport using the operation subspacesubspaceThreeDimensionalViewport. You can think of the space component as the object inside the window (here, the ribbon). Let's call it . To add the second ribbon, you draw the second ribbon using the option .
Extract the space component of .
Add the ribbon for alongside that for .
Unless you moved the original viewport, the new viewport covers the old one. You might want to check that the old object is still there by moving the top window.
Let's show quadrilateral polygon outlines on the ribbons and then enclose the ribbons in a box.
Show quadrilateral polygon outlines.
Enclose the ribbons in a box.
This process has become tedious! If we had to add two or three more ribbons, we would have to repeat the above steps several more times. It is time to write an interpreter program to help us take care of the details.
10.2 A Ribbon Program
The above approach creates a new viewport for each additional ribbon. A better approach is to build one object composed of all ribbons before creating a viewport. To do this, use makeObject rather than draw. The operations have similar formats, but draw returns a viewport and makeObject returns a space object.
We now create a function drawRibbons of two arguments: , a list of formulas for the ribbons you want to draw, and , the range over which you want them drawn. Using this function, you can just say
drawRibbons([x**2, x**3], x=-1..1)
to do all of the work required in the last section. Here is the drawRibbons program. Invoke your favorite editor and create a file called ribbon.input containing the following program.
drawRibbons(flist, xrange) ==}{}
sp := createThreeSpace() Create empty space $sp$.
y0 := 0 The initial ribbon position.
for f in flist repeat For each function $f$,
makeObject(f, xrange, y=y0..y0+1, create and add a ribbon
space==sp, var2Steps == 1) for $f$ to the space $sp$.
y0 := y0 + 1 The next ribbon position.
vp := makeViewport3D(sp, "Ribbons") Create viewport.
drawStyle(vp, "shade") Select shading style.
outlineRender(vp, "on") Show polygon outlines.
showRegion(vp,"on") Enclose in a box.
n := # flist The number of ribbons
zoom(vp,n,1,n) Zoom in x- and z-directions.
rotate(vp,0,75) Change the angle of view.
vp Return the viewport.
Here are some remarks on the syntax used in the drawRibbons function (consult Chapter ugUser for more details). Unlike most other programming languages which use semicolons, parentheses, or begin--end brackets to delineate the structure of programs, the structure of an Axiom program is determined by indentation. The first line of the function definition always begins in column 1. All other lines of the function are indented with respect to the first line and form a pile (see ugLangBlocks ).
The definition of drawRibbons consists of a pile of expressions to be executed one after another. Each expression of the pile is indented at the same level. Lines 4-7 designate one single expression: since lines 5-7 are indented with respect to the others, these lines are treated as a continuation of line 4. Also since lines 5 and 7 have the same indentation level, these lines designate a pile within the outer pile.
The last line of a pile usually gives the value returned by the pile. Here it is also the value returned by the function. Axiom knows this is the last line of the function because it is the last line of the file. In other cases, a new expression beginning in column one signals the end of a function.
The line drawStyle(vp,"shade") is given after the viewport has been created to select the draw style. We have also used the zoomzoomThreeDimensionalViewport option. Without the zoom, the viewport region would be scaled equally in all three coordinate directions.
Let's try the function drawRibbons. First you must read the file to give Axiom the function definition.
Read the input file.
Draw ribbons for for
10.3 Coloring and Positioning Ribbons
Before leaving the ribbon example, we make two improvements. Normally, the color given to each point in the space is a function of its height within a bounding box. The points at the bottom of the box are red, those at the top are purple.
To change the normal coloring, you can give an option . When Axiom goes about displaying the data, it determines the range of colors used for all points within the box. Axiom then distributes these numbers uniformly over the number of hues. Here we use the simple color function for the -th ribbon.
Also, we add an argument so you can give the range of occupied by the ribbons. For example, if the is given as and there are ribbons to be displayed, each ribbon would have width and would appear in the range .
Refer to lines 4-9. Line 4 assigns to the variable part of the (after all, it need not be ). Suppose that is given as where and have numerical values. Then line 5 assigns the value of to the variable . Line 6 computes the width of the ribbon by dividing the difference of and by the number, , of ribbons. The result is assigned to the variable . Note that in the for-loop in line 7, we are iterating in parallel; it is not a nested loop.
drawRibbons(flist, xrange, yrange) ==}{}
sp := createThreeSpace() Create empty space $sp$.
num := # flist The number of ribbons.
yVar := variable yrange The ribbon variable.
y0:Float := lo segment yrange The first ribbon coordinate.
width:Float := (hi segment yrange - y0)/num The width of a ribbon.
for f in flist for color in 1..num repeat For each function $f$,
makeObject(f, xrange, yVar = y0..y0+width, create and add ribbon to
var2Steps == 1, colorFunction == (x,y) +-> color, _
space == sp) $sp$ of a different color.
y0 := y0 + width The next ribbon coordinate.
vp := makeViewport3D(sp, "Ribbons") Create viewport.
drawStyle(vp, "shade") Select shading style.
outlineRender(vp, "on") Show polygon outlines.
showRegion(vp, "on") Enclose in a box.
vp Return the viewport.
10.4 Points, Lines, and Curves
What you have seen so far is a high-level program using the graphics facility. We now turn to the more basic notions of points, lines, and curves in three-dimensional graphs. These facilities use small floats (objects of type DoubleFloat) for data. Let us first give names to the small float values and .
The small float 0.
The small float 1.
The @ sign means ``of the type.'' Thus is of the type DoubleFloat. You can also say .
Points can have four small float components: coordinates and an optional color. A ``curve'' is simply a list of points connected by straight line segments.
Create the point with color zero, that is, the lowest color on the color map.
Create the point with color zero.
Create the curve (well, here, a line) from to .
We make this line segment into an arrow by adding an arrowhead. The arrowhead extends to, say, on the left, and to, say, on the right. To describe an arrow, you tell Axiom to draw the two curves and We also decide through experimentation on values for , the ratio of the size of the arrowhead to the stem of the arrow, and , the angle between the arrowhead and the arrow.
Invoke your favorite editor and create an input file called arrows.input. This input file first defines the values of
and , then defines the function makeArrow to draw an arrow from point to .
arrowAngle := %pi-%pi/10.0@DFLOAT The angle of the arrowhead.
arrowScale := 0.2@DFLOAT The size of the arrowhead
relative to the stem.
makeArrow(p1, p2) ==
delta := p2 - p1 The arrow.
len := arrowScale * length delta The length of the arrowhead.
theta := atan(delta.1, delta.2) The angle from the x-axis
c1 := len*cos(theta + arrowAngle) The x-coord of left endpoint
s1 := len*sin(theta + arrowAngle) The y-coord of left endpoint
c2 := len*cos(theta - arrowAngle) The x-coord of right endpoint
s2 := len*sin(theta - arrowAngle) The y-coord of right endpoint
z := p2.3*(1 - arrowScale) The z-coord of both endpoints
p3 := point [p2.1 + c1, p2.2 + s1, z, p2.4] The left endpoint of head
p4 := point [p2.1 + c2, p2.2 + s2, z, p2.4] The right endpoint of head
[ [p1, p2, p3], [p2, p4] ] The arrow as a list of curves
Read the file and then create an arrow from the point to the point .
Read the input file defining makeArrow.
Construct the arrow (a list of two curves).
Create an empty object of type .
Add each curve of the arrow to the space .
Create a three-dimensional viewport containing that space.
Here is a better viewing angle.
10.5 A Bouquet of Arrows
Let's draw a ``bouquet'' of arrows. Each arrow is identical. The arrowheads are uniformly placed on a circle parallel to the -plane. Thus the position of each arrow differs only by the angle , , between the arrow and the -axis on the -plane.
Our bouquet is rather special: each arrow has a different color (which won't be evident here, unfortunately). This is arranged by letting the color of each successive arrow be denoted by . In this way, the color of arrows ranges from red to green to violet. Here is a program to draw a bouquet of arrows.
drawBouquet(n,title) ==}{}
angle := 0.0@DFLOAT The initial angle
sp := createThreeSpace() Create empty space $sp$
for i in 0..n-1 repeat For each index i, create:
start := point [0.0@DFLOAT,0.0@DFLOAT,0.0@DFLOAT,angle]
the point at base of arrow;
end := point [cos angle, sin angle, 1.0@DFLOAT, angle]
the point at tip of arrow;
arrow := makeArrow(start,end) the $i$th arrow
for a in makeArrow(start,end) repeat For each arrow component,
curve(sp,a) add the component to $sp$
angle := angle + 2*%pi/n The next angle
makeViewport3D(sp,title) Create the viewport from $sp$
Read the input file.
A bouquet of a dozen arrows.
10.6 Diversion: When Things Go Wrong
Up to now, if you have typed in all the programs exactly as they are in the book, you have encountered no errors. In practice, however, it is easy to make mistakes. Computers are unforgiving: your program must be letter-for-letter correct or you will encounter some error.
One thing that can go wrong is that you can create a syntactically incorrect program. As pointed out in Diversion 1 the meaning of Axiom programs is affected by indentation.
The Axiom parser will ensure that all parentheses, brackets, and braces balance, and that commas and operators appear in the correct context.
A common mistake is to misspell an identifier or operation name. These are generally easy to spot since the interpreter will tell you the name of the operation together with the type and number of arguments which it is trying to find.
Another mistake is to either to omit an argument or to give too many. Again Axiom will notify you of the offending operation.
Indentation makes your programs more readable. However there are several ways to create a syntactically valid program. A most common problem occurs when a line is indented improperly. If this is a first line of a pile then all the other lines will act as an inner pile to the first line. If it is a line of the pile other than the first line Axiom then thinks that this line is a continuation of the previous line. More frequently than not a syntactically correct expression is created. Almost never however will this be a semantically correct. Only when the program is run will an error be discovered.
10.7 Drawing Complex Vector Fields
We now put our arrows to good use drawing complex vector fields. These vector fields give a representation of complex-valued functions of complex variables. Consider a Cartesian coordinate grid of points in the plane, and some complex-valued function defined on this grid. At every point on this grid, compute the value of and call it . Since has both a real and imaginary value for a given grid point, there are four dimensions to plot. What do we do? We represent the values of by arrows planted at each grid point. Each arrow represents the value of in polar coordinates . The length of the arrow is proportional to . Its direction is given by .
The code for drawing vector fields is in the file vectors.input. We discuss its contents from top to bottom.
Before showing you the code, we have two small matters to take care of. First, what if the function has large spikes, say, ones that go off to infinity? We define a variable for this purpose. When exceeds the value of , then the value of is used instead of that for . For convenience, we define a function which uses to ``clip'' the value of .
clipValue : DFLOAT := 6 Maximum value allowed
clipFun(x) == min(max(x,-clipValue),clipValue)
Notice that we identify as a small float but do not declare the type of the function clipFun. As it turns out, clipFun is called with a small float value. This declaration ensures that clipFun never does a conversion when it is called.
The second matter concerns the possible ``poles'' of a function, the actual points where the spikes have infinite values. Axiom uses normal DoubleFloat arithmetic which does not directly handle infinite values. If your function has poles, you must adjust your step size to avoid landing directly on them (Axiom calls error when asked to divide a value by , for example).
We set the variables and to hold the number of steps taken in the real and imaginary directions, respectively. Most examples will have ranges centered around the origin. To avoid a pole at the origin, the number of points is taken to be odd.
realSteps: INT := 25 Number of real steps
imagSteps: INT := 25 Number of imaginary steps
)read arrows
Now define the function drawComplexVectorField to draw the arrows. It is good practice to declare the type of the main function in the file. This one declaration is usually sufficient to ensure that other lower-level functions are compiled with the correct types.
C := Complex DoubleFloat
S := Segment DoubleFloat
drawComplexVectorField: (C -> C, S, S) -> VIEW3D
The first argument is a function mapping complex small floats into complex small floats. The second and third arguments give the range of real and imaginary values as segments like . The result is a three-dimensional viewport. Here is the full function definition:
drawComplexVectorField(f, realRange,imagRange) ==
delReal := (hi(realRange)-lo(realRange))/realSteps The real step size
delImag := (hi(imagRange)-lo(imagRange))/imagSteps The imaginary step size
sp := createThreeSpace() Create empty space $sp$
real := lo(realRange) The initial real value
for i in 1..realSteps+1 repeat Begin real iteration
imag := lo(imagRange) initial imaginary value
for j in 1..imagSteps+1 repeat Begin imaginary iteration
z := f complex(real,imag) value of $f$ at the point
arg := argument z direction of the arrow
len := clipFun sqrt norm z length of the arrow
p1 := point [real, imag, 0.0@DFLOAT, arg] base point of the arrow
scaleLen := delReal * len scaled length of the arrow
p2 := point [p1.1 + scaleLen*cos(arg), tip point of the arrow
p1.2 + scaleLen*sin(arg),0.0@DFLOAT, arg]
arrow := makeArrow(p1, p2) Create the arrow
for a in arrow repeat curve(sp, a) Add arrow to space $sp$
imag := imag + delImag The next imaginary value
real := real + delReal The next real value
makeViewport3D(sp, "Complex Vector Field") Draw it
As a first example, let us draw . There is no need to create a user function: just pass the sinsinComplex DoubleFloat from Complex DoubleFloat.
Read the file.
Draw the complex vector field of .
10.8 Drawing Complex Functions
Here is another way to graph a complex function of complex arguments. For each complex value , compute , again expressing the value in polar coordinates . We draw the complex valued function, again considering the -plane as the complex plane, using as the height (or -coordinate) and as the color. This is a standard plot---we learned how to do this in Chapter ugGraph --- but here we write a new program to illustrate the creation of polygon meshes, or grids.
Call this function drawComplex. It displays the points using the ``mesh'' of points. The function definition is in three parts.
drawComplex: (C -> C, S, S) -> VIEW3D
drawComplex(f, realRange, imagRange) == The first part
delReal := (hi(realRange)-lo(realRange))/realSteps The real step size
delImag := (hi(imagRange)-lo(imagRange))/imagSteps The imaginary step size
Initial list of list of points $llp$
llp:List List Point DFLOAT := []
Variables and give the step sizes along the real and imaginary directions as computed by the values of the global variables and . The mesh is represented by a list of lists of points , initially empty. Now alone is ambiguous, so to set this initial value you have to tell Axiom what type of empty list it is. Next comes the loop which builds .
real := lo(realRange) The initial real value
for i in 1..realSteps+1 repeat Begin real iteration
imag := lo(imagRange) initial imaginary value
lp := []$(List Point DFLOAT) initial list of points $lp$
for j in 1..imagSteps+1 repeat Begin imaginary iteration
z := f complex(real,imag) value of $f$ at the point
pt := point [real,imag, clipFun sqrt norm z, Create a point
argument z]
lp := cons(pt,lp) Add the point to $lp$
imag := imag + delImag The next imaginary value
real := real + delReal The next real value
llp := cons(lp, llp) Add $lp$ to $llp$
The code consists of both an inner and outer loop. Each pass through the inner loop adds one list of points to the list of lists of points . The elements of are collected in reverse order.
makeViewport3D(mesh(llp), "Complex Function") Create a mesh and display
The operation mesh then creates an object of type ThreeSpace(DoubleFloat) from the list of lists of points. This is then passed to makeViewport3D to display the image.
Now add this function directly to your vectors.input file and re-read the file using read vectors. We try drawComplex using a user-defined function .
Read the file.
This one has a pole at .
Draw it with an odd number of steps to avoid the pole.
10.9 Functions Producing Functions
In ugUserMake , you learned how to use the operation function to create a function from symbolic formulas. Here we introduce a similar operation which not only creates functions, but functions from functions.
The facility we need is provided by the package MakeUnaryCompiledFunction(E,S,T). MakeUnaryCompiledFunction This package produces a unary (one-argument) compiled function from some symbolic data generated by a previous computation. MakeBinaryCompiledFunction is available for binary functions. MakeBinaryCompiledFunction The tells where the symbolic data comes from; the and give Axiom the source and target type of the function, respectively. The compiled function produced has type -> . To produce a compiled function with definition , call from this package. The function you get has no name. You must to assign the function to the variable to give it that name.
Do some computation.
Convert this to an anonymous function of . Assign it to the variable to give the function a name.
Apply the function.
For a more sophisticated application, read on.
10.10 Automatic Newton Iteration Formulas
This setting is needed to get Newton's iterations to converge.
We resume our continuing saga of arrows and complex functions. Suppose we want to investigate the behavior of Newton's iteration function Newton iteration in the complex plane. Given a function , we want to find the complex values such that .
The first step is to produce a Newton iteration formula for a given : We represent this formula by a function that performs the computation on the right-hand side, that is, .
The type Expression Integer (abbreviated EXPR INT) is used to represent general symbolic expressions in Axiom. Expression To make our facility as general as possible, we assume has this type. Given , we want to produce a Newton iteration function which, given a complex point , delivers the next Newton iteration point .
This time we write an input file called newton.input. We need to import MakeUnaryCompiledFunction (discussed in the last section), call it with appropriate types, and then define the function which references it. Here is the function :
C := Complex DoubleFloat The complex numbers
complexFunPack:=MakeUnaryCompiledFunction(EXPR INT,C,C)
Package for making functions
newtonStep(f) == Newton's iteration function
fun := complexNumericFunction f Function for $f$
deriv := complexDerivativeFunction(f,1) Function for $f'$
(x:C):C +-> Return the iterator function
x - fun(x)/deriv(x)
complexNumericFunction f == Turn an expression $f$ into a
v := theVariableIn f function
compiledFunction(f, v)$complexFunPack
complexDerivativeFunction(f,n) == Create an nth derivative
v := theVariableIn f function
df := D(f,v,n)
compiledFunction(df, v)$complexFunPack
theVariableIn f == Returns the variable in $f$
vl := variables f The list of variables
nv := # vl The number of variables
nv > 1 => error "Expression is not univariate."
nv = 0 => 'x Return a dummy variable
first vl
Do you see what is going on here? A formula is passed into the function newtonStep. First, the function turns into a compiled program mapping complex numbers into complex numbers. Next, it does the same thing for the derivative of . Finally, it returns a function which computes a single step of Newton's iteration.
The function complexNumericFunction extracts the variable from the expression and then turns into a function which maps complex numbers into complex numbers. The function complexDerivativeFunction does the same thing for the derivative of . The function theVariableIn extracts the variable from the expression , calling the function error if has more than one variable. It returns the dummy variable if has no variables.
Let's now apply newtonStep to the formula for computing cube roots of two.
Read the input file with the definitions.
The cube root of two.
Get Newton's iteration formula.
Let denote the result of applying Newton's iteration once to the complex number .
Now apply it repeatedly. How fast does it converge?
Check the accuracy of the last iterate.
In MappingPackage1, we show how functions can be manipulated as objects in Axiom. A useful operation to consider here is , which means composition. For example causes the Newton iteration formula to be applied twice. Correspondingly, means to apply the iteration formula times.
Apply twice to the point .
Apply 11 times.
Look now at the vector field and surface generated after two steps of Newton's formula for the cube root of two. The poles in these pictures represent bad starting values, and the flat areas are the regions of convergence to the three roots.
The vector field.
The surface.
